Generate Sorcerean Primes

'Generate Sorcerean Primes
'-------------------------
'Numbers whose residues are prime after being reduced modulo 2^n, for any n.
'1 is considered prime, for consistency. :)

'eg, 107 --> 43 (mod 64) --> 11 (mod 32 mod 16) --> 3 (mod 8 mod 4) --> 1 (mod 2)

_DEFINE A-Z AS _UNSIGNED _INTEGER64

PRINT "Generating powers..."
DIM pow(63)
pow(0) = 1: FOR i = 1 TO 63: pow(i) = pow(i - 1) * 2: NEXT

PRINT "Generating small primes...";
DIM pr(99999)
pr(0) = 2: pr(1) = 3: pr(2) = 5
i = 2: n = 5
WHILE i <= 99999
  GOSUB ptest
  IF f = 1 THEN pr(i) = n: i = i + 1
  IF n MOD 6 = 1 THEN n = n + 4 ELSE n = n + 2
WEND
PRINT pr(99999)

PRINT "Generating Sorcerean primes..."
PRINT 1; 3;
a = 5
DO
  b = INT(LOG(a) / LOG(2) + 1.5) 'estimate bits needed to represent n
  FOR i = 1 TO b
    n = a MOD pow(i)
    GOSUB ptest
    IF f = 0 THEN i = b
  NEXT
  IF f = 1 THEN PRINT a;
  IF a MOD 6 = 1 THEN a = a + 4 ELSE a = a + 2
LOOP

ptest:
f = 1
s = SQR(n)
j = 0: DO
  IF n MOD pr(j) = 0 THEN f = 0
  j = j + 1
LOOP WHILE pr(j) <= s AND f = 1
RETURN