如何实现维基百科的素数测试算法？

Question

On Wikipedia the algorithm below 应该通过概率 Rabin-Miller 素数检验来测试奇数 n 是否是合数。

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_

我在下面BigIntANSForth on GitHub 中的算法实现可能是错误的。前缀“b”代表“大”；有一个大数字的参数堆栈和一个额外的大数字堆栈'bx'。 'ys' 也是一个单元格整数的额外堆栈。

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'

bar 前缀代表“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

8356813658135713571357310573157138571385713057130111111111111111111111112429 由于任何精度的实现都是复合的，但由于Wolfram Alpha 是素数。

我不确定我是否正确解释了算法。有什么建议吗？

Answer 1

这是实现 Miller-Rabin 方法的方法。

// 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.

希望这会有所帮助。祝你好运！