jennico Posted September 27, 2008 Share Posted September 27, 2008 (edited) 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. expandcollapse popup#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 Spoiler I actively support Wikileaks | Freedom for Julian Assange ! | Defend freedom of speech ! | Fight censorship ! | I will not silence.Don't forget this IP: 213.251.145.96 Link to comment Share on other sites More sharing options...
jennico Posted September 27, 2008 Author Share Posted September 27, 2008 (edited) this is part of a bigger script, that's why you can't find the related functions. j Edited September 28, 2008 by jennico Spoiler I actively support Wikileaks | Freedom for Julian Assange ! | Defend freedom of speech ! | Fight censorship ! | I will not silence.Don't forget this IP: 213.251.145.96 Link to comment Share on other sites More sharing options...
Malkey Posted September 28, 2008 Share Posted September 28, 2008 (edited) 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. expandcollapse popup$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 Link to comment Share on other sites More sharing options...
jennico Posted September 28, 2008 Author Share Posted September 28, 2008 (edited) 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 Spoiler I actively support Wikileaks | Freedom for Julian Assange ! | Defend freedom of speech ! | Fight censorship ! | I will not silence.Don't forget this IP: 213.251.145.96 Link to comment Share on other sites More sharing options...
jennico Posted September 28, 2008 Author Share Posted September 28, 2008 (edited) 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 Spoiler I actively support Wikileaks | Freedom for Julian Assange ! | Defend freedom of speech ! | Fight censorship ! | I will not silence.Don't forget this IP: 213.251.145.96 Link to comment Share on other sites More sharing options...
jennico Posted September 28, 2008 Author Share Posted September 28, 2008 ain't nobody here who knows about modular arithmetic ? j. Spoiler I actively support Wikileaks | Freedom for Julian Assange ! | Defend freedom of speech ! | Fight censorship ! | I will not silence.Don't forget this IP: 213.251.145.96 Link to comment Share on other sites More sharing options...
enaiman Posted September 28, 2008 Share Posted 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). expandcollapse popupHotKeySet("{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 SNMP_UDF ... for SNMPv1 and v2c so far, GetBulk and a new example script wannabe "Unbeatable" Tic-Tac-Toe Paper-Scissor-Rock ... try to beat it anyway :) Link to comment Share on other sites More sharing options...
jennico Posted September 29, 2008 Author Share Posted September 29, 2008 (edited) 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 > 0 do begin if Exponent and 1 <> 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 Spoiler I actively support Wikileaks | Freedom for Julian Assange ! | Defend freedom of speech ! | Fight censorship ! | I will not silence.Don't forget this IP: 213.251.145.96 Link to comment Share on other sites More sharing options...
AdmiralAlkex Posted September 29, 2008 Share Posted September 29, 2008 (edited) @jennicoIf you managed to build MulMod in machine code you could use THIS, that would be very fast. Edited September 29, 2008 by AdmiralAlkex .Some of my scripts: ShiftER, Codec-Control, Resolution switcher for HTC ShiftSome of my UDFs: SDL UDF, SetDefaultDllDirectories, Converting GDI+ Bitmap/Image to SDL Surface Link to comment Share on other sites More sharing options...
jennico Posted October 1, 2008 Author Share Posted 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 readyj. Spoiler I actively support Wikileaks | Freedom for Julian Assange ! | Defend freedom of speech ! | Fight censorship ! | I will not silence.Don't forget this IP: 213.251.145.96 Link to comment Share on other sites More sharing options...
Recommended Posts
Create an account or sign in to comment
You need to be a member in order to leave a comment
Create an account
Sign up for a new account in our community. It's easy!
Register a new accountSign in
Already have an account? Sign in here.
Sign In Now