Fast prime finder?? - Namir - 07-28-2016 08:22 PM
Hi All,
I am looking for the fastest algorithm to:
1) Determine if a number is a prime.
2) Find the next prime after a given number.
I plan to include the best answer in my HHC2016 presentation.
The keyword here is fast computational-wise.
Thanks!
Namir
RE: Fast prime finder?? - Claudio L. - 07-29-2016 02:50 AM
(07-28-2016 08:22 PM)Namir Wrote: Hi All,
I am looking for the fastest algorithm to:
1) Determine if a number is a prime.
2) Find the next prime after a given number.
I plan to include the best answer in my HHC2016 presentation.
The keyword here is fast computational-wise.
Thanks!
Namir
This problem has been studied a lot, and the algorithms are well known. I'm not sure what you are aiming at: calculator code? any code? or just the algorithm? what magnitude of numbers do you need to test?
I wrote a few routines for newRPL that use a bitmap table for primes <120000. This allows you to both test for primality and get the next prime number (given a number <120000) extremely fast, which speeds up the brute force primality test for numbers up to 120000^2 = 1.44e10.
For larger numbers it falls back to brute force (test the remainder after division by every prime <= sqrt(N)).
In the future it will use Miller-Rabin or similar pseudo-prime algorithm for large numbers, but that's not implemented yet.
I recently tested the number 9,999,999,967 per this thread:
http://www.hpmuseum.org/forum/thread-6567.html
Since it's < 1.44e10 I got a very fast result at 0.057s.
While the tables are encoded in a creative way, there's nothing to be proud of, as it's basically a brute force primality test accelerated with look-up tables.
RE: Fast prime finder?? - Namir - 07-29-2016 07:27 AM
I am looking for the algorithm(s).
Namir
RE: Fast prime finder?? - Namir - 07-29-2016 12:03 PM
I found what I was looking for on the Internet.
:-)
Namir
RE: Fast prime finder?? - Namir - 08-01-2016 05:24 PM
(08-01-2016 01:45 PM)Mike (Stgt) Wrote: (07-29-2016 12:03 PM)Namir Wrote: I found what I was looking for on the Internet.
:-)
Namir
Not showing what you found useful on the Internet implies you assume it could not be of general interest. Or all but you knew it already.
What makes you so nosa about prime numbers??
Ciao.....Mike
I was hoping to find a clever equation that calculates primes. I did find some code that uses a loop (iterating over odd numbers) that looks for primes.
Namir
RE: Fast prime finder?? - rprosperi - 08-01-2016 07:00 PM
(08-01-2016 05:24 PM)Namir Wrote: I was hoping to find a clever equation that calculates primes. I did find some code that uses a loop (iterating over odd numbers) that looks for primes.
I think Mike was suggesting that readers of this thread likely would also be interested in such a solution, and it would be courteous to share the link you found.
Unless that steals thunder for your HHC presentation, which I am looking forward to as always.
RE: Fast prime finder?? - ttw - 08-01-2016 09:59 PM
There are really no formulas for primes. The primes seem to occur at random (with a known density of ln(N)/N for the number of primes up to N). One can do clever sieving and use some pseudo-prime tests (like Fermat's) to check.
For searching (not on calculators), I've often used a storage of 30 numbers in 8 bits. Checking all numbers not divisible by 30, one get that primes have remainder of (1,7,11,13,17,18,23,29). I pack 30 numbers into one byte and set the appropriate bit if that number is prime. One can go higher with more complicated packing. This is still linear. One can do a bit better.
A more efficient packing is to store the difference between primes as primes get large. This is a bit complicated.
None of this tells whether a number is prime, but most prime checking methods work better with a bunch of stored primes.
https://en.wikipedia.org/wiki/Primality_test gives some methods.
RE: Fast prime finder?? - Namir - 08-02-2016 12:01 AM
(08-01-2016 09:59 PM)ttw Wrote: There are really no formulas for primes. The primes seem to occur at random (with a known density of ln(N)/N for the number of primes up to N). One can do clever sieving and use some pseudo-prime tests (like Fermat's) to check.
For searching (not on calculators), I've often used a storage of 30 numbers in 8 bits. Checking all numbers not divisible by 30, one get that primes have remainder of (1,7,11,13,17,18,23,29). I pack 30 numbers into one byte and set the appropriate bit if that number is prime. One can go higher with more complicated packing. This is still linear. One can do a bit better.
A more efficient packing is to store the difference between primes as primes get large. This is a bit complicated.
None of this tells whether a number is prime, but most prime checking methods work better with a bunch of stored primes.
https://en.wikipedia.org/wiki/Primality_test gives some methods.
Thanks ttw for the link. That Wikipedia page has p-code that seems to be a very good prime tester, compared to the few other similar functions I have found so far.
Namir
RE: Fast prime finder?? - Namir - 08-02-2016 08:28 PM
(08-02-2016 03:34 PM)Mike (Stgt) Wrote: I have here a paper titeled "Die Methode der häufigen Zeugen", with the subtitle "Der randomisierte Primzahltest von Solovay und Strassen". I got it from a very nice lad - for me only, not for distribution. Bing for 'The Solovay-Strassen Algorithm' and you'll find some hits about it.
Ciao.....Mike
Thanks for the tip. I am looking for the topic AND source code implementation to see how coders deal with raising numbers to powers of even larger numbers!!
:-)
Namir
RE: Fast prime finder?? - EdS2 - 08-03-2016 06:30 AM
(08-02-2016 08:28 PM)Namir Wrote: I am looking for the topic AND source code implementation to see how coders deal with raising numbers to powers of even larger numbers!!
It's not just exponentiation, which means lots of digits, it's modular exponentiation, which means you don't need to keep all the digits. Rosetta code can be a good resource:
https://rosettacode.org/wiki/Modular_exponentiation
RE: Fast prime finder?? - ttw - 08-03-2016 08:36 AM
For exponentiation of large numbers, a reasonable (if not optimal) strategy is to use binary exponentiation followed by a modular reduction at each step. Intermediates can be held to no more that twice the number of digits as the numbers of interest. Optimal orders of multiplication (factorization of the exponent or a tree structure, see Knuth vol. 2) are only useful if the same exponent is used repeatedly.
RE: Fast prime finder?? - Namir - 08-03-2016 11:07 AM
I found on the Internet Java, C++, and Python listings for the Solovay algorithm. I translated the C++ code into Excel VBA code. The code worked fine except it would let a few composite numbers slip as primes. I added a loop to filter for the first few primes (from 2 to 19) and that make the code much better. I will post the VBA code later today.
Namir
RE: Fast prime finder?? - Namir - 08-03-2016 11:50 AM
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
RE: Fast prime finder?? - Namir - 08-04-2016 11:54 PM
I have working code for the HP-Prime which implements the Solovay algorithm. I will make a brief mention of this algorithm in one of my HHC2016 talk.s
Namir
RE: Fast prime finder?? - Claudio L. - 08-05-2016 12:46 AM
(08-04-2016 11:54 PM)Namir Wrote: I have working code for the HP-Prime which implements the Solovay algorithm. I will make a brief mention of this algorithm in one of my HHC2016 talk.s
Namir
According to Wikipedia, Solovay-Strassen was superseded by Miller-Rabin and Baillie-PSW methods.
Miller-Rabin is one of the most widely used and studied.
https://en.wikipedia.org/wiki/Miller–Rabin_primality_test
Depending on the number of primes and the values of the primes you test against, it works better or worse. Some people spent a lot of time looking for the best values to maximize its performance:
https://miller-rabin.appspot.com
The best cases are quite amazing: with only one test it can find all primes up to 341531 without letting anything slip, while only 7 tests guarantee perfect results up to 64-bit integers.
Baillie-PSW seems like a cocktail of 4 other filters, one after another to get even better results, at the expense of slightly higher complexity.
|