restart: with(numtheory); printlevel:=0: prmstart:=29: expstart:=200: expend:=300: p:=nextprime(prmstart): test:=false: for k from expstart to expend while evalb(test=false) do n:=p*(2^k)+1; test:=isprime(n): end do; cat(`k = `,k); cat(`n = `,n); printlevel:=1: if evalb(k>=expend) then print(`The exponent of 2 left your chosen range without finding a prime; choose new starting values.`); else cat(` n = `,p,`*2^(`,k,`) may be prime, it has `,ilog10(n),` digits. Checking: `); printlevel:=0: m:=2: while evalb(m&^((n-1)/2) mod n =1) do m:=m+1: end do; s:=2: while evalb(s&^((n-1)/p) mod n = 1) do s:=s+1: end do; printlevel:=1: cat(` `); cat(`n is prime. The base `,m,` verifies Lucas' condition for the prime 2 ;`); cat(`the base `,s,` verifies Lucas' condition for the prime `,p,` . Specifically:`); cat(` `); cat(m,`^`,(n-1)/2,` = `,m&^((n-1)/2) mod n,` mod `,n); cat(m,`^`,(n-1),` = `,m&^((n-1)) mod n,` mod `,n); cat(s,`^`,(n-1)/p,` = `,s&^((n-1)/p) mod n,` mod `,n); cat(s,`^`,(n-1),` = `,s&^((n-1)) mod n,` mod `,n); end if; printlevel:=0: