# Mintoris Forum

### Author Topic: Is N Prime? Test for 1-15 digit integerss  (Read 3220 times)

#### Laszlo

• Full Member
•   • Posts: 34 ##### Is N Prime? Test for 1-15 digit integerss
« on: Feb 08, 2011, 04:59 PM »
For large integers trying all possible divisors for testing whether a number is prime, is very slow. Up to15 digit numbers the program below works faster, performing the Deterministic Miller-Rabin test, after first excluding numbers divisible with primes less than 200. (The main point is to demonstrate the algorithms. Testing with trial divisions do work for 15 digit numbers.)

The included function ModPow(x,p,m) computes x^p mod m, for all positive parameters up to 2^50-1. The result is computed exactly, even though x^p could have trillions of digits. It uses a version of modular multiplication, which is exact also for all positive parameters up to 2^50-1.

This is quite complicated code, so give it a try, and post feedback here.

' Is N Prime? Deterministic Miller-Rabin test for 1-15 digit integers (Laszlo Hars)
' - Primes: 2011 11 101 1009 10007 100003 1000003 10000019 100000007 1000000007 (not 10000049000057)

Rounding Off                            ' For full precision of double floating points:P
FixDecimal 15,1
Global P(),Two()
Dim P(45)=2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,127,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199
Dim Two(127) = 1
For j = 1 To 127
Two(j) = 2*Two(j-1)
Next j

Do
Input "Enter a number for primality test",n

If n > 1125899906842623 Then          ' 2^50-1 (we can probably allow [2,2^51-1])
Print "Error - too large number"
ElseIf n < 2 Then
Print LI\$(n), "is NOT prime."
Else
For i = 0 To 45                     ' test small primes
If n % P(i) = 0 Then
prm = (n = P(i))                ' divisible: we can answer
i = 9999                        ' break
EndIf
Next i
' no small factor: not sure yet
If i < 9999 Then
If n < 44521 Then
prm = 1                         ' less than first untested prime^2
Else
prm = MillerRabin(n)
EndIf
EndIf

If prm=0 Then Print LI\$(n), "is NOT prime."
If prm>0 Then Print LI\$(n), "is prime."
EndIf
Loop

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub MillerRabin(n)                      ' The Deterministic Miller Rabin test
If n < 1373653 Then                   ' http://oeis.org/A014233
L = 2
ElseIf n < 25326001 Then
L = 3
ElseIf n < 3215031751 Then
L = 4
ElseIf n < 2152302898747 Then
L = 5
ElseIf n < 3474749660383 Then
L = 6
ElseIf n < 341550071728321 Then
L = 7
Else ' if n < 3825123056546413051 ~ 3.8e18 (62 bits), Else D > 11, not exactly known
L = 9
EndIf

d = n-1
s = 0
Do While d&1 = 0
d = d>>1                            ' remove trailing 0 bits from n-1
s = s+1                             ' length of trailing 00‥0 bits
Loop
For i = 0 To L
x = ModPow(P(i),d,n)
If x<>1 AND x<>n-1 Then
For j = 2 To s
x = ModPow(x,2,n)
If x = 1 Then
return 0
ElseIf x = n-1 Then
j = s+2                       ' break
EndIf
Next j
If j < s+2 Then return 0          ' continue i-loop after break out (j>s+1)
EndIf
Next i
return 1
End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub ModPow(x,p,m)                       ' x^p mod m with binary right->left method
If p < 1 Then Return 1
L = Len2(m)
r = floor(Two(2*L)/m)                 ' L-bit integer reciprocal
y = 1
Do While p > 1
If p%2 Then y = ModMult(y,x,m,L,r)  ' y = y*x mod m
x = ModMult(x,x,m,L,r)
p = floor(p/2)                      ' no (32-bit) bitwise operations
Loop
Return ModMult(y,x,m,L,r)
End Sub
' Barrett's multiplication for 1‥50 bit ints
Sub ModMult(a,b,m,L,r)                  ' a*b mod m (L=Len2(m), r=2^(2L)/m)
MLmult(&abM,&abL, a,b,L+1,L+2)        ' 2 extra bits in abL
MLmult(&qM,&qL, abM,r,L-1,0)          ' qM ~ quotient, error = 0‥2
s = abL - LSmult(qM, m, L+2)          ' need 2 extra bits to correct qM error
If s < 0 Then s = s + Two(L+2)        ' s can only be a couple bits longer than m
return s % m
End Sub

Sub MLmult(M,L, a,b, m,n)       ' M: bits of a*b with m LS bits discarded, L: LS n bits
w  = Two(26)
aL = a % w
bL = b % w
aM = floor(a/w)
bM = floor(b/w)

y = aM*bL + aL*bM
x = aL*bL + (y%w)*w                   ' bit 53 could be set: carry

L = x % Two(n)
M = (aM*bM + floor(y/w)) * Two(52-m) + floor(x/Two(m))
End Sub

Sub LSmult(a,b,n)                       ' LS n bits of the product ab
w  = Two(ceil(n/2))
a0 = a % w                            ' LS part
b0 = b % w
return (((a0*floor(b/w)+floor(a/w)*b0) % w) * w + a0*b0) % Two(n)
End Sub

Sub Len2(n)                             ' #significant bits in n
return Len(Trim\$(Replace\$(bin\$(n),"0"," ")+"a")) - 1
End Sub

Sub LI\$(z)                              ' Long Integer output format
s\$ = Str\$(abs(z))                     ' ASSUME: FixDecimal 15,1 (scientific notation)
s\$ = Left\$(s\$,1)+Right\$(s\$,Len(s\$)-2) ' Remove chr2 (decimal point)
n  = Val(Right\$(s\$,3))                ' exponent: +12
If n < 0 Then return "0"
If z < 0 Then return "-" + Left\$(s\$,n+1)
return Left\$(s\$,n+1)
End Sub

Edit20110213: removed initializations for return parameters (to match Basic vers. 3.4.3)
« Last Edit: Feb 13, 2011, 04:38 PM by Laszlo »