How to implement the primality test algorithm from Wikipedia?

Question

On Wikipedia the algorithm below is supposed to test if an odd integer, n, is a composite by the probabilistic Rabin-Miller primality test.

Input: n > 3, an odd integer to be tested for primality;
Input: k, a parameter that determines the accuracy of the test
Output: _composite_ if n is composite, otherwise _probably prime_
write n − 1 as 2^r·d with d odd by factoring powers of 2 from n − 1
WitnessLoop: repeat k times:
   pick a random integer a in the range [2, n − 2]
   x ← a^d mod n
   if x = 1 or x = n − 1 then do next WitnessLoop
      repeat r − 1 times:
      x ← x^2 mod n
      if x = 1 then return _composite_
      if x = n − 1 then do next WitnessLoop
   return _composite_
return _probably prime_

My implementation of the algorithm in BigIntANSForth on GitHub below might be erroneous. The prefix 'b' stands for 'big'; there is a parameterstack for big numbers and an extra big number stack 'bx'. Also 'ys' is an extra stack for one cell integers.

barmod ( w -- v ) is a faster variant of
bmod ( w u -- v ) where u is stored in a big variable 'den' used in the word 'barmod'

The bar prefix stands for 'Barrett reduction'.

5 value accuracy
: bmiller~ \ u -- | k -- f   false=composite, u odd >3, k is accuracy.
  bdup >bx rs true >ys accuracy 0          \ d   | r k 0
  do bx xsi bover bx b**mod                \ d x | r
     bdup 1 digit= b1+ bx b= or 0=         \ d   | r f
     if dup 1- 0                           \ d   | r r-1 0
        do bx bdup b* barmod               \ d x | r
           bdup 1 digit= b1+ bx b= 0= or   \ d   | r f
           if false ys! leave
           then
        loop ys@ 0= if leave then
     then
  loop drop bdrop xdrop ys> ;
\ b**mod ( b e m -- x ) calculates x=b^e(mod m)
\ The word rs (u -- d | -- r ) above calculates d (big) and r (cell) as described in the algorithm.
\ 'bx xsi' gives the random (big) integer a in the algorithm
\ '1 digit=' compare top of big stack with the single cell number '1'
\ 'xdrop' drop the bx stack and 'bx' correspond to 'r@' etc

83568136581357135713573105731571385713857130571301111111111111111111111112429 is a composite due to the implementation for any accuracy, but a prime due to Wolfram Alpha.

I'm not sure that I have interpreted the algorithm correctly. Any suggestion?

Answer 1

Here's how you can implement Miller-Rabin method.

// This procedure is called for all k trials. (explained later)
// It returns false if n is composite and
// returns true if n is probably prime
// d is an odd number such that d * 2ʳ = n - 1 for some r >= 1

Procedure millerTest(int n, int d):
1. Pick a random number 'a' in range [2, n-2]
2. Compute: x = (aᵈ) mod n
3. If x is 1 or n-1, Return true
4. while d doesn't equal to n - 1
    a) x = (x*x) mod n
    b) if x is equal to 1 Return false.
    c) if x is equal to n-1 Return true


// The procedure returns false if n is composite and 
// returns true if n is probably prime.
// k is an input parameter that determines the accuracy level.
// Higher value of k indicates more accuracy

Procedure isPrime(int n, int k):
1. Handle base cases for n < 3
2. If n is even Return false
3. Find an odd number d such that n-1 can be written as d * 2ʳ
   Since n can't be even here, so n-1 must be even
   and r must be greater than 0.
4. Do the following for k times
   if millerTest(n, d) returns false
        Return false.
5. Return true.

Hope this helps. Good luck!