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.
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:
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.