 Fast prime finder?? - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: Fast prime finder?? (/thread-6619.html) 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
"          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.