Miller-Rabin Primality Testing

APL-related discussions - a stream of APL consciousness.
Not sure where to start a discussion ? Here's the place to be
Forum rules
This forum is for discussing APL-related issues. If you think that the subject is off-topic, then the Chat forum is probably a better place for your thoughts !

Miller-Rabin Primality Testing

Postby Roger|Dyalog on Fri Aug 21, 2020 3:02 am

⍺ 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.
Roger|Dyalog
 
Posts: 207
Joined: Thu Jul 28, 2011 10:53 am

Return to APL Chat

Who is online

Users browsing this forum: No registered users and 1 guest