Miller-Rabin primality test

[danbaron]

I tried to implement the Miller-Rabin primality test, in Racket (Scheme).

http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test

Why did I use Racket? Because, I like it.

The test is probabilistic, not deterministic. It doesn't absolutely determine if an integer is prime. You set the probability, P, and if the test says that a particular integer is prime, then, the likelihood is at least P, that it is actually prime.

My experience with my implementation of it, is that it sometimes makes errors. So, maybe I did something wrong.

There are 109 primes less than 600, and the test said that 24 of them are not prime. However, it seems that as the candidate integer gets bigger, the procedure becomes more accurate.

Below is shown the Racket code, and some tests I did on integers in the Racket REPL (read-eval-print loop).

For each test, the output is either #t (true), or #f (false), and the time in seconds to do the calculation, is also indicated.

For the tests, I used primes which I got from,

http://en.wikipedia.org/wiki/List_of_prime_numbers .

I did 8 tests, on primes from the web page.

The first 6 are Carol primes, and the next 2 are Bell number primes.

The procedure correctly indicated that all 8 are prime.

Then, to make sure that the procedure wasn't just answering true for every large integer, I did the 8 tests again.

This time, I added the value 2, to each of the original 8 primes, and tested those integers.

For the second set of 8 tests, the procedure indicated that the first 7 were not prime, as I hoped it would.

However, for the 8th test in the second set, the procedure indicated that the integer is prime.

At that point I thought that the procedure was no good.

But, I decided to try one more thing.

I went to,

http://www.wolframalpha.com/ .

I executed the command,

"PrimeQ[359334085968622831041960188598043661065388726959079839]".

"Super-Mathematica" indicated, that the number is indeed prime.

So, I felt better.

In all of the tests, I set the probability, P, equal to, 0.999.

:oops: :x :twisted:


; code ---------------------------------------------------------------------------------------------

#lang racket

(define ran-max 4294967087)

(define (square x) (* x x))

(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp) (remainder (square (expmod base (/ exp 2) m)) m))
(else (remainder (* base (expmod base (- exp 1 ) m)) m))))

(define (power-of-two-factor x)
(define (power-of-two-internal x n)
(if (= (modulo x 2) 1)
n
(power-of-two-internal (/ x 2) (+ n 1))))
(power-of-two-internal x 0))

(define (factor-even-number x)
(let ((po2 (power-of-two-factor x)))
(values po2 (/ x (expt 2 po2)))))

(define (big-random n)
(let ((r (round (/ (* n (random ran-max)) ran-max))))
(if (= r 1) (big-random n) r)))

(define (num-repetitions p)
(define 1-p (- 1 p))
(ceiling (- (/ (log 1-p) (log 4)))))

(define (is-prime? n prob)
(cond
((< n 5) #f)
((even? n) #f)
(else
(define-values (s d) (factor-even-number (- n 1)))
(define k (num-repetitions prob))
(define (outer-loop k)
(define a (big-random (- n 2)))
(define x (expmod a d n))
(cond
((= k 0) #t)
((or (= x 1) (= x (- n 1))) (outer-loop (- k 1)))
(else (inner-loop (- s 1) x))))
(define (inner-loop i y)
(set! y (modulo (* y y) n))
(cond
((or (= y 1) (= i 0)) #f)
((= y (- n 1)) (outer-loop (- k 1)))
(else (inner-loop (- i 1) y))))
(outer-loop k))))

(define (time-is-prime? n p)
(let ((t1 (current-milliseconds)))
(let ((prime? (is-prime? n p)))
(let ((t2 (current-milliseconds)))
(let ((tt (/ (- t2 t1) 1000.0)))
(values prime? tt))))))

; REPL interactions --------------------------------------------------------------------------------

Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (time-is-prime? 68718952447 0.999)
#t
0.141
> (time-is-prime? 274876858367 0.999)
#t
0.122
> (time-is-prime? 4398042316799 0.999)
#t
0.108
> (time-is-prime? 1125899839733759 0.999)
#t
0.113
> (time-is-prime? 18014398241046527 0.999)
#t
0.097
> (time-is-prime? 1298074214633706835075030044377087 0.999)
#t
0.274
> (time-is-prime? 35742549198872617291353508656626642567 0.999)
#t
0.17
> (time-is-prime? 359334085968622831041960188598043661065388726959079837 0.999)
#t
1.488
> (time-is-prime? 68718952449 0.999)
#f
0
> (time-is-prime? 274876858369 0.999)
#f
0.001
> (time-is-prime? 4398042316801 0.999)
#f
0.123
> (time-is-prime? 1125899839733761 0.999)
#f
0.001
> (time-is-prime? 18014398241046529 0.999)
#f
0.001
> (time-is-prime? 1298074214633706835075030044377089 0.999)
#f
0.002
> (time-is-prime? 35742549198872617291353508656626642569 0.999)
#f
0.118
> (time-is-prime? 359334085968622831041960188598043661065388726959079839 0.999)
#t
3.527
>

[danbaron]

I tried a 100 digit prime, which I got from,

http://primes.utm.edu/lists/small/small.html#80 .

Then, I added 2 to it, and tried that.

It seems to work OK.


Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (time-is-prime? 2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550077 0.999)
#t
5.463
> (time-is-prime? 2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550079 0.999)
#f
0.086
>

[danbaron]

I tried a 200 digit prime, which I got from,

http://primes.utm.edu/lists/small/small2.html .

Then, I added 2 to it, and tried that.

It still seems to work OK.


Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (time-is-prime? 58021664585639791181184025950440248398226136069516938232493687505822471836536824298822733710342250697739996825938232641940670857624514103125986134050997697160127301547995788468137887651823707102007839 0.999)
#t
3.899
> (time-is-prime? 58021664585639791181184025950440248398226136069516938232493687505822471836536824298822733710342250697739996825938232641940670857624514103125986134050997697160127301547995788468137887651823707102007841 0.999)
#f
0.165
>

[danbaron]

I tried a 1050 digit prime, which I got from,

http://primes.utm.edu/primes/lists/all.txt .

It is #5706 in the list, 11^1008 + 998672782.

Then, since the last digit of the number is 3, I tried the number plus 4.

It still seems to work OK.


Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (define a (+ (expt 11 1008) 998672782))
> a
529452056448793684463171701371965165599024891268044535141171238938714197123500025362132544081526670388715712690296898616337461300706986447973191761836855880775521428076605793975641811081746047876506709541924966850486767835110575769384906730749667777016471469177953052696916660058190262098621182600138065372500625688141864409509834051913316703890710975152285824832755202082018739474065720249911588366808071879632022058226980343003483610162880483389826691629327795587270718562153925742749978926314647876232385694302214489674918510983781514194149348746847611591937897776475529798214185181022673653771123333515766636914861846761384172063717773371541993327403277856399026738370062644554264826051111760730911233180455686309777627985916941132711706748812406408935146311575975928657195294322425757348388741593204454372595463565373254895887022065283560617002190260565624775467971401594792170355772678679741987491858658211032034438449802653787330079926121284263030291537500873819738703978995623111675987883052680605464876084482213424458493450872074280702291663
> (define a+4 (+ a 4))
> (time-is-prime? a 0.999)
#t
15.988
> (time-is-prime? a+4 0.999)
#f
2.402
>

[Johannes]

Dan,

My program (below) agrees with all of your checks. What is striking is that my script is much faster for the first examples, performs about equally for the 200-digit numbers, and is painfullly slower for the 1050-digit numbers. At a guess Racket is using FFT multiplication for its big integers.

One question: when you calculate the number of loops, the -log(1-p)/log(4) comes out as 4.98. Is that result rounded or truncated? Because if it is truncated your program is incorrect with respect to the probability. ;)

' > BigInt_PrimeMillerRabin

Uses "Console"
Module "BigInt"
Alias String As BigInt

Randomize

TimeMR("68718952447")
TimeMR("274876858367")
TimeMR("4398042316799")
TimeMR("1125899839733759")
TimeMR("18014398241046527")
TimeMR("1298074214633706835075030044377087")
TimeMR("35742549198872617291353508656626642567")
TimeMR("359334085968622831041960188598043661065388726959079837")

TimeMR("68718952449")
TimeMR("274876858369")
TimeMR("4398042316801")
TimeMR("1125899839733761")
TimeMR("18014398241046529")
TimeMR("1298074214633706835075030044377089")
TimeMR("35742549198872617291353508656626642569")
TimeMR("359334085968622831041960188598043661065388726959079839")

TimeMR("2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550077")
TimeMR("2074722246773485207821695222107608587480996474721117292752992589912196684750549658310084416732550079")

TimeMR("58021664585639791181184025950440248398226136069516938232493687505822471836536824298822733710342250697739996825938232641940670857624514103125986134050997697160127301547995788468137887651823707102007839")
TimeMR("58021664585639791181184025950440248398226136069516938232493687505822471836536824298822733710342250697739996825938232641940670857624514103125986134050997697160127301547995788468137887651823707102007841")

TimeMR("529452056448793684463171701371965165599024891268044535141171238938714197123500025362132544081526670388715712690296898616337461300706986447973191761836855880775521428076605793975641811081746047876506709541924966850486767835110575769384906730749667777016471469177953052696916660058190262098621182600138065372500625688141864409509834051913316703890710975152285824832755202082018739474065720249911588366808071879632022058226980343003483610162880483389826691629327795587270718562153925742749978926314647876232385694302214489674918510983781514194149348746847611591937897776475529798214185181022673653771123333515766636914861846761384172063717773371541993327403277856399026738370062644554264826051111760730911233180455686309777627985916941132711706748812406408935146311575975928657195294322425757348388741593204454372595463565373254895887022065283560617002190260565624775467971401594792170355772678679741987491858658211032034438449802653787330079926121284263030291537500873819738703978995623111675987883052680605464876084482213424458493450872074280702291663")
TimeMR("529452056448793684463171701371965165599024891268044535141171238938714197123500025362132544081526670388715712690296898616337461300706986447973191761836855880775521428076605793975641811081746047876506709541924966850486767835110575769384906730749667777016471469177953052696916660058190262098621182600138065372500625688141864409509834051913316703890710975152285824832755202082018739474065720249911588366808071879632022058226980343003483610162880483389826691629327795587270718562153925742749978926314647876232385694302214489674918510983781514194149348746847611591937897776475529798214185181022673653771123333515766636914861846761384172063717773371541993327403277856399026738370062644554264826051111760730911233180455686309777627985916941132711706748812406408935146311575975928657195294322425757348388741593204454372595463565373254895887022065283560617002190260565624775467971401594792170355772678679741987491858658211032034438449802653787330079926121284263030291537500873819738703978995623111675987883052680605464876084482213424458493450872074280702291667")

WaitKey

Sub TimeMR(n As BigInt,Optional p As Ext)
Local f As Boolean
Local t As Quad
n = BigInt_FromString(n)
If p=0 Then p=.999
HiResTimer_Init
HiResTimer_Get
f = MillerRabin(n,p)
t=HiResTimer_Delta
PrintL BigInt_ToString(n)
If f Then
Print "PRIME "
Else
Print " "
EndIf
PrintL Format$(t/1000000,"0.000")
End Sub

Function MillerRabin(n As BigInt,p As Ext) As Boolean
Local a,d,m,n1,n2,t,w,x As BigInt
Local i,j,s As DWord
Local b As Byte
Local q As Ext = 1
n1 = BigInt_Dec(n)
n2 = BigInt_Dec(n1)
d = n1
While BigInt_IfEven(d)
s += 1
d = BigInt_Div(d,2)
Wend
w = BigInt_FromInteger(1)
t = BigInt_FromInteger(2)
m = BigInt_FromInteger(1000000000)
p = 1-p
Do
q /= 4
a = BigInt_Div(BigInt_Mul(n2,BigInt_FromInteger(Rnd(1,1000000000))),m)
a = BigInt_Max(a,t)
x = w
For i=2 To Len(d)
b = Asc(Mid$(d,i,1))
For j=1 To 8
If (b And 1)=1 Then x = BigInt_Mod(BigInt_Mul(x,a),n)
a = BigInt_Mod(BigInt_Mul(a,a),n)
SHIFT RIGHT b,1
Next j
Next i
If BigInt_IfEQ(x,w) Then Iterate Do
If BigInt_IfEQ(x,n1) Then Iterate Do
If s>1 Then
For i=2 To s
x = BigInt_Mod(BigInt_Mul(x,x),n)
If BigInt_IfEQ(x,w) Then Return FALSE
If BigInt_IfEQ(x,n1) Then Iterate Do
Next i
EndIf
Return FALSE
Loop Until q<=p
Return TRUE
End Function

[Johannes]

Dan, remember this one?

Factor the following composite integer.

21130819027319923178509622632727980622577522491335986162202750512019525136688895741633235387039306968582318693021114825649130164135868576427684517893

Factoring is still "impossible" but with the Miller-Rabin test it takes about half a second on my system to determine that that number is, in fact, composite.

[danbaron]

(I don't think anyone will call you lazy, Johannes.:rolleyes:)

"ceiling", so, hopefully, 4.98 --> 5.

(But, it should do at least one repetition no matter what the value of p is. As it is currently written, if you enter p as 0, then it will do no repetitions.)


(define (num-repetitions p)
(define 1-p (- 1 p))
(ceiling (- (/ (log 1-p) (log 4)))))

I remember it, so, I'll do it (or, more accurately, Racket will do it).


Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (time-is-prime? 21130819027319923178509622632727980622577522491335986162202750512019525136688895741633235387039306968582318693021114825649130164135868576427684517893 0.999)
#f
0.133
>

I ran two more numbers from the list,

http://primes.utm.edu/primes/lists/all.txt .

#5680, 546! - 1, 1260 digits

#5647, 83 * 2^5318 - 1, 1603 digits


Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (define a (- (fact 546) 1))
> (time-is-prime? a 0.999)
#t
29.135
> a
141302009261418325453762621393526286698658617796079910071353871488907640504160743565501241236107733762080837422133568434484213467457560621538131284091550631968816592748754460044234605778793037038634578162615257746896342308546754776364849349150715175577496877924111556740124186753117294177162303783405503829723023736326990734833969973121886593758921817486418079092325492863765948445727040946819291138129804200523981866521636537165676001267152151758174473660767772438139761519901508251733325317410707450984627544775891183776230678501637671689726049915871033588824272610663202323719389931283164345628668801970555064368945390634202638535930076803929062184736870310291943624402875877497305289749392357778468456405062145336924250773407066015135454416984550131804906976797622941301057536164271657477854025104530247400241587703885280095204104624235632766825184645847127845152855411530623429061372120607469651273815795437429759478958062539197279629575104299984571584633649528868101236905793259525917436877798780909275667980803203926354424450148098596386681966959342498581232330994696677186100809764161773397785695845695608487026255462399999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999
> (define b (+ a 2))
> (time-is-prime? b 0.999)
#f
3.927
> (define c (- (* 83 (expt 2 5318)) 1))
> (time-is-prime? c 0.999)
#t
66.261
> c
6260298383829657225286112584739828350961566071813614817090655661726702203000727443155352137677415727639270005152424260690857553900368341450346970741308511417934442673066866644403367228220364439878425897186245253177691894131530314490537670719045204854370814516293642618537389174767012520242616529405673247696080934123952725656223766315340714614855806548932055139425179151479881206038701763802779138916668883306784294260180704008909313195142729581852896420768539776228014321360869785451249902588785914955877269128146778963486558120668354631060167501141743903122509454201202547688715196727040523077806030531988397663538853730963809344267934674614590712227267805424158243982282344265009362357321981931239055386336840472883326309867734816135908681990888503697042094651571968230153780164618679215781152726244096578224696053246041697606541689435068727931476353894324623910695079566421040585311985103805812421650786510353427759898861776000121326090856058979212077591696177260616037165555751507285878767478858118625647588877281294325434409067338807313030610213165442239480105674376617197846030146545262330219853080619480300992036754949135625216566071983425953282414885683202259216717086223201216107231214643841453821800479705443657770734535275127870398933661350313638594922371398078999763881759388037191499672886496845677727495900732005269358724467804503279937265707718167486322051615118497652040645286382390708500068385222336755432189454238226603697941580242718059483611416057744941499700978593848565252669806462644516538399273148998665777617816962249185512503379530702903585164261671957932841026458847073533951
> (define d (+ c 2))
> (time-is-prime? d 0.999)
#f
4.931
>

When I tried numbers bigger than #5647, something went wrong, because Racket calculated forever.

We probably agree that this way is much much better than the deterministic way.

I agree, Racket probably uses FFT to multiply integers. It is amazingly fast.

:eek:

[Johannes]

"ceiling", so, hopefully, 4.98 --> 5.

Oh snap! :D

Should have seen that.

I was wondering what you were doing there and then I thought to check Wikipedia to see what the reliability was. Amazing that with only 5 random tries you can get less than a one-in-a-thousand error. Ten tries and it's less than one in a million.

I did the repetitive division for the probability because it is probably (ahem) faster than the two LOGs. But you get the prize for cryptic code on that one. :p

[Johannes]

When I tried numbers bigger than #5647, something went wrong, because Racket calculated forever.

I'm running #5551 ((2^10211-1)/306772303457009724362047724636324707614338377) with 3,030 digts now. The calculation of the number itself took almost no time at all. The Miller-Rabin test is... running.

[Johannes]

The Miller-Rabin test is... running.
Two thousand seven hundred and nineteen seconds. And the number is prime with a certainty of at least 99.9%.

Dan, most likely Racket is running into memory problems with the larger numbers. I think I remember you having to define the amount of memory you want to work with.

[danbaron]

You made a mistake, it's #5552.

http://primes.utm.edu/primes/lists/all.txt ,

#5552, (2^10211-1)/306772303457009724362047724636324707614338377, 3030 digits.

I got it to work.

Apparently, it runs OK if I don't switch to another application while it is doing the calculation.

(I have to hurry to catch the train now.)

:bom:


Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (define a (/ (- (expt 2 10211) 1 ) 306772303457009724362047724636324707614338377))
> a
214027499482579132858887371940194720410174682692622338138373623371045292557815520045653509057682304478636294725744974976651789238833093921482920608867250818705789405057130008142306132064285305101940344523042763475816275438082461956957398106525878173567591967013410748523711930321497448154784619301982087427578086554013257185656786503572050932921341886940131084021943282522786653523825328688654928819340230360303370235518970181520387289033802457247214947248909978175882415141727984730164903804581970647000146953735297504267417832419440377208919473243565195543421718881615334948535134367639133447911232324594688215243535951946986566871991661304273052089123602813503876435919896292793639825090607457533352599118935819769435330713282137648515196027235926175883961948381508991013196439250801269563828297872094509176352753919162909010872619413584794333321812231202624131410710447566208583425922710075602111102768713021645326912585208390378351221956400677988044957906609043803535938903440761559694878054670730994693110669486310601066121391218142222553185800515119789608586891469241549939835278334745107024756491673675553607864812841413553324370643743265431950334045684031272956974012748846724367603353083165299542015407437000268554013108211049235352998867154599061172192902450631443211480297377333143398551487132200694659548889381600557163997968481288650381905936567631275227525822541606126444648951697619975596603427395471954729846036842186118055564496770804555953478592189018724100044613736718497543153016886206067742990963208754687720902252465272901815886304613986568059025443766818288836826074056630570373719268157466465160045255275463910326044933411062421320264359276822436022832349134640148605752270225362622198182098596036208743649181322997561896587617927374339975372695608841653175424143558970361575061430153730163149808158987657772324110044792872185726114110281650288934568186511187455600975197045290162116876673028341111299058673294322244163689345516056432560629003122541290834177133909469395927712284413194375985543117058027782493055150822973385039712720636424382466068254658498514918545694852209364073734159037840243544843487710434969363092770382602844208930542830780248981183044496174381844740848649903560513474126998310143686324244817292846455474996479224175994893041768202847456735670094121013565841410980001573955813472922755138802997745066765925772580380989672231239919104007287958321284025825440751798625677362350384659913180420068782328239411776395251225032411741984214433501639075120947846010020464050390907477193424987293869650698033887742390819717032303359032709898401966208279135177949329785326787567171415583335186049210079801121276570981524462732903623371608282999828629679269050841202372148497104343119800490071432719173702610963797657960238515665613543949340138714843391144085749180340221931803107038697823516338722144619716665907125031046563646482592564279123116110346022854424083008537664787417445537128879620470522557656209209655022093331866484571760614730583954562451341070311239543902959465356555362216711
> (time-is-prime? a 0.999)
#t
246.243
> (define b (+ a 2))
> (time-is-prime? b 0.999)
#f
41.386
>

[Johannes]

You made a mistake, it's #5552.

Not when I looked at the page.

The time stamp for the current list, with "my" number being #5552 is

(Tue Apr 19 14:50:52 CDT 2011)That looks like a time for the US continent. When I looked at the list it was 11am in Paris. ;)

[danbaron]

Tue Apr 19 21:30 PDT

I thought about that on the train - that you were probably correct at the time you did it.

They don't wait long to change the list.

Now, the time stamp is,
(Tue Apr 19 20:50:52 CDT 2011).

That was approximately 2.5 hours ago.

I therefore retract my accusation - you are exonerated, and restored to your previous position. All fines, levies and fees are rescinded. All records shall be expunged, and all properties returned to your name. I hereby declare it to be so, and so shall it be. Etc. Tip top. Cheerio.

:!:

[danbaron]

http://primes.utm.edu/primes/lists/all.txt

(Tue Apr 19 20:50:52 CDT 2011)

#5393, 22971 * 2^22971 - 1, 6920 digits.


; new code --------------------------------------------------------------------------

(define (fact n)
(cond ((= n 1) 1)
(else (* n (fact (- n 1))))))

(define (last-digit n)
(- n (* 10 (truncate (/ n 10)))))

(define (prime-neighbor n)
(cond
((= 3 (last-digit n)) (- n 2))
(else (+ n 2))))

; REPL interactions -----------------------------------------------------------------

Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (define a (- (* 22971 (expt 2 22971)) 1))
> (time-is-prime? a 0.999)
#t
2376.215
> (time-is-prime? (prime-neighbor a) 0.999)
#f
176.898
>

[danbaron]

http://primes.utm.edu/primes/lists/all.txt

(Wed Apr 20 22:50:58 CDT 2011)

#5348, 10^9999 + 33603, 10000 digits.


; REPL interactions ---------------------------

Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (define a (+ (expt 10 9999) 33603))
> (time-is-prime? a 0.999)
#t
3629.146
> (time-is-prime? (prime-neighbor a) 0.999)
#f
610.178
>

[danbaron]

http://primes.utm.edu/primes/lists/all.txt

(Thu Apr 21 19:20:52 CDT 2011)

#5302, (2^42737 + 1) / 3, 12865 digits.


; REPL interactions ---------------------------

Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (define a (/ (+ (expt 2 42737) 1) 3))
> (time-is-prime? a 0.999)
#t
7187.778
> (time-is-prime? (prime-neighbor a) 0.999)
#f
1198.876
>