Resurrecting a 2000 Year Old Thread

[czardas]

Some time ago a question came up in General Help and Support about primes. I've always had a fascination with this apparently chaotic number sequence and thought this is a worthwhile challenge. Without reading up on the subject, I set about trying to figure out how I might approach the problem of testing for primality using AutoIt.
 
I started by writing a function which eliminated even numbers, and multiples of 3 and 5. Everything seemed fine at this point; but when I tried to scale up my idea, I took an interesting wrong turn. The approach I had adopted involved testing the primes in the range 7 to 31, and then incrementing them by steps of 30. I had calculated that this would produce a set of numbers guaranteed to include all higher primes. Testing for division using this smaller set of numbers would be faster than testing all the odd numbers ending 1, 3, 7 and 9.

At this stage the only thing that seemed obvious to me was - the more prime multiples I skip, the faster my function would run. Here's what I was thinking at the time: to scale this up I need to eliminate a few more prime multiples from a larger sequence and increment by much larger steps. All the numbers in the larger sequences will be made up exclusively of all the higher primes. BUT NOT WITH MY METHOD! This was a fatal oversight on my part.

For a while I was distracted away from coding and this bug went undetected. After a slightly embarrassing moment in the Chat forum I discovered that some prime factors were missing from my sequence ==> 3673. My error was very basic - I had made an assumption without any proof.
 

1 * 3 * 5 * 7 * 9 ** 11 ** 13 ** 15 ** 17 ** 19 ** 21 ** 23 ** 25 ** 27 ** 29 ** 31 Odd numbers
1 2 * 4 5 * 7 8 * 10 11 ** 13 14 ** 16 17 ** 19 20 ** 22 23 ** 25 26 ** 28 29 ** ** Factors of 3
1 2 3 4 * 6 7 8 9 ** 11 12 13 14 ** 16 17 18 19 ** 21 22 23 24 ** 26 27 28 29 ** 31 Factors of 5
1 _________ 7 ______ 11 __ 13 ________ 17 __ 19 ________ 23 ______________ 29 __ 31 Sequence

;

The sequence created above consists exclusively of prime numbers but the larger sequence generated below contains five extra numbers which are not prime - 121, 143, 169, 187 and 209.;

1 * 3 * 5 * 7 * 9 ** 11 ** 13 ** 15 ** 17 ** 19 ** 21 ** 23 ** 25 ** 27 ** 29 ** 31 ** 33 ** 35 ** 37 ** 39 ** 41 ** 43 ** 45 ** 47 ** 49 ** 51 ** 53 ** 55 ** 57 ** 59 ** 61 ** 63 ** 65 ** 67 ** 69 ** 71 ** 73 ** 75 ** 77 ** 79 ** 81 ** 83 ** 85 ** 87 ** 89 ** 91 ** 93 ** 95 ** 97 ** 99 *** 101 *** 103 *** 105 *** 107 *** 109 *** 111 *** 113 *** 115 *** 117 *** 119 *** 121 *** 123 *** 125 *** 127 *** 129 *** 131 *** 133 *** 135 *** 137 *** 139 *** 141 *** 143 *** 145 *** 147 *** 149 *** 151 *** 153 *** 155 *** 157 *** 159 *** 161 *** 163 *** 165 *** 167 *** 169 *** 171 *** 173 *** 175 *** 177 *** 179 *** 181 *** 183 *** 185 *** 187 *** 189 *** 191 *** 193 *** 195 *** 197 *** 199 *** 201 *** 203 *** 205 *** 207 *** 209 *** 
1 2 * 4 5 * 7 8 * 10 11 ** 13 14 ** 16 17 ** 19 20 ** 22 23 ** 25 26 ** 28 29 ** 31 32 ** 34 35 ** 37 38 ** 40 41 ** 43 44 ** 46 47 ** 49 50 ** 52 53 ** 55 56 ** 58 59 ** 61 62 ** 64 65 ** 67 68 ** 70 71 ** 73 74 ** 76 77 ** 79 80 ** 82 83 ** 85 86 ** 88 89 ** 91 92 ** 94 95 ** 97 98 ** 100 101 *** 103 104 *** 106 107 *** 109 110 *** 112 113 *** 115 116 *** 118 119 *** 121 122 *** 124 125 *** 127 128 *** 130 131 *** 133 134 *** 136 137 *** 139 140 *** 142 143 *** 145 146 *** 148 149 *** 151 152 *** 154 155 *** 157 158 *** 160 161 *** 163 164 *** 166 167 *** 169 170 *** 172 173 *** 175 176 *** 178 179 *** 181 182 *** 184 185 *** 187 188 *** 190 191 *** 193 194 *** 196 197 *** 199 200 *** 202 203 *** 205 206 *** 208 209 *** 
1 2 3 4 * 6 7 8 9 ** 11 12 13 14 ** 16 17 18 19 ** 21 22 23 24 ** 26 27 28 29 ** 31 32 33 34 ** 36 37 38 39 ** 41 42 43 44 ** 46 47 48 49 ** 51 52 53 54 ** 56 57 58 59 ** 61 62 63 64 ** 66 67 68 69 ** 71 72 73 74 ** 76 77 78 79 ** 81 82 83 84 ** 86 87 88 89 ** 91 92 93 94 ** 96 97 98 99 *** 101 102 103 104 *** 106 107 108 109 *** 111 112 113 114 *** 116 117 118 119 *** 121 122 123 124 *** 126 127 128 129 *** 131 132 133 134 *** 136 137 138 139 *** 141 142 143 144 *** 146 147 148 149 *** 151 152 153 154 *** 156 157 158 159 *** 161 162 163 164 *** 166 167 168 169 *** 171 172 173 174 *** 176 177 178 179 *** 181 182 183 184 *** 186 187 188 189 *** 191 192 193 194 *** 196 197 198 199 *** 201 202 203 204 *** 206 207 208 209 *** 
1 2 3 4 5 6 * 8 9 10 11 12 13 ** 15 16 17 18 19 20 ** 22 23 24 25 26 27 ** 29 30 31 32 33 34 ** 36 37 38 39 40 41 ** 43 44 45 46 47 48 ** 50 51 52 53 54 55 ** 57 58 59 60 61 62 ** 64 65 66 67 68 69 ** 71 72 73 74 75 76 ** 78 79 80 81 82 83 ** 85 86 87 88 89 90 ** 92 93 94 95 96 97 ** 99 100 101 102 103 104 *** 106 107 108 109 110 111 *** 113 114 115 116 117 118 *** 120 121 122 123 124 125 *** 127 128 129 130 131 132 *** 134 135 136 137 138 139 *** 141 142 143 144 145 146 *** 148 149 150 151 152 153 *** 155 156 157 158 159 160 *** 162 163 164 165 166 167 *** 169 170 171 172 173 174 *** 176 177 178 179 180 181 *** 183 184 185 186 187 188 *** 190 191 192 193 194 195 *** 197 198 199 200 201 202 *** 204 205 206 207 208 209 *** 
1 __________________ 11 __ 13 ________ 17 __ 19 ________ 23 ______________ 29 __ 31 ______________ 37 ________ 41 __ 43 ________ 47 ______________ 53 ______________ 59 __ 61 ______________ 67 ________ 71 __ 73 ______________ 79 ________ 83 ______________ 89 ____________________ 97 _________ 101 ___ 103 ____________107 ___ 109 ___________ 113 ___________________________ 121 ___________________ 127 ___________ 131 ___________________ 137 ___ 139 ___________ 143 ___________________ 149 ___ 151 ___________________ 157 ___________________ 163 ___________ 167 ___ 169 ___________ 173 ___________________ 179 ___ 181 ___________________ 187 ___________ 191 ___ 193 ___________ 197 ___ 199 ___________________________________ 209 ___

;

So there was a big flaw in my reasoning. My earlier assumption that the sequence in my function had to be comprised only of prime numbers was totally wrong. I imagine this may be the cause of amusement to some, but it also goes to show how easy it can be to overlook something fundamental like this. Always check your assumptions are correct and save yourself any future embarrassment!
 

Currently untested on an infinite series.

ConsoleWrite(_IsPrime(999999937) & @LF & _IsPrime(999999961) & @LF)

Func _IsPrime($iNumber)
    If Not IsInt($iNumber) Then Return SetError(1, 0, False) ; Only integers are allowed

    $iNumber = Abs($iNumber)
    If $iNumber = 2 Or $iNumber = 3 Or $iNumber = 5 Or $iNumber = 7 Then Return True
    If $iNumber < 2 Or Not (Mod($iNumber, 2) And Mod($iNumber, 3) And Mod($iNumber, 5) And Mod($iNumber, 7)) Then Return False

    Local $nRoot = Sqrt($iNumber)
    If IsInt($nRoot) Then Return False
    $nRoot = Floor($nRoot)

    Static $aFactors[48] = [11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,101,103,107,109,113,121,127,131,137,139,143,149,151,157,163,167,169,173,179,181,187,191,193,197,199,209,211]

    Local $bPrime = True
    For $i = 0 To 47
        For $j = $aFactors[$i] To $nRoot Step 210 ; = 2 *3 *5 *7
            If Mod($iNumber, $j) Then ContinueLoop ; $j is not a factor
            $bPrime = False
            ExitLoop 2
        Next
    Next

    Return $bPrime
EndFunc ; _IsPrime

Edited January 31, 2014 by czardas

[jchd]

Welcome to the wonderland of arithmetic. Primes, their repartition, their fantastic properties have always fascinated countless people, mathematicians and laymen alltogether.

Simple as that may seem at first look, characterizing them in general is often far from trivial and many pitfalls are awaiting the unsuspecting amateur.

Every primality test, compositeness test, factoring program, younameit will start with a trial division step to remove small factors. As you have experienced, a loop step of 30 without further checking is not enough to guarantee scaling. Don't by shy: there's nothing wrong by being fooled by intuition. It just needs to be backed by a more in-depth analysis.

BTW there is some relation between primes and guitar.

[czardas]

What an excellent and appropriate article jchd. This guy has obviously spent time studying the guitar and knows how demanding it can be. I very much appreciate your responce.

I posted my little story above in the hope that it may serve as an illustration to some younger readers who may be just starting out - I consider it to be an instructive mistake. By the time I sat down and wrote out the sequences, I already knew what I was looking for because I had eliminated all other possible causes for this error that I could think of. This was a valuable lesson - for me at least.

Now I should get back to studying my tremolo!

Edited January 31, 2014 by czardas

[jchd]

That or program a prime nuimber generator based on Jones, Wada & al. polynomial:

P(a,...,z) =
    (k+2) (1 - 
    [wz + h + j - q]^2 - 
    [(gk + 2g + k + 1)(h + j) + h - z]^2 - 
    [16(k + 1)^3(k + 2)(n + 1)^2 + 1 - f^2]^2 - 
    [2n + p + q + z - e]^2 - 
    [e^3(e + 2)(a + 1)^2 + 1 - o^2]^2 - 
    [(a^2 - 1)y^2 + 1 - x^2]^2 - 
    [16r^2y^4(a^2 - 1) + 1 - u^2]^2 - 
    [n + l + v - y]^2 - 
    [(a^2 - 1)l^2 + 1 - m^2]^2 - 
    [ai + k + 1 - l - i]^2 - 
    [((a + u^2(u^2 - a))^2 - 1)(n + 4dy)^2 + 1 - (x + cu)^2]^2 - 
    [p + l(a - n - 1) + b(2an + 2a - n^2 - 2n - 2) - m]^2 - 
    [q + y(a - p - 1) + s(2ap + 2a - p^2 - 2p - 2) - x]^2 - 
    [z + pl(a - p) + t(2ap - p^2 - 1) - pm]^2)

Plug any non-negative values into variables a...z and when the result is positive, then it is a prime. Surprisingly enough, you can generate all primes this way, but in totally unpractical order.

More here.

Edited January 31, 2014 by jchd

[czardas]

OMG - that looks really tricky.

[jchd]

Less difficult than finding it in the first place and actually proving that the set of positive results it generates is exactly the set of primes.

[czardas]

Less difficult than finding it in the first place and actually proving that the set of positive results it generates is exactly the set of primes.

 

Of course. Of the few basic (Euclidean) proofs I have read and managed to understand, I have always been struck by the beauty and rigor they contain. It is often an inspiration when the truth dawns on you.

[jchd]

A simple (sometimes not so) introduction which you can read in your spare time.

[czardas]

It's good to have this info in one place, however there is some rather advanced stuff in some of the links posted here. Learning maths is a great asset to any would-be coder, but you don't have to be an expert in the subject to proceed. It is possible to understand how I made my mistake (described above) without necessarily understanding all the details. This is still of value to anyone learning to write code.

Knowing what pitfalls to avoid is just as relevant as learning how to write good code. There are some fantastic examples in this forum to follow, and we all learn things in different ways. Some people tend to visualize things more, and a group of people have photographic memories. I was led down the wrong path by false intuition. It was still a step towards my ultimate goal - albeit one which fell short of the target in the first instance.

Edited February 1, 2014 by czardas

[czardas]

I have taken this a little further, but quickly hit limits with optimization as I expected. Well there are other methods as jchd has shown, but if you use any of these functions, make sure to read the descriptions carefully. I have added a couple of examples at the end of the code and a new function to return a random prime of a specified number of digits (between 1 to 15). Searching large ranges will be slow and greedy on RAM: use at your own risk. Don't forget to clear memory as shown in the examples.

;

expandcollapse popup

; #FUNCTION# ====================================================================================================================
; Name...........: _ClearDivisors
; Description ...: Used to clear memory after calling _IsPrime() or _RandomPrime()
; Syntax.........: _ClearDivisors()
; Author ........: czardas
; Comments ......; Call this function if you no longer need to use _IsPrime() or _RandomPrime() and want to clear memory.
; ===============================================================================================================================

Func _ClearDivisors()
    _IsPrime("", True)
EndFunc ;==> _ClearDivisors


; #FUNCTION# ====================================================================================================================
; Name...........: _IsPrime
; Description ...: Tests a number for primality.
; Syntax.........: _IsPrime($iInt [, _ClearDivisors = False])
; Parameters ....: $iInt           - The number to test
;                : $bClearDivisors - DO NOT USE THIS PARAMETER
; Return values .: Success   - Returns True, or False if not prime.
;                  Failure   - Returns False and sets @error to 1 if $iInt is not an integer or is out of range:
; Author ........: czardas
; Comments ......; This function is designed to be called multiple times but uses lots of memory.
;                ; It is advisable to use Sleep() after any batch processes to allow the CPU to cool off a little.
;                ; After using this function you should free memory by calling the function _ClearDivisors()
; ===============================================================================================================================

Func _IsPrime($iInt, $bClearDivisors = False)

    Static $aDivisors = __GetDivisors()

    If $bClearDivisors Then
        $aDivisors = 0
        Return
    ElseIf Not IsArray($aDivisors) Then
        $aDivisors = __GetDivisors()
    EndIf

    If Not IsInt($iInt) Then Return SetError(1, 0, False) ; Only integers are allowed
    $iInt = Abs($iInt)

    If $iInt = 2 Or $iInt = 3 Or $iInt = 5 Or $iInt = 7 Or $iInt = 11 Or $iInt = 13 Then Return True
    If $iInt < 2 Or Not (Mod($iInt, 2) And Mod($iInt, 3) And Mod($iInt, 5) And Mod($iInt, 7) And Mod($iInt, 11) And Mod($iInt, 13)) Then Return False

    Local $bPrime = True, $iRoot = Sqrt($iInt)
    If IsInt($iRoot) Then Return False

    $iRoot = Floor($iRoot)
    For $i = 0 To 5759
        For $j = $aDivisors[$i] To $iRoot Step 30030 ; 2 *3 *5 *7 *11 *13
            If Mod($iInt, $j) Then ContinueLoop ; $j is not a factor
            $bPrime = False
            ExitLoop 2
        Next
    Next

    Return $bPrime
EndFunc ;==> _IsPrime


; #FUNCTION# ====================================================================================================================
; Name...........: _RandomPrime
; Description ...: Returns a random prime number of between 1 and 15 digits in length (length must be specified).
; Syntax.........: _RandomPrime($iDigits)
; Parameters ....: $iDigits - The number of digits in the random prime number
; Return values .: Success   - Returns a random prime number with the required number of digits.
;                  Failure   - Sets @error to 1 if $iDigits is not an integer or is out of range:
; Author ........: czardas
; Comments ......; This function is designed to be called multiple times but uses lots of memory.
;                ; It is advisable to use Sleep() after any batch processes to allow the CPU to cool off a little.
;                ; After using this function you should free memory by calling the function _ClearDivisors()
; ===============================================================================================================================

Func _RandomPrime($iDigits)
    If Not IsInt($iDigits) Or $iDigits < 1 Or $iDigits > 15 Then Return SetError(1)
    $iDigits -= 1

    Local $iNumber = Number(Random(1, 9, 1) & StringMid(Random(), 3, $iDigits))
    Local $iStep = 1 -2*Random(0, 1, 1)

    Local $aRange[15][2] = [[2,7],[11,97],[101,997],[1009,9973],[10007,99991],[100003,999983], _
    [1000003,9999991],[10000019,99999989],[100000007,999999937],[1000000007,9999999967], _
    [10000000019,99999999977],[100000000003,999999999989],[1000000000039,9999999999971], _
    [10000000000037,99999999999973],[100000000000031,999999999999989]]

    If $iNumber > $aRange[$iDigits][1] And $iStep = 1 Then
        $iNumber = $aRange[$iDigits][0]
    ElseIf $iNumber < $aRange[$iDigits][0] And $iStep = -1 Then
        $iNumber = $aRange[$iDigits][1]
    EndIf

    While 1 ; Do not increase the size of the next loop.
        For $iNumber = $iNumber To $iNumber + 10*$iStep Step $iStep
            If _IsPrime($iNumber) Then ExitLoop 2
        Next
        Sleep($iDigits + 16) ; CPU needs to cool down.
    WEnd
    Return $iNumber
EndFunc ;==> _RandomPrime


; #INTERNAL_USE_ONLY# ===========================================================================================================
; Name...........: __GetDivisors
; Description ...: Sieve to remove prime multiples from the number range 17 to 30045.
; Syntax.........: __GetDivisors()
; Return values .: Success   - Returns an array of the remaining numbers within the range.
; Author ........: czardas
; Comments ......; This function is called internally by _IsPrime
; ===============================================================================================================================

Func __GetDivisors()
    Local $aPrimes[5] = [3,5,7,11,13], $aNumbers[15015], $iTemp
    For $i = 0 To 15014
        $iTemp = $i*2 +17 ; Produces odd numbers from 17 to 30045
        For $j = 0 To 4
            If Not Mod($iTemp, $aPrimes[$j]) Then ContinueLoop 2 ; Skip multiples of 3, 5, 7, 11 and 13
        Next
        $aNumbers[$i] = $iTemp
    Next
    $iTemp = 0

    Local $aDivisors[5760] ; Will only contain the factors that made it through the sieve.
    For $i = 0 To 15014
        If $aNumbers[$i] Then
            $aDivisors[$iTemp] = $aNumbers[$i]
            $iTemp += 1
        EndIf
    Next
    Return $aDivisors
EndFunc ;==> __GetDivisors


; ****************************************************************************************
; EXAMPLE 1

Global $iMax = 0

ConsoleWrite("Highest primes from 1 to 12 digits" & @LF)

For $i = 1 To 12
    $iMax &= 9 ; $iMax now becomes a string.
    For $j = $iMax to 10^($i -1) Step -2 ; $iMax is still interpreted as a number here.
        If _IsPrime($j) Then
            ConsoleWrite($j & @LF)
            Sleep(50) ; The larger the search range, the more time the CPU needs to cool off.
            ExitLoop
        EndIf
    Next
Next

_ClearDivisors() ; Free up memory

; ****************************************************************************************
; EXAMPLE 2

ConsoleWrite(@LF & "Random primes from 1 to 15 digits" & @LF)

For $i = 1 to 15
    ConsoleWrite(_RandomPrime($i) & @LF)
    Sleep(10)
Next

_ClearDivisors() ; Free up memory

Edited February 4, 2014 by czardas

[jchd]

Good job.

Not being picky, but I wanted to test: >my bignum implementation of _IsPrime (or its complement) favorably competes in both speed and memory usage over a large range of input values. We can even increase the number of iterations of Miller-Rabin to increase confidence on the result by a huge margin. 10 iterations provides probability 0.999999046325684, 15 iterations give 0.999999999068677, 20 iterations provide 0.999999999999091 and 25 iterations give 0.999999999999999

Edited February 3, 2014 by jchd

[czardas]

Thanks jchd - your words mean a lot to me. I don't have a great deal of use for my function personally, but it might be possible to make use of it with music composition somehow - I'm not sure. I don't need it for encryption purposes and I imagine (my) 15 digit version would be too weak for this. I do like the method of cycling through a sequence - it's similar to generating rhythms.