_IsPrime - need (very advanced !) help on Primes

jennico · September 27, 2008

hello world !

i am working hard on an optimized _IsPrime() function and i can't find a way.

the conditions: an algorithm to determine primality for every integer within AutoIt limitation ( 1,024,000,000,000,000 or 2^64 or ~ 10^15 ) within max. ~ 0.1 seconds and without using any prime list.

i made an example with three different algorithms which i sorted out to be best. their results are identical and correct. the codes are tested and correct, though maybe not optimized.

the first one (_dividedown) is easy to understand and definitely the fastest for integers up to about 10^7, but then it looses speed and the calculation time grows rapidly up to some minutes which is not acceptable. the time grows exponentially depending on the size of the result of the first division. in spite of this hazard the function can handle all integers up to the limit.

the other two functions are far more complicated and hardly understandable although they are very short and interesting. i went deep into arithmetic theory for them. the basic ideas were developed by such great men as de Fermat, Goldbach, Jacobi, Riemann. i added links to informational web pages. these two functions are designed to work with REALLY big integers far beyond 10^16, maybe 10^200, and they solve this incredible work within some thousand parts of a second.

my problem is that _SPRP and _Miller_Rabin fail with AutoIt because they use temporary results beyond the limits, for example large exponents with 6 or 7 digits length. i found quite a genial solution for that which is called Modular (Binary) Exponentiation (_ModEx). but the codes still fail on certain circumstances. i put messageboxes in the function where it comes to overflow situations. this happens on multiplication of big integers. generally speaking these problems arise with integers > Sqrt(2^64) = 32,000,000 which means that you potentially encounter false negative results with numbers > 32,000,000.

i hope anyone is able to show me a way how to calculate beyond this limit, maybe by splitting the big integers, maybe thru binary operations on them or a different trick.

any help will be appreciated.

thx j.

#cs
;#=#Title#==============================================================================================================================#
;#  Name ..........: _IsPrime.au3                                                                                                       #
;#  Description....: Examples for included functions about 'IsPrime($s_number)'                                                         #
;#  Author ........: jennico (jennicoattminusonlinedotde)                                                                               #
;#  Date ..........: 27.9.08                                                                                                            #
;#======================================================================================================================================#
#ce

$off = 1
$limit = 32000000   ;   enter your range and change this part for your needs

For $i = $off To $limit
    $timer = TimerInit()
    $x1 = _DivideDown($i)
    $t1 = Round(TimerDiff($timer) / 1000, 6)
    $timer = TimerInit()
    $x2 = _SPRP($i)
    $t2 = Round(TimerDiff($timer) / 1000, 6)
    $timer = TimerInit()
    $x3 = _Miller_Rabin($i)
    $t3 = Round(TimerDiff($timer) / 1000, 6)
    MsgBox(0, " Number " & $i & " is Prime ?", "Divide Down" & @TAB & "= " & $x1 & "  Time = " & _
            $t1 & " sec" & @CRLF & _
            "SPRP Test" & @TAB & "= " & $x2 & "  Time = " & $t2 & " sec" & @CRLF & _
            "Miller Rabin" & @TAB & "= " & $x3 & "  Time = " & $t3 & " sec")
Next

#cs
;#=#Function#===========================================================================================================================#
;#  Name ..........: _DivideDown ( $s_number )                                                                                          #
;#  Description....: Checks if Number is Prime, Result Probability 100 %                                                                #
;#  Parameters.....: $s_number = must be a positive Integer > 1                                                                         #
;#  Return Value ..: Success:   Returns 1 if Number is Prime                                                                            #
;#                   Failure:   Returns 0 and sets @error to 1 if Number is Composite                                                   #
;#                              Returns -1 and sets @error to -1 if Number is not an Integer or < 2                                     #
;#  Author ........: jennico (jennicoattminusonlinedotde)                                                                               #
;#  Date ..........: 18.9.08                                                                                                            #
;#  Remarks .......: Integers < 2 are excluded by Definition                                                                            #
;#                   _DivideDown is not a probability test, results are 100 % true.                                                     #
;#                   Unlike most common tests, _DivideDown is very fast, highly efficient, always true and barely consumes pc ressources    #
;#  Related .......: _PrimeNext, _PrimePrevious                                                                                         #
;#  Example........: yes                                                                                                                #
;#======================================================================================================================================#
#ce
Func _DivideDown($s_number)
    If IsInt($s_number) = 0 Or $s_number < 2 Then Return SetError(-1, 0, -1)
    If $s_number = 2 Then Return 1
    If Mod($s_number, 2) = 0 Then Return SetError(1, 0, 0)
    Local $i = 3, $s_limit = $s_number
    While $s_limit > $i
        If Mod($s_number, $i) = 0 Then Return SetError(1, 0, 0)
        $s_limit = $s_number / $i
        $i += 2
    WEnd
    Return 1
EndFunc   ;==>_DivideDown
#cs
;#=#Function#===========================================================================================================================#
;#  Name ..........: _SPRP ( $s_number )                                                                                                #
;#  Description....: Checks if a Number is a Strong Probable Prime (SPRP)                                                               #
;#  Parameters.....: $s_number = Number to be checked, must be a positive odd Integer > 2                                               #
;#  Return Value ..: Success:   Returns 1 if $s_number is an SPRP                                                                       #
;#                   Failure:   Returns 0 if $s_number is not an SPRP (and therefore not prime either)                                  #
;#                              Returns -1 and sets @error to -1 if Number is not an Integer or < 2                                     #
;#  Author ........: jennico (jennicoattminusonlinedotde)                                                                               #
;#  Date ..........: 18.9.08                                                                                                            #
;#  Remarks .......: Integers < 2 are excluded by Definition                                                                            #
;#                   True if N is a strong probable prime (SPRP) base B. False if N is composite, or if either N or B < 2.              #
;#                   Find d and s that satisfy N - 1 = 2^s * d, where d is odd and s >= 0.                                              #
;#                   If i = 1 or N - 1, then N is a SPRP base B.                                                                        #
;#                   Otherwise, see if i^(2^r) mod N = N - 1, where 0 <= r < s.                                                         #
;#                   It is based on Pierre de Fermat's "Theorem on Sums of two Squares (1634) and succeeded to factorize 2027651281     #
;#                   The single SPRP test itself is a weak primality test, but the combination of several subsequent tests              #
;#                   on the same number with different distinctive bases will PROVE primality (100% true).                              #
;#                   "If n < 10^16 is a 2, 3, 7, 61 and 24251 = SPRP, then n is 46,856,248,255,981 or n is prime." [Greathouse]         #
;#                   SPRP base B - Write:                                                                                               #
;#                   N - 1 = 2^s * d, with s >= 0 and d odd, then N is probably prime if either:                                        #
;#                   B^d = 1 (mod N), or (B^d)^(2^r) = -1 (mod N), for 0 <= r < s.                                                      #
;#                   http://primes.utm.edu/prove/prove2_3.html                                                                          #
;#  Related .......: _Factorize, __Trial, __Pollard_p, __Pollard_Rho, __Pollard_Strassen _IsPrime, _PrimeNext, _PrimePrevious, _gcd     #
;#  Example........: yes                                                                                                                #
;#======================================================================================================================================#
#ce
Func _SPRP($s_number)
    If IsInt($s_number) = 0 Or $s_number < 2 Then Return SetError(-1, 0, -1)
    If $s_number = 2 Then Return 1
    If Mod($s_number, 2) = 0 Then Return 0
    Local $B, $N = $s_number, $s = 0, $d, $i, $j, $r, $s_base[5] = [24251, 61, 7, 3, 2]
    Do
        $s += 1
        $d = BitShift($N - 1, $s)
    Until Mod($d, 2)
    For $j = 0 To 4
        $B = $s_base[$j]
        If $N = $B Then Return 1
If $B > 10 ^ 14 Then MsgBox(0, "SPRP first", "overflow" & @CRLF & $B);catch overflow
        $i = _ModEx($B, $d, $N)
        If $i = 1 Or $i = $N - 1 Then ContinueLoop
        For $r = 0 To $s - 1
If $i > 10 ^ 14 Then MsgBox(0, "SPRP second", "overflow" & @CRLF & $i);catch overflow
            $i = _ModEx($i, 2, $N)
            If $i = $N - 1 Then ExitLoop
            If $r = $s - 1 Then Return 0
        Next
    Next
    Return 1
EndFunc   ;==>_SPRP
#cs
;#=#Function#===========================================================================================================================#
;#  Name ..........: _Miller_Rabin ( $s_number )                                                                                        #
;#  Description....: Checks if a Number is probable Prime (PRP)                                                                         #
;#  Parameters.....: $s_number = Number to be checked (Integer)                                                                         #
;#  Return Value ..: Success:   Returns 1 if $s_number is probable prime                                                                #
;#                   Failure:   Returns 0 if $s_number is definitely not prime                                                          #
;#                              Returns -1 and sets @error to -1 if Number is not an Integer or < 2                                     #
;#  Author ........: jennico (jennicoattminusonlinedotde)                                                                               #
;#  Date ..........: 27.9.08                                                                                                            #
;#  Remarks .......: Integers < 2 are excluded by Definition                                                                            #
;#                   The Miller-Rabin Test (1974) is derived from $PRP test and uses a similar algorithm, but is much faster.           #
;#                   It bases on the extended Riemann hypothesis which is not yet proved theoretically.                                 #
;#                   On the other hand the result is weak for all results = True (but 100% true for result=False).                      #
;#                   Instead of using destinctive bases like SPRP, the MRT uses 10-30 RANDOM bases.                                     #
;#                   In this way the trustness is raised to almost 100% (maybe you find an error once in your life).                    #
;#                   The MRT is so trustworthy that is was used for encryption of your VISA cards, until the RSA algorithm was found.   #
;#                   "Multiple MRTs using the first 7 primes (using 8 gives no improvement) are valid for every number up to 10^14."    #
;#                   http://mathworld.wolfram.com/Rabin-MillerStrongPseudoprimeTest.html                                                #
;#  Related .......: _Factorize, __Trial, __Pollard_p, __Pollard_Rho, __Pollard_Strassen _IsPrime, _PrimeNext, _PrimePrevious, _gcd     #
;#  Example........: yes                                                                                                                #
;#======================================================================================================================================#
#ce
Func _Miller_Rabin($s_number)
    If IsInt($s_number) = 0 Or $s_number < 2 Then Return SetError(-1, 0, -1)
    If $s_number = 2 Then Return 1
    If Mod($s_number, 2) = 0 Then Return 0
    Local $t, $u, $a, $j, $i, $r, $v, $N = $s_number, $b[7] = [2, 3, 5, 7, 11, 13, 17]
    Do
        $t += 1
        $u = ($N - 1) / 2 ^ $t
    Until Mod($u, 2)
    For $j = 0 To 6
        $a = $b[$j]
        If $a = $N Then Return 1
        If _ModEx($a, $u, $N) = 1 Then ContinueLoop
        For $i = 0 To $t - 1
            $v = _ModEx(2, $i, $N)
If $v * $u > 10 ^ 14 Then MsgBox(0, "miller-rabin", "overflow" & _
        @CRLF & $v & @CRLF & $u & @CRLF & $t - 1);catch overflow
            $v = _ModEx($a, $v * $u, $N)
            If $v = -1 Or $v = $N - 1 Then ExitLoop
        Next
        If $i = $t Then Return 0
    Next
    Return 1
EndFunc   ;==>_Miller_Rabin
#cs
;#=#Function#===========================================================================================================================#
;#  Name ..........: _ModEx ( $s_base, $s_exp, $s_modulus )                                                                             #
;#  Description....: Modular_right-to-left_binary_exponentiation_algorithm, Model: ==> Mod($s_base^($s_exp),$s_modulus)                 #
;#  Parameters.....: $s_base = Base of exponentiation                                                                                   #
;#                   $s_exp = Exponent of exponentiation                                                                                #
;#                   $s_modulus = the divisor for the modular process                                                                   #
;#  Return Value ..: Returns Modulus (Remainder of Division), an Integer                                                                #
;#  Author ........: jennico (jennicoattminusonlinedotde)                                                                               #
;#  Date ..........: 18.9.08                                                                                                            #
;#  Remarks .......: No error handling in this subfunction !                                                                            #
;#                   This function should only be used internally as an auxiliary subfunction to factorize big Integers                 #
;#                   Use right-to-left binary exponentiation to Calculate B^d mod N, result in i.                                       #
;#                   Example: Mod(3^45678786756453,45676545656) cannot be calculated due to System limitations (stack overflow)         #
;#                   Therefore, the exponent will be processed as a binary string in a loop from right (Low) to left (High Bit)         #
;#                   The Modulo operation is commutative                                                                                #
;#  Related .......: _Factorize, __Trial, __Pollard_p, __Pollard_Rho, __Pollard_Strassen _IsPrime, _PrimeNext, _PrimePrevious, _gcd     #
;#  Example........: no                                                                                                                 #
;#======================================================================================================================================#
#ce
Func _ModEx($s_base, $s_exp, $s_modulo, $ret = 1)
If $s_base > 10 ^ 14 Then MsgBox(0, "basefirst", "overflow" & @CRLF & $s_base);catch overflow
    While 1
If $s_base > 10 ^ 14 Then MsgBox(0, "base", "overflow" & @CRLF & $s_base);catch overflow
        If Mod($s_exp, 2) Then $ret = Mod($ret * $s_base, $s_modulo)
If $ret > 10 ^ 14 Then MsgBox(0, "ret", "overflow" & @CRLF & $ret);catch overflow
        If $s_exp < 2 Then Return $ret
        $s_exp = BitShift($s_exp, 1)
        $s_base = Mod($s_base * $s_base, $s_modulo)
    WEnd
EndFunc   ;==>_ModEx

Edited September 27, 2008 by jennico

jennico · September 27, 2008

this is part of a bigger script, that's why you can't find the related functions.

j

Edited September 28, 2008 by jennico

Malkey · September 28, 2008

hello world !

i am working hard on an optimized _IsPrime() function and i can't find a way.

For a long time , I have used this function to test for a prime number. If the input number is the same as the output number, the number is prime.

It is slower than your examples.

I would like you to report back about this function. Tell me what you think. You may be albe to optimize it for your use if you think it is worth spending time doing so.

I am hoping it will be helpful.

$i = 32000000
$timer = TimerInit()
$x1 = Prim($i)
$t1 = Round(TimerDiff($timer) / 1000, 6)

MsgBox(0, " Number " & $i & " Factors are ", "Factors" & @TAB & "= " & $x1 & "  Time = " & _
        $t1 & " sec" & @CRLF)

Func Prim($a)
    $factors = ""
    Dim $b, $c
    If $a > 0 Then
        If IsNumber($a) Then
            $a = Int($a)
            While $a / 2 - Int($a / 2) = 0
                $a = $a / 2
                $factors = $factors & "2 "
            WEnd
            $b = 3
            $c = ($a) * ($a) + 1
            Do
                If $a / $b - Int($a / $B) = 0 Then
                    If $a / $b * $b - $a = 0 Then
                        $factors = $factors & $b & " "
                        $a = $a / $b
                    EndIf
                EndIf
                If $a / $b - Int($a / $B) <> 0 Then $b = $b + 2
                $c = ($a) * ($a) + 1
            Until $b > $c Or $b = $c

            If $a <> 1 Then $factors = $factors & $a & " "
        EndIf
    EndIf
    ;MsgBox(0,"Prim",  $factors)
    Return $factors
EndFunc   ;==>Prim

Edit: I think it may be the same as your _DivideDown() example.

Edited September 28, 2008 by Malkey

jennico · September 28, 2008

thank you for your contribution !

it is a very interesting approach and uses trial division starting with 2,3 and then +2. you can reduce the work by additionally skipping the integers that end with 5. there is an even more efficient way in repeatedly taking the steps: 2-6-4-2-4-6-2. but trial division will always be toooooooooo slow on large integers like below.

it was very fast in finding out that 737373737373737 is prime (my favourite prime). bravo

but it is still working since some minutes on:

100000000003

1000000000039

10000000000037

100000000000031

1000000000000037

i think your code can be optimized, you do not have to calculate the same process more than once, better make a variant, and IsInt is better than the difference of number and integer, mod() is faster than dividing, bitshift is faster than deviding by 2 and so on. in recursive functions one has to reduce the calculating steps for the cpu.

and, of course, the factoring part is great. this is my next aim: finding a factoring algorithm, very fast and with the range of 10^15.

cheers j.

Edited September 28, 2008 by jennico

jennico · September 28, 2008

ja, might be the same as divide down ...... not sure.

i am going to implement the elliptic curve factorization (ECM) and the Quadratic sieve (QS) maybe they do it better...

other languages like delphi and python have very fast built in functions for _isPrime and _getfactors but i can't find a scource code.

j.

ah, my _dividedown needs 4 minutes for 737373737373737 while your code gets it instantly. i have to find out why !

Edited September 28, 2008 by jennico

jennico · September 28, 2008

ain't nobody here who knows about modular arithmetic ?

j.

enaiman · September 28, 2008

Here is my code for testing prime numbers:

$number = 737373737373737 
_isPrime($number)

Func _isPrime($nb)
    $j=2
    $not_prime = 0
    Do 
        If Mod($nb, $j) = 0 Then 
            $not_prime = 1
            ExitLoop
        EndIf
        If $j < 3 Then
            $j +=1
        Else
            $j +=2  
        EndIf
    Until $j > $nb/$j   
    If $not_prime = 1 Then 
        ConsoleWrite("Number "&$nb&" is not a prime number" & @CRLF)
    Else
        ConsoleWrite("Number "&$nb&" is a prime number" & @CRLF)
    EndIf
EndFunc

It needs 137 seconds on my machine to return the answer.

Another script will find all prime numbers to a specified number and write them on the console (based on the same algorythm).

HotKeySet("{end}", "Terminate")
HotKeySet("{ins}", "calc");Shift-Alt-r
$i=0
$upNbr = 100000
ConsoleWrite(@HOUR & ':' & @MIN & ':' & @SEC & @CRLF)
calc()
While 1
    Sleep(350)
WEnd
Func calc()
    While 1
        $i += 1
        $j=2
        $not_prime = 0
        Do 
            If Mod($i, $j) = 0 Then 
                $not_prime = 1
                ExitLoop
            EndIf
            If $j < 3 Then
                $j +=1
            Else
                $j +=2  
            EndIf
        Until $j > $i/$j    
        If $i > $upNbr Then Terminate()
        If $not_prime = 1 Then 
            ContinueLoop
        Else
            ConsoleWrite($i & @CRLF)
        EndIf
    WEnd
EndFunc ;==>calc

Func Terminate()
    ConsoleWrite(@HOUR & ':' & @MIN & ':' & @SEC & @CRLF)
    Exit 0
EndFunc ;==>Terminate

jennico · September 29, 2008

if you want to list all primes within a limit there is a much faster method called the Sieve of Eratosthenes (ca 240 b.c).

you dim an empty array of so many indices as your limit. then you begin with 2 and set all arrays whose indices are multiples of 2 to 1: array[4]=1,array[6]=1, array[8]=1 and so on. then you take 3 and do the same. if an array is already set to 1 you can skip it.

you take 4, 5, 6 and so on up to the square root of your limit +1, setting all multiples to 1. you do not need to calculate above the squareroot, because all numbers before already cover the complete array.

when done, you go thru your "sieved" array, and the indices that are not set to 1 are primes. voila !

this is so fast, you won't believe it.

you can speed it up by leaving all even indices away because 2 is the only even prime. in this way you can extend the array to 32,000,000 entries - again, the magic number from above, the square root of 2^64. but this still does not cover all integers up to 2^64 like i want it to.

and my basic problem is still the modular multiplication:

Func _ModEx($s_base, $s_exp, $s_modulo, $ret = 1)
    While 1
        If Mod($s_exp, 2) Then $ret = Mod($ret * $s_base, $s_modulo)
        If $s_exp < 2 Then Return $ret
        $s_exp = BitShift($s_exp, 1)
        $s_base = Mod($s_base * $s_base, $s_modulo)
    WEnd
EndFunc   ;==>_ModEx

the problem is that the lines

If Mod($s_exp, 2) Then $ret = Mod($ret * $s_base, $s_modulo) and
$s_base = Mod($s_base * $s_base, $s_modulo)

can exceed the 10^15 limit and i need a method to bypass the limit

other languages know the _MulMod function which enables calculation beyond the limit: (MulMod = modular multiplication)

// Compute a*b modulo _m
inline ModInt MulMod(ModInt a, ModInt B )
{
// classical trick to bypass the 64-bit limitation, when a*b does not fit into the ModInt type.
// Works whenever a*b/_m is less than 2^52 (double type maximal precision)
ModInt q = (ModInt) (_invm*(double)a*(double)B );
return a*b-q*_m;
}

i need help to implement this trick in my code, like it has been done for other languages.

this is the bottleneck. when it is solved, i can write UDFs that factorizes and proves primality within 5 thousands of a second, no matter how large the integer is.

who can help ?

j.

edit: there's another scource i found for it:

// Calculates Src1^Src2 mod Modulus
// Assembly needed to prevent overflow without Int64
function MulMod(Src1, Src2, Modulus: LongWord): LongWord;
asm
    mul edx
    div ecx
    mov eax,   edx
end;
 
// Calculates Base^Exponent mod Modulus
function PowerMod(Base, Exponent, Modulus: LongWord): LongWord;
begin
    Result := 1;
    while Exponent &gt; 0 do
    begin
        if Exponent and 1 &lt;&gt; 0 then
            Result := MulMod(Result, Base, Modulus);
        Exponent := Exponent shr 1;
        Base := MulMod(Base, Base, Modulus);
    end;
end;

the second function is the same as my _ModEx. what i need is the MulMod above that. in this example directly assembler is used, but i think this can't be done with autoit.

Edited September 29, 2008 by jennico

AdmiralAlkex · September 29, 2008

@jennico

If you managed to build MulMod in machine code you could use THIS, that would be very fast.

Edited September 29, 2008 by AdmiralAlkex

jennico · October 1, 2008

Yeah, i did it !

@ first:

i managed MulMod the AutoIt way (without assembler). these modular arithmetic things are very fast, i would even prefer them to "normal" division and multiplication. i will make a complete UDF for that when i finished optimizing the codes.

BUT: the reason why my SPRP and RMT examples failed on large numbers, was a completely different issue. From the help file (well hidden):

Bit operations are performed as 32-bit integers.

so, when it comes to calculate integers > 4,294,967,296 you cannot use BitShift, BitOr, BitAnd anymore !!!

haha. well i tested my codes over and over and i could not find a mistake until i suddenly found this restriction.

in fact, that means that i have to implement a manual IntToBin function for integers. of course this costs calculation speed and it's kind of stupid to manually convert an int to bin, since the pc does this anyway with every input. at least my _isprime function now works up to 16 digits, but it has slowed down to about half a second which is still not optimal (my aim was 0.1 sec) but kind of nearly acceptable.

so i would like to ask the autoit team to remove the restrictions on bit operations. please !!!!

@AdmiralAlkex :

thank you for this link. seems that this is what i need. i will test if it is a suitable and fast assembler interface. it sounds great so far. thank you indeed.

i will post all my results when ready

j.

_IsPrime - need (very advanced !) help on Primes

Recommended Posts

Link to comment

Share on other sites

Link to comment

Share on other sites

Link to comment

Share on other sites

Link to comment

Share on other sites

Link to comment

Share on other sites

Link to comment

Share on other sites

Link to comment

Share on other sites

Link to comment

Share on other sites

Link to comment

Share on other sites

Link to comment

Share on other sites

Create an account or sign in to comment

Create an account

Sign in

Recently Browsing 0 members