Post Reply 
Fast prime finder??
08-03-2016, 11:50 AM (This post was last modified: 08-03-2016 08:48 PM by Namir.)
Post: #13
RE: Fast prime finder??
OK here is the VBA code that tests a set of prime finder routines. I frequently use Excel VBA first to write new code. The spreadsheet is very useful to view intermediate and final results. From there I convert code to Matlab, HP Prime, True Basic, and so on.

The spreadsheet needs the following input:

1) First number in the range of primes in cell B1.
2) Last number in the range of primes in cell B2.
3) Number of iterations used in the Solovoy algorithm in Cell B3.

Cells A1, A2, and A3 contain the text "First", "Last", and "Iters".

To run the test invoke the macro (i.e. SUB) Go which tests five different methods. I included the number of iterations used to determine the primes (which is a fixed number for the Solovoy algorithm):

Code:

Option Explicit

Sub Go()
  Range("C1:Z1000").Value = ""
  Call NextPrime([B1].Value, [B2].Value)
  Call NextPrime2([B1].Value, [B2].Value)
  Call NextPrime3([B1].Value, [B2].Value)
  Call NextPrime4([B1].Value, [B2].Value)
  Call NextPrime5([B1].Value, [B2].Value)
End Sub
  

Sub NextPrime(ByVal Floor As Long, ByVal Limit As Long)
  Dim Count As Long, I As Long, testdiv As Long
  Dim LoopCounter As Long, OuterLoopCounter As Long
  Dim to_test As Long, to_test_root As Long
  Dim was_prime As Boolean
  
  Count = 0
  OuterLoopCounter = 0
  to_test = Floor
 
  ''i' is used as an increment of either 2 or 4,
  ' depending on the previous prime.
  I = 0
  Do While to_test < Limit
    testdiv = 2
    to_test_root = Sqr(to_test)
    was_prime = False
    LoopCounter = 0
    
    Do
      OuterLoopCounter = OuterLoopCounter + 1
      LoopCounter = LoopCounter + 1
      If testdiv > to_test_root Then
         'Print "count - to_test\n<br>"
         Count = Count + 1
         Cells(Count, 3) = to_test
         Cells(Count, 4) = LoopCounter
         Cells(Count, 5) = OuterLoopCounter
         was_prime = True
         Exit Do
      End If
      If to_test Mod testdiv = 0 Then Exit Do
      testdiv = testdiv + 1
    Loop
 
    ' primes 2, 3 and 5 are excluded from this test:
    If was_prime And to_test > 5 Then
      I = IIf((to_test Mod 6) = 5, 2, 4)
      to_test = to_test + I
    Else
    ' not a prime, so the test doesn't apply.
       to_test = to_test + 1
    End If
   Loop
   
End Sub

Sub NextPrime2(ByVal Floor As Long, ByVal Limit As Long)
  Dim Sqrt As Long, Count As Long
  Dim I As Long, j As Integer
  Dim LoopCounter As Long, OuterLoopCounter As Long
  Dim Divisor As Long
  
  On Error GoTo HandleErr
  
  Count = 0
  LoopCounter = 0
  
  If Floor = 2 Then
    Count = 1
    Cells(Count, 7) = Floor
    Floor = Floor + 1
  End If
    
  If Floor Mod 2 = 0 Then
    Floor = Floor + 1
  End If

  For I = Floor To Limit Step 2
    If IsPrime(I, LoopCounter) Then
      Count = Count + 1
      Cells(Count, 7) = I
      Cells(Count, 8) = LoopCounter
      'Cells(Count, 9) = OuterLoopCounter
    End If
  Next I
  
ExitProc:

HandleErr:
  Resume Next
End Sub

Function IsPrime(Num As Long, ByRef LoopCounter As Long) As Boolean
  Dim I As Long
  
  IsPrime = False
  If Num < 2 Or (Num <> 2 And Num Mod 2 = 0) Then Exit Function
  For I = 3 To Sqr(Num) Step 2
    LoopCounter = LoopCounter + 1
    If Num Mod I = 0 Then Exit Function
  Next I
  IsPrime = True
End Function

Sub NextPrime3(ByVal Floor As Long, ByVal Limit As Long)
  Dim Sqrt As Long, Count As Long
  Dim I As Long, j As Integer
  Dim LoopCounter As Long, OuterLoopCounter As Long
  Dim Divisor As Long
  
  On Error GoTo HandleErr
  
  Count = 0
  LoopCounter = 0
  
  If Floor = 2 Then
    Count = 1
    Cells(Count, 10) = Floor
    Floor = Floor + 1
  End If
    
  If Floor Mod 2 = 0 Then
    Floor = Floor + 1
  End If

  
  For I = Floor To Limit Step 2
    If IsPrime2(I, LoopCounter) Then
      Count = Count + 1
      Cells(Count, 10) = I
      Cells(Count, 11) = LoopCounter
      'Cells(Count, 9) = OuterLoopCounter
    End If
  Next I
  
ExitProc:

HandleErr:
  Resume Next
End Sub

Function IsPrime2(Num As Long, ByRef LoopCounter As Long) As Boolean
  Dim I As Long
  
  IsPrime2 = False
  If Num < 4 Or Num Mod 2 = 0 Or Num Mod 3 = 0 Then Exit Function
  I = 5
  Do While I * I <= Num
    LoopCounter = LoopCounter + 1
    If Num Mod I = 0 Or Num Mod (I + 2) = 0 Then Exit Function
    I = I + 6
  Loop
  IsPrime2 = True
End Function

Sub NextPrime4(ByVal Floor As Long, ByVal Limit As Long)
  Dim Sqrt As Long, Count As Long
  Dim I As Long, j As Integer
  Dim LoopCounter As Long, OuterLoopCounter As Long
  Dim Divisor As Long
  
  On Error GoTo HandleErr
  
  Count = 0
  LoopCounter = 0
  
  If Floor = 2 Then
    Count = 1
    Cells(Count, 13) = Floor
    Floor = Floor + 1
  End If
    
  If Floor Mod 2 = 0 Then
    Floor = Floor + 1
  End If
    
  For I = Floor To Limit Step 2
    If Solovoy(I, [B3].Value, LoopCounter) Then
      Count = Count + 1
      Cells(Count, 13) = I
      Cells(Count, 14) = LoopCounter
      'Cells(Count, 9) = OuterLoopCounter
    End If
  Next I
  
ExitProc:

HandleErr:
  Resume Next
End Sub


' Next three routines are transalted from C++ code found in:
' http://www.sanfoundry.com/cpp-program-implement-solovay-strassen-primality-test/
'
Function Modulo(ByVal base As Long, ByVal exponent As Long, ByVal ModVal As Long) As Long
  Dim x As Long, y As Long
  
  x = 1
  y = base
  Do While exponent > 0
    If exponent Mod 2 = 1 Then x = (x * y) Mod ModVal
    y = (y * y) Mod ModVal
    exponent = exponent \ 2
  Loop
  Modulo = x Mod ModVal
End Function

Function calculateJacobian(ByVal a As Long, ByVal n As Long) As Integer
  Dim ans As Integer
  Dim temp As Long
  
  If a = 0 Then
    calculateJacobian = 0
    Exit Function
  End If
  
  ans = 1
  If a < 0 Then
    a = -a
    If n Mod 4 = 3 Then ans = -ans
  End If
  
  If a = 1 Then
    calculateJacobian = ans
    Exit Function
  End If
  
  Do While a <> 0
    If a < 0 Then
      a = -a
      If n Mod 4 = 3 Then ans = -ans
    End If
    
    Do While a Mod 2 = 0
      a = a \ 2
      If n Mod 8 = 3 Or n Mod 8 = 5 Then ans = -ans
    Loop
    
    temp = a
    a = n
    n = temp
    
    If a Mod 4 = 3 And n Mod 4 = 3 Then ans = -ans
    a = a Mod n
    If a > n \ 2 Then a = a - n
  Loop
  
  If n = 1 Then
    calculateJacobian = ans
  Else
    calculateJacobian = 0
  End If
End Function

Function Solovoy(ByVal p As Long, ByVal iteration As Integer, ByRef LoopCounter As Long) As Boolean
  Dim I As Integer
  Dim a As Long, Jacobian As Integer, ModVal As Long
  Dim vPrimes As Variant
  
  vPrimes = Array(0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
  Solovoy = False
  If p < 2 Then Exit Function
  
  For I = 1 To UBound(vPrimes)
    If p = vPrimes(I) Then
      Solovoy = True
      Exit Function
    End If
    
    If p > vPrimes(I) And p Mod vPrimes(I) = 0 Then
      Solovoy = False
      Exit Function
    End If
  Next I
  
  Solovoy = False
  Randomize Timer
  LoopCounter = 0
  For I = 1 To iteration
    LoopCounter = LoopCounter + 1
    a = 1 + CLng(Rnd(1) * 2 ^ 31) Mod (p - 1)
    Jacobian = (p + calculateJacobian(a, p)) Mod p
    ModVal = Modulo(a, (p - 1) \ 2, p)
    If Jacobian = 0 Or ModVal <> Jacobian Then Exit Function
  Next I
  Solovoy = True
End Function

Sub NextPrime5(ByVal Floor As Long, ByVal Limit As Long)
  Dim Sqrt As Long, Count As Long
  Dim I As Long, j As Integer
  Dim LoopCounter As Long, OuterLoopCounter As Long
  Dim Divisor As Long
  
  On Error GoTo HandleErr
  
  Count = 0
  LoopCounter = 0
  
  If Floor = 2 Then
    Count = 1
    Cells(Count, 16) = Floor
    Floor = Floor + 1
  End If
    
  If Floor Mod 2 = 0 Then
    Floor = Floor + 1
  End If

  
  For I = Floor To Limit Step 2
    If IsPrime3(I, LoopCounter) Then
      Count = Count + 1
      Cells(Count, 16) = I
      Cells(Count, 17) = LoopCounter
      'Cells(Count, 9) = OuterLoopCounter
    End If
  Next I
  
ExitProc:

HandleErr:
  Resume Next
End Sub

Function IsPrime3(ByVal Number, ByRef LoopCounter As Long) As Boolean
' Optimised isPrime function
' from
' http://www.vbforums.com/showthread.php?104395-Prime-number-generation-rapid-algorithm-required
'
 Dim check As Long, SR As Double, j As Long
 
 If Number = 2 Or Number = 3 Or Number = 5 Or Number = 7 Or Number = 11 Or Number = 13 Then
    ' if 2,3,5,7,11,13 then must be a prime
    IsPrime3 = True
    Exit Function
 End If
 
 If Number Mod 2 = 0 Or Number Mod 3 = 0 Or Number Mod 5 = 0 Or Number Mod 7 = 0 _
    Or Number Mod 11 = 0 Or Number Mod 13 = 0 Then
    ' if not 2,3,5,7,11,13, yet divisable by these nums then can't be prime
    Exit Function
 End If
    
 SR = Sqr(Number)
 check = 15
 
 ' Now instead of checking the mod of each integer less than a number to see if it is
 ' prime, we can skip checks for numbers divisable by 2,3 and 5 meaning the algorithm is
 ' mucho optimised.  Also the alogorithm terminates when the checking integer increases
 ' passed the square root of number, which further optimises the algo.
 Do
    For j = 1 To 8
        LoopCounter = LoopCounter + 1
        
        If j = 1 Then
            check = check + 2
        ElseIf j = 2 Then check = check + 2
        ElseIf j = 3 Then check = check + 4
        ElseIf j = 4 Then check = check + 6
        ElseIf j = 5 Then check = check + 2
        ElseIf j = 6 Then check = check + 6
        ElseIf j = 7 Then check = check + 4
        ElseIf j = 8 Then check = check + 2
        End If
                        
                       
        ' If check > SR, as you can't have a factor greater than a number's
        ' square root, it must be prime.
        If check > SR Then
            IsPrime3 = True
            Exit Function
        End If
            
        ' This number has a factor, so can't be prime
        If Number Mod check = 0 Then Exit Function
            
    Next j
    check = check + 2
     
 Loop
 
End Function

In the case of the Solovoy function, I added the following VBA code to the translated C++ code to filter out composite numbers that were coming out as primes:

Code:

 vPrimes = Array(0, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
  Solovoy = False
  If p < 2 Then Exit Function
  
  For I = 1 To UBound(vPrimes)
    If p = vPrimes(I) Then
      Solovoy = True
      Exit Function
    End If
    
    If p > vPrimes(I) And p Mod vPrimes(I) = 0 Then
      Solovoy = False
      Exit Function
    End If
  Next I

The prime finder results imrpove from NextPrime to NextPrim1 ... to NextPrime4 and NextPrime5.

I included VBA Excel code since Excel is very popular, allowing folks who read this message to easily copy and paste the VBA code, setup the spreadsheet, and run the macro Go.

Namir
Find all posts by this user
Quote this message in a reply
Post Reply 


Messages In This Thread
Fast prime finder?? - Namir - 07-28-2016, 08:22 PM
RE: Fast prime finder?? - Claudio L. - 07-29-2016, 02:50 AM
RE: Fast prime finder?? - Namir - 07-29-2016, 07:27 AM
RE: Fast prime finder?? - Namir - 07-29-2016, 12:03 PM
RE: Fast prime finder?? - ttw - 08-01-2016, 09:59 PM
RE: Fast prime finder?? - Namir - 08-02-2016, 12:01 AM
RE: Fast prime finder?? - ttw - 08-03-2016, 08:36 AM
RE: Fast prime finder?? - Namir - 08-03-2016, 11:07 AM
RE: Fast prime finder?? - Namir - 08-03-2016 11:50 AM
RE: Fast prime finder?? - Namir - 08-04-2016, 11:54 PM
RE: Fast prime finder?? - Claudio L. - 08-05-2016, 12:46 AM



User(s) browsing this thread: 1 Guest(s)