Prime number tests

A simple pseudo prime number test is using Fermat's little theorem: 

for any prime number p and integer a such that 0<a<p. To select a first there is a need for a pseudo random number generator.

variable rnd base @ hex 

: reset_seed  0ABCDEF1 rnd ! ; reset_seed 

: rand ( -- u )
  rnd @ 08088405 * 1+ dup rnd ! ;

: random ( u1 -- u2 ) 
  rand um* nip ;

base ! 

A good thing with fermat is that flag always is true if p is a prime.

: fermat ( numb -- flag ) 
  dup 2 - random 2 + 
  over 1- rot u**mod 1 = ; 

It could therefore be used as a rather fast first primality test to see whether it's necessary to take a closer look.

Rabin-Miller strong pseudoprime test

is an efficient probabilistic algorithm for determining if a given number is prime.

: 2/mod \ n -- r q
  dup 1 and swap 1 rshift ;

: get-rs 0 loc{ numb r -- r s }
  numb 1-
  begin dup 2/mod swap 0=            \ n qout rest=0
  while nip r 1+ to r
  repeat 2drop 
  r numb 1- r rshift ;

: get-a ( numb -- a )
  2 - random 2 + ;

: rabin-miller1 ( numb -- flag )
  dup dup get-rs rot get-a false loc{ numb  r s a flag }
  a s numb u**mod 1 = 
  if true exit
  then false to flag r 0 
  ?do a s i lshift numb u**mod numb 1- = 
     if true to flag leave then
  loop flag ;

Both fermat and rabin-miller1 depends on pseudo random numbers that varies from test to test and are not even combined accurate enough. A repeated test with all odd numbers between 3 and 1000000 gave the series 12 12 7 17 13 of number of errors. A second test series with fermat plus 2 times rabin-miller gave the series 1 1 2 1 3.

Due to OEIS A014233 it is enough to repeat the Rabin-Miller test a few times, depending of the word length, to be sure. Instead of random numbers one use the first prime numbers: 

2 3 5 7 11 13 17 19 23 29 31 37 ...

Repeated Rabin-Miller will reveal all composite numbers less than:

         2027 if a=2
      1373653 if a=2,3
     25326001 if a=2,3,5

   3215031751 if a=2,3,5,7

2152302898747 if a=2,3,5,7,11

Since 2^32 < 2152302898747 a=2,3,5,7,11 suffice for 32 bit integers. For an accurate test of all 64 bit integers it's enough to test for a=2,3,5,7,11,13,17,19,23,29,31,37.

create pnr 
2 c, 3 c, 5 c, 7 c, 11 c, 13 c, 17 c, 19 c, 23 c, 29 c, 31 c, 37 c,

The word create defines the word pnr which leaves an address on tos, the address where the bytes 2 to 37 are allocated.

: rabin-miller2 loc{ numb a | s r flag -- flag } \ odd numb>a>1
  numb get-rs to s to r
  a s numb u**mod 1 =
  if true exit then
  false to flag r 0
  ?do a s i lshift numb u**mod numb 1- =
      if true to flag leave then
  loop flag ;

: rmx ( numb -- n )   \ 32 bit integers
  dup     2047 u< if drop 1 exit then
  dup  1373653 u< if drop 2 exit then
  dup 25326001 u< if drop 3 exit then
    3215031751 u< if 4 else 5 then ;

: rmx ( numb -- n )   \ 64 bit integers
  dup            2047 u< if drop 1 exit then
  dup         1373653 u< if drop 2 exit then
  dup        25326001 u< if drop 3 exit then
  dup      3215031751 u< if drop 4 exit then
  dup   2152302898747 u< if drop 5 exit then
  dup   3474749660383 u< if drop 6 exit then
  dup 341550071728321 u< if drop 8 exit then
  3825123056546413051 u< if 11 else 12 then ;

: rabin-miller ( numb -- flag )
  dup rmx 0
  do dup i pnr + c@ rabin-miller2 0=
     if 0= leave then
  loop 0= 0= ;

Or faster when composite:

: fermat-rabin-miller ( numb -- flag )   \ numb odd
  dup fermat
  if rabin-miller
  else 0=
  then ;

suitable to define 

: nextprime ( numb -- prime )  \ numb unsigned integer
  dup 3 u< if drop 2 exit then
  1 or
  begin dup fermat-rabin-miller 0=
  while 2 +
  repeat ;

: prime ( numb -- flag )
  dup 3 u< if 2 = exit then
  dup 1 and 0= 
  if drop false exit 
  then rabin-miller ;