Miller-Rabin Primality Testing
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 !
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 !
1 post
• Page 1 of 1
Miller-Rabin Primality Testing
⍺ 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.
(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.
Translated from the Miller-Rabin section of the J Wiki essay Primality Tests.
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: 238
- Joined: Thu Jul 28, 2011 10:53 am
1 post
• Page 1 of 1
Who is online
Users browsing this forum: No registered users and 1 guest
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group