Miller-Rabin Primality Testing

[Roger|Dyalog]

⍺ MillerRabin ⍵ tests for primality on rational number ⍵. When ⍵<3.4155e14 or if the result is 0 (composite), the function is deterministic; when ⍵≥3.4155e14 and the result is 1 (prime), the function is probabilistic with error probability ≤ 0.25*⍺. No factors are provided.

First, some examples. We use function pco from the dfns workspace, in particular 1 pco ⍵ tests for primality for 32-bit integer ⍵, and pco ⍵ gives the ⍵-th prime.

      )copy dfns pco
      ⊢ c←0 1+⍤0 1⊢1 100∘.×2 3 4 5 6 
  2   3   4   5   6
201 301 401 501 601
      1 pco c
1 1 0 1 0
0 0 1 0 1

      qi c    ⍝ convert ordinary numbers to rational numbers
┌───────────┬───────────┬───────────┬───────────┬───────────┐
│┌─┬─┬─┐    │┌─┬─┬─┐    │┌─┬─┬─┐    │┌─┬─┬─┐    │┌─┬─┬─┐    │
││2│1│1│    ││3│1│1│    ││4│1│1│    ││5│1│1│    ││6│1│1│    │
│└─┴─┴─┘    │└─┴─┴─┘    │└─┴─┴─┘    │└─┴─┴─┘    │└─┴─┴─┘    │
├───────────┼───────────┼───────────┼───────────┼───────────┤
│┌─────┬─┬─┐│┌─────┬─┬─┐│┌─────┬─┬─┐│┌─────┬─┬─┐│┌─────┬─┬─┐│
││2 0 1│1│1│││3 0 1│1│1│││4 0 1│1│1│││5 0 1│1│1│││6 0 1│1│1││
│└─────┴─┴─┘│└─────┴─┴─┘│└─────┴─┴─┘│└─────┴─┴─┘│└─────┴─┴─┘│
└───────────┴───────────┴───────────┴───────────┴───────────┘
      MillerRabin⍤0 qi c
1 1 0 1 0
0 0 1 0 1

      pco 1e8+0 1e5
2038074751 2040217609

      ⊢ d←×Q/ qi pco 1e8+0 1e5
┌───────────────────────────────────────────┐
│┌─────────────────────────────────────┬─┬─┐│
││4 1 5 8 1 1 5 9 9 5 4 4 8 4 9 0 3 5 9│1│1││
│└─────────────────────────────────────┴─┴─┘│
└───────────────────────────────────────────┘
      ⍕Q d
4158115995448490359
      MillerRabin d
0

(Misalignments in the display are due to defects in the APL Chat Forum software.)

The number 19-digit number d would be difficult to test by conventional means such as trial division, because it only has two large prime factors, neither of which is close to the square root of d.

The Code

Function qi and monadic operator Q are from the APL Chat Forum post Rational Arithmetic. qi converts ordinary numbers to rational numbers. f Q works on rational numbers. Monadic operator qmodpower is defined and described in the APL Chat Forum post modpower.

The magic numbers in the witnesses function are from http://primes.utm.edu/prove/prove2_3.html.

      ⍝
assert←{⍺←'assertion failure' ⋄ 0∊⍵:⍺ ⎕SIGNAL 8 ⋄ shy←0}

huo← {1=2|⊃⌽⊃⊃⌽⍵:,⍵ ⋄ ⍵,∇(¯1↑⍵)÷Q qi 2}  ⍝ halve until odd

witnesses←{
 r←qi' '(≠⊆⊢)' 9080191 4759123141 2152302898747 3474749660383 341550071728321'
 s←(31 73)(2 7 61)(2 3 5 7 11)(2 3 5 7 11 13)(2 3 5 7 11 13 17)
 (≢r)>i←(r>Q ⍵)⍳qi 1:i⊃s
 (qi 1)+Q¨?Q ⍺⍴⍵-Q qi 1
}

qmodpower←{   ⍝ ⍺⍺|⍺*⍵
 m←⍺⍺
 z←{⊃⌽⍵:m|Q×Q/⍺ ⋄ 0⌷⍺}
 ((qi 1),m|Q ⍺){⍬≡⍵:0⌷⍺ ⋄ ((⍺ z ⍵),m|Q×Q⍨1⌷⍺)∇ ¯1↓⍵}⊤Q ⍵
}

MillerRabin←{
 assert 0=≢⍴⍵:
 assert((⊂,1),1)≡1↓⊃⍵:
 ⍺←40
 0=2|⊃⌽⊃⊃⍵:⍵≡qi 2
 ws←⍺ witnesses ⍵
 (⍵<Q qi 100)∧⍵∊ws:1
 q1←qi 1
 e←huo ⍵-Q q1
 ws{0=≢⍺:1 ⋄ (∨/c=Q ⍵-Q q1)⍱q1=Q ¯1↑c←(⍬⍴⌽⍺)(⍵ qmodpower)⍤0⊢e:0 ⋄ (1↓⍺)∇ ⍵}⍵
 1
}

Translated from the Miller-Rabin section of the J Wiki essay Primality Tests.