All about Primes - The Prime Suite - _primes.au3 UDF lots of math functions

[jennico]

All about Primes - The Primes Suite

everything began when i tried to find an algorithm for prime factor(iz)ing integers of any size within a reasonable time (let's say 0.1 seconds). i had to learn that this issue is still unsolved and will be unsolved forever. at last, the problem depends on the recognition of the prime factor itself.

The Fundamental Theorem of Arithmetic: "Every natural number is either prime or can be uniquely factored as a product of primes in a unique way."

my first attempt was a prime list based approach. if you create a list of all primes within the AutoIt limitation, you could easily create a factorizing algorithm out of it. so i implemented a fascinating fast ancient algorithm, the "Sieve of Eratosthenes", ca. 300 BC. since you only have to calculate primes up to the square root of the limitation, the numbers below 32,000,000 should be sufficient.

the list, saved to a textfile, resulted in incredible 16 mb and it took some 15 minutes to be created. okay, once you have the list, the rest is easy. but still not good enough, because only reading in the list takes about three seconds.

so i kept searching for better results and made my way deep into math science. i came across the Euclidean Algorithm, the Fermat theorem, the Goldbach conjecture, the Riemann hypothesis, Gauss and Carmichael numbers, the Zeta and Beta function and modern encryption algorithms, all of them deal with the same problem: how to determine that a number is a prime ?

the algorithms i finally implemented in my UDF are very strange: they are nearly impossible to explain and to calculate by hand, but transposed to computer language they result incredibly simple, just one or two loops and hardly 10 lines of codes. nevertheless the implementation is complicated. the numbers are looped some thousand times through, you have to deal with huge exponents that break every limit. in doing so you have to use a very interesting modern part of arithmetics called modular arithmetics.

please observe that you will, of course, not find any kind of primes list within the UDF, this would have been unfair.

on my way i came across a whole bunch of algorithms representing 2500 years of algebra, that i put together in one file. many of them i would judge very useful for other implementations.

i am not sure if they are all optimized, so be free to make suggestions if you find an improvement. on making prime functions, all that matters in the end is the calculation speed.

the zip file (Primes.zip) contains:

- _Primes.au3 : the UDFs

- _PrimesExamples.au3 : Examples for the UDFs

- PrimeSpiral.au3 : an implementation that visualizes primes in a very clever way based on Ulam's Prime Spiral

Primes.au3 contents:

Main Functions :

_GetFactors ( $s_number )

_IsPrime ( $s_number )

_PrimeNext ( $s_number )

_PrimePrevious ( $s_number )

_RandomPrime ( $s_max [, $s_min] )

_IsTwinPrime ( $s_number )

_PrimeGap ( $s_number [, $s_mode] )

_PrimesInRange ( $s_max [, $s_min] )

_Triangular ( $s_count )

_IsTriangular ( $s_number )

_TriangularsInRange($s_max [, $s_min] )

Prime Related Math Functions :

_GCD ( $s_number1, $s_number2 = "-", $s_number3 = "-", $s_number4 = "-", $s_number5 = "-" )

(Greatest Common Divisor)

_ggT ( $s_number1, $s_number2 = "", $s_number3 = "", $s_number4 = "", $s_number5 = "" )

_LCM ( $s_number1, $s_number2 = "-", $s_number3 = "-", $s_number4 = "-", $s_number5 = "-" )

(Least Common Multiple)

_kgV ( $s_number1, $s_number2 = "", $s_number3 = "", $s_number4 = "", $s_number5 = "" )

_IsEven ( $s_number )

_IsOdd ( $s_number )

_Cross ( $s_number, $s_mode=0 )

_Factorial ( $s_number )

_LogX ( $s_Result, $s_Base )

_LogE ( $s_Result )

_Log10 ( $s_Result )

_Log2 ( $s_Result )

_RootX ( $s_Number, $s_Base )

_Swap ( ByRef $s_value1, ByRef $s_value2 )

Bit Operations :

_BitAND ( $s_value1 [, $s_value2 = 0] )

_BitLen ( $s_value )

_BitMid ( $s_value, $s_start, $s_count = 1000000000000000 )

_BitNOT ( $s_value )

_BitOR ( $s_value1 [, $s_value2 = 0] )

_BitShift ( $s_value, $s_shift )

_BitXOR ( $s_value1 [, $s_value2 = 0] )

_ModEx ( $s_base, $s_exp, $s_modulus )

_MulMod ( $x, $y, $N )

_Binary ( $s_number )

Listbased :

_FactorizeL ( $s_number, $s_mode = 0, $s_separator = "", $s_limit = 32000000 )

_CreatePrimeList ( [$s_limit = 32000000] )

_LoadPrimeList ( [$mode = 0, $file = "Primes-32000000.txt"] )

_GetFactorsL ( $s_number [,$a_Prime = 0] )

_CreatePrimeListEx ( ) NOT working

_CreatePrimeListMax ( ) NOT working

Numeric Conversions :

_IntToBin ( $s_number )

_BinToInt ( $s_string )

_IntToHex ( $hx_dec [, $hx_length] )

_HexToInt ( $hx_hex )

_IntToAny ( $s_number, $s_base [, $s_length] )

_AnyToInt ( $s_string, $s_base )

_AnyToAny ( $s_string, $s_base1, $s_base2 [, $s_length] )

_Arabian ( $input )

_Roman ( $input )

_1000 ( $x, $s_separator = "" )

Subfunctions :

_Eratosthenes ( $s_max )

_Euclid ( $a, $b )

_DivideDown ( $s_number )

_SPRP ( $s_number )

_Miller_Rabin ( $s_number )

_Fermat ( $s_number )

_Pollard_p ( $s_number )

_Pollard_Rho ( $s_number )

_Pollard_Strassen ( $s_number )

_Atkins ()

_GetFactorsOld ( $s_number )

_IsPrimeOld ( $s_number )

_PrimesExamples.au3

Maybe anyone can tell me why "PrimeSpiral" is always started twice ? i can't find out why.

PrimeSpiral.au3

In 1963, Stanislav Ulam discovered a fascinating way to visualize primes. beginning with number 1 in the center, write down all consequent numbers in a spiral.

17---*---*---*--13   ...
 |             |   |
 *   5---*---3   *   29
 |   |     |   |   |
19  *   *---2  11   *
 |   |         |   |
 *   7---*---*---*   *
 |                 |
 *---*--23---*---*---*

mark all primes in a different color and the result is a suprising fractal-like structure. you will see big numbers of primes lined up in diagonal lines, and others vertically or horizontally. it is surprising because it seems to implement some kind of rule in prime distribution. But still, there has not been found one.

you can choose other base numbers, but the result will always be alike.

my program can also mark square and triangular numbers in this matrix. to see them clearly, choose a blue color like the background for the primes. you can also mark twin primes or differentiate the primes by color. or you can set a fixed diagonal and let the bases cycle. there's a zoom function and some other tools.

the algorithm to calculate the primes is the Sieve of Eratosthenes, so this is also a demonstration of how this works. the first ten numbers are quite slow (because the lower the number, the more often it has to strike out its multiples), but then the spiral explodes and will fill your screen within 50 seconds. the program is faster if you click on the child guis, because the graphic will not be refreshed then. when the flickering begins, it means that the refreshing needs more time than calculating a new square of numbers. so this is impressingly fast !

everytime you start the program you should leave it build at least once a complete full screen run, so that the prime matrix can be filled. otherwise you will find false primes on the next runs.

as an example exploration, set offset (=base) to primes 17 or 41. you can find a long continous prime diagonal (NE-SW) through the center, that ends on both sides right at the starting points of the two square number diagonals (NW-SE), so that these 3 diagonals build a perfect, right angled, continuous Z. (still this is not the famous Zeta function, but it's close !) There should be more offset primes that build the same pattern, with an increasing central prime diagonal size.

okay, have fun and interesting conclusions with my programs.

jennico

Edit: 23.10. : updated _Primes.au3 : v1.2 (corrected _Factorial, improved _IsPrime, _PrimeNext and PrimePrevious, corrected loop declarations) prev downloads : 23

Primes.zip

Edited October 23, 2008 by jennico

[trancexx]

I was expecting something big, but this is huge! Impressive!

Please correct _Factorial() function and definition.

[jennico]

thank you. seems i was confused.

j.

corrected and updated.

Edited October 22, 2008 by jennico

[JRowe]

http://www.claymath.org/millennium/Riemann_Hypothesis/

If you can correlate the distribution of primes with the Riemann zeta function, and offer a proof of it, it's worth a cool million. And math history. Probably a Fields prize, as well (think Nobel prize in mathematics.)

[madflame991]

Good job!

I'm also interested in fast algorithms for prime numbers

While $s_limit > $i
    If Mod($s_number, $i) = 0 Then Return
    $s_limit = $s_number / $i
    $i += 2
WEnd

this can be optimized into this:

$s_limit = Round(Sqrt($s_number))
While $s_limit > $i
    If Mod($s_number, $i) = 0 Then Return
    $i += 2
WEnd

this way you avoid doing $s_limit = $s_number / $i in the loop

[jennico]

well, no...

just try it with a big prime, let's say 31,999,939. without the division it will calculate about two minutes.

you are right for smaller numbers and negative results, but when it comes to a high prime, it looses.

in fact, the division limits the amount of necessary tests to a minimum, though it might slow down the function for other cases. but my attempt was to achieve a TRUE result in any case within some parts of a second.

thank you for this try anyway ! i am pleased of people trying to help !

j.

[jennico]

here's a list of some bigger primes for testing:

1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063

909091

858577

8353021

37956673

71717117

2497145999

737373737373737

100000000003

1000000000039

10000000000037

100000000000031

j.

[jennico]

to those who want to follow my way into the math universe, i recommend:

http://mathworld.wolfram.com

http://primes.utm.edu

wikipedia

these are the by far best and most competent sources in the web.

j.

we all know that AutoIt is not made for mathematics, and other languages would do it much faster. it would be a help if AutoIt allowed direct raw machine code calculations. i tried some of the attempts that have been made here in the forum. there is a way to calculate divisions and multiplications via a dll. but when it comes to mod, there was no solution yet.

Edited October 22, 2008 by jennico

[madflame991]

Have a look at the results

expandcollapse popup

$number = Int(InputBox('Input','Enter a number',31999939))
$timer = TimerInit()
For $i = 1 to 100
    _IsPrime1($number)
Next
$tdiff1 = TimerDiff($timer)
For $i = 1 to 100
    _IsPrime2($number)
Next
$tdiff2 = TimerDiff($timer)
MsgBox(0,'','Function 1: executed in '&$tdiff1&' miliseconds'&@CRLF&'Function 2: executed in '&$tdiff2-$tdiff1&' miliseconds')

Func _IsPrime1($s_number)
    If $s_number < 1 Then Return 
    If $s_number = 2 Then Return 1
    If $s_number = 1 Or Mod($s_number, 2) = 0 Then Return
    Local $i = 3, $s_limit = $s_number
    While $s_limit > $i
        If Mod($s_number, $i) = 0 Then Return
        $s_limit = $s_number / $i
        $i += 2
    WEnd
    Return 1
EndFunc  ;==>_IsPrime

Func _IsPrime2($n);optimized one
    If $n < 1 Then Return
    If $n = 2 Then Return 1
    If $n = 1 Or Mod($n,2) = 0 Then Return
    Local $i = 3, $l = Round(Sqrt($n))
    While $i < $l
        If Mod($n,$i) = 0 Then Return 0
        $i += 2
        If Mod($n,$i) = 0 Then Return 0
        $i += 4
    WEnd
    Return 1
EndFunc

[jennico]

on first sight it might be okay, i will try it later. what you propose is simple trial division. might be faster indeed. your idea to just calculate up to the square root is right. thanks

j.

Edited October 22, 2008 by jennico

[jennico]

there is a logical error in your function:

you begin with 3 and add 2 = 5

then you add 4 = 9

then 2 = 11

then 4 = 15 - 17 - 21 - 25 - 29 - 31 - 35 - 37 - 41 - 43 - 47 -49

the number 9 or 49 for example will result prime in your algorithm. this needs some more more attention.

expandcollapse popup

#include <_Primes.au3>

While 1
    $number = Int(InputBox('Input','Enter a number',31999939))
    $timer = TimerInit()
;   For $i = 1 to 100
        $x1=_IsPrime1($number)
;   Next
    $tdiff1 = TimerDiff($timer)
    $timer = TimerInit()
;   For $i = 1 to 100
        $x2=_IsPrime2($number)
;   Next
    $tdiff2 = TimerDiff($timer)
    $timer = TimerInit()
;   For $i = 1 to 100
        $x3=_IsPrime($number)
;   Next
    $tdiff3 = TimerDiff($timer)
    MsgBox(0,'','Function 1: executed in '&$tdiff1&' miliseconds'&$x1&@CRLF&'Function 2: executed in '&$tdiff2&' miliseconds'&$x2&@CRLF&'Function 3: executed in '&$tdiff3&' miliseconds'&$x3)
WEnd

Func _IsPrime1($s_number)
    If $s_number < 1 Then Return 
    If $s_number = 2 Then Return 1
    If $s_number = 1 Or Mod($s_number, 2) = 0 Then Return
    Local $i = 3, $s_limit = $s_number
    While $s_limit > $i
        If Mod($s_number, $i) = 0 Then Return
        $s_limit = $s_number / $i
        $i += 2
    WEnd
    Return 1
EndFunc  ;==>_IsPrime

Func _IsPrime2($n);optimized one
    If $n < 1 Then Return
    If $n = 2 Then Return 1
    If $n = 1 Or Mod($n,2) = 0 Then Return
    Local $i = 3, $l = Round(Sqrt($n))
    While $i < $l
        If Mod($n,$i) = 0 Then Return 0
        $i += 2
        If Mod($n,$i) = 0 Then Return 0
        $i += 4
    WEnd
    Return 1
EndFunc

Edited October 22, 2008 by jennico

[Malkey]

we all know that AutoIt is not made for mathematics, and other languages would do it much faster. it would be a help if AutoIt allowed direct raw machine code calculations. i tried some of the attempts that have been made here in the forum. there is a way to calculate divisions and multiplications via a dll. but when it comes to mod, there was no solution yet.

I was looking through the other math functions in msvcrt.dll file which trancexx used here

http://www.autoitscript.com/forum/index.ph...st&p=592196

I saw the fmod function and thought of your problem in the other thread about calculating with double, 64bit(8bytes) floating point, within a structure of double or float type.

I did these two wrappers for the fmod and pow functions in the msvcrt.dll file because they use double type. I just wanted to see if they worked. And gave up when I could not get modf to work.

Your comments in the above post reminded me.

I don't know if fmod() will be of any use to you jennico. Here they are.

expandcollapse popup

$a = 1023
MsgBox(4096, "", "fmod() = " & fmod(2 ^ 1019, $a) & @CRLF & "pow() = " & pow(2, $a) & @CRLF)


Func fmod(Const $nY, Const $nX)
    Local $aResult
    Local $h_DLL = DllOpen("msvcrt.dll")
    $aResult = DllCall($h_DLL, "double:cdecl", "fmod", "double", $nY, "double", $nX, "double*", 0)
    ;local $res
    ;for $x = 0 to 2
    ;   $res &= $aResult[$x] & @CRLF
    ;next
    ;ConsoleWrite($res & @CRLF)
    If @error Then
        DllClose($h_DLL)
        Return SetError(2, 0, "-1.#err")
    EndIf

    DllClose($h_DLL)
    Return SetError(0, 0, $aResult[0])

EndFunc   ;==>fmod


Func pow(Const $nY, Const $nX)
    Local $aResult
    Local $h_DLL = DllOpen("msvcrt.dll")
    $aResult = DllCall($h_DLL, "double:cdecl", "pow", "double", $nY, "double", $nX, "double*", 0)
    ;local $res
    ;for $x = 0 to 2
    ;   $res &= $aResult[$x] & @CRLF
    ;next
    ;ConsoleWrite($res & @CRLF)
    If @error Then
        DllClose($h_DLL)
        Return SetError(2, 0, "-1.#err")
    EndIf
    DllClose($h_DLL)
    Return SetError(0, 0, $aResult[0])
EndFunc   ;==>pow

[jennico]

i corrected and improved madflames idea. i think function 3 is the fastest for primes < 100000000003, for higher primes, my implemented one becomes faster. i will update this later.

expandcollapse popup

#include <_Primes.au3>

While 1
    $number = Int(InputBox('Input','Enter a number',31999939))
    $timer = TimerInit()
;   For $i = 1 to 100
        $x1=_IsPrime1($number)
;   Next
    $tdiff1 = TimerDiff($timer)
    $timer = TimerInit()
;   For $i = 1 to 100
        $x2=_IsPrime2($number)
;   Next
    $tdiff2 = TimerDiff($timer)
    $timer = TimerInit()
;   For $i = 1 to 100
        $x3=_IsPrime3($number)
;   Next
    $tdiff3 = TimerDiff($timer)
    $timer = TimerInit()
;   For $i = 1 to 100
        $x4=_IsPrime($number)
;   Next
    $tdiff4 = TimerDiff($timer)
;   $timer = TimerInit()
;   For $i = 1 to 100
;       $x5=_IsPrime($number)
;   Next
;   $tdiff5 = TimerDiff($timer)
    MsgBox(0,'','Function 1: executed in '&$tdiff1&' miliseconds'&$x1&@CRLF&'Function 2: executed in '&$tdiff2&' miliseconds'&$x2&@CRLF&'Function 3: executed in '&$tdiff3&' miliseconds'&$x3&@CRLF&'Function 4: executed in '&$tdiff4&' miliseconds'&$x4)
WEnd

Func _IsPrime1($s_number)
    If $s_number < 1 Then Return 
    If $s_number = 2 Then Return 1
    If $s_number = 1 Or Mod($s_number, 2) = 0 Then Return
    Local $i = 3, $s_limit = $s_number
    While $s_limit > $i
        If Mod($s_number, $i) = 0 Then Return
        $s_limit = $s_number / $i
        $i += 2
    WEnd
    Return 1
EndFunc  ;==>_IsPrime

Func _IsPrime2($n);optimized one
    If $n < 1 Then Return
    If $n = 2 Then Return 1
    If $n = 1 Or Mod($n,2) = 0 Then Return
    Local $i = 3, $l = Round(Sqrt($n))
    While $i < $l
        If Mod($n,$i) = 0 Then Return 0
        $i += 2
        If Mod($n,$i) = 0 Then Return 0
        $i += 4
    WEnd
    Return 1
EndFunc

Func _IsPrime3($n);optimized one
    If $n < 1 Then Return
    If $n = 2 Or $n = 3 Or $n = 5 Then Return 1
    If $n = 1 Or Mod($n,2) = 0 Or Mod($n,3) = 0 Or Mod($n,5) = 0 Then Return
    Local $i = 7, $l = Int(Sqrt($n))+1
    While $i < $l
        If Mod($n,$i) = 0 Then Return 0
        $i += 4
        If Mod($n,$i) = 0 Then Return 0
        $i += 2
        If Mod($n,$i) = 0 Then Return 0
        $i += 4
        If Mod($n,$i) = 0 Then Return 0
        $i += 2
        If Mod($n,$i) = 0 Then Return 0
        $i += 4
        If Mod($n,$i) = 0 Then Return 0
        $i += 6
        If Mod($n,$i) = 0 Then Return 0
        $i += 2
        If Mod($n,$i) = 0 Then Return 0
        $i += 6
    WEnd
    Return 1
EndFunc

j.

first you sort out all numbers dividable by 2,3 and 5. from 7 on, you test the sequence +4+2+4+2+4+6+2+6. the limit has to be int(sqrt($n))+1 so that e.g. 49 is covered by 7.

Edited October 22, 2008 by jennico

[YellowLab]

This looks extremely fascinating. Thanks for sharing.

Bob

[engine]

Hi jennico

Great work!

Why did you make a "_LogE" function? It's already implemented in AutoIt as "Log".

I did a _Factorial function in the last two days. It's here -> #593401

You can make yours simpler by removing one or two lines, this way:

Func _Factorial(Const $s_number)
    If Not IsInt($s_number) Or $s_number < 0 Then Return SetError(-1, 0, -1)
    Local $ret = 1
    For $i = 2 To $s_number
        $ret *= $i
    Next
    Return $ret
EndFunc   ;==>_Factorial

You don't need to declare $i variable Local. As it is inside a For, To, Next loop.

Unfortunately AutoIt can only calculate factorial up to _Factorial(20).

Kind regards.

[jennico]

thank you all. i will update the codes later today.

i found an even faster _IsPrime (two times faster than above):

edit: no, this was false !

it seems that for..next is faster than while...wend. for what reason ever, but it's interesting.

didn't know i don't have to declare loop variants. (are you sure on that ?)

i made a loge function to make clear the difference. i think the built in log function is a little bit - hmm, unconventional. for example, the windows calculator does not use loge for log.

j.

Edited October 22, 2008 by jennico

[engine]

(...)

didn't know i don't have to declare loop variants. (are you sure on that ?)

(...)

It's stated in AutoIt's documentation:

For...To...Step...Next

(...)

Remarks

The Variable will be created automatically with a LOCAL scope, even when MustDeclareVars is on.

(...)

Also, you can make a little script to test for that.

Regards.

[jennico]

@malkey : thanks for your idea. the function names are exactly what i need (pow, fmod taken from c, i think). i will see if this corresponds to what i need. i hope so.

j.

edit: where did you find a connection between fmod / pow and msvcrt.dll ? i only find some trigonometric functions in it.

Edited October 22, 2008 by jennico

[jennico]

okay, i found something:

pow pow(x as double, y as double) as double math.bi Returns x to the power y.

fmod fmod(x as double, y as double) as double math.bi Calculates the remainder of x divided by y.

but this is not what i need. i could use fmod instead of mod() in order to get an double precision output. but i cannot use pow because Autoit cannot handle the 64bit result.

(in fact, that is why your attempts in your functions failed.)

what i need is

mulmod = mod(x*y,z) (modular multiplication) (Montgomery reduction) with double precision and

modex (or powmod) = mod(x^y,z) (modular exponentiation) with double precision

i am afraid, this is not implemented in msvcrt.dll. but it seems a good idea to search the dlls for it.

j.

Edited October 22, 2008 by jennico

[Malkey]

@malkey : thanks for your idea. the function names are exactly what i need (pow, fmod taken from c, i think). i will see if this corresponds to what i need. i hope so.

j.

edit: where did you find a connection between fmod / pow and msvcrt.dll ? i only find some trigonometric functions in it.

I google searched msvcrt.dll fmod,, and found

http://msdn.microsoft.com/en-us/library/aa246388(VS.60).aspx

qsort is another function from msvcrt.dll in the list of functions that has already been used with AutoIt. See here

http://www.autoitscript.com/forum/index.ph...st&p=474072

Edited October 22, 2008 by Malkey