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
Bookmarks