## 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

⍺ 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   6201 301 401 501 601      1 pco c1 1 0 1 00 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 c1 1 0 1 00 0 1 0 1      pco 1e8+0 1e52038074751 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 d4158115995448490359      MillerRabin d0`

(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 oddwitnesses←{ 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: 236
Joined: Thu Jul 28, 2011 10:53 am