## Monday, October 26, 2015

### 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 ;