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 »