Fast prime finder??

[Namir]

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

[Claudio L.]

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

[Namir]

I am looking for the algorithm(s).

Namir

[Namir]

I found what I was looking for on the Internet.

:-)

Namir

[Namir]

(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

[rprosperi]

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

[ttw]

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.

[Namir]

(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

[Namir]

(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

[EdS2]

(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

[ttw]

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.

[Namir]

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

[Namir]

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

[Namir]

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

[Claudio L.]

(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–Rab...ality_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.