Counting Primes

[sneakyimp]

I've been reading Bruce Schneier's book Applied Cryptography, which is pretty awesome.

I'm interested in doing a little experimentation to help me better understand some of the mathematical concepts pertaining to trap door functions and key generation. In particular, I was thinking I might write a little PHP script to locate all the primes of bit length N to help me make some graphs so I understand a bit more about the distribution of prime numbers in a particular numeric range. Were I a better mathematician I could probably just take a look at Prime Number Theorem or perhaps the Riemann Zeta Function, but I figured I'd start small.

That said, anyone have suggestions for a good algorithm for locating all the Primes for an N-bit integer? I'm looking at the Sieve of Eratosthenes and the the Sieve of Atkin.

[dalecosp]

I'm not much of a mathematician; Weedpacket is quite a good one, though. Perhaps he'll chime in.

Of some interest, perhaps historical, probably in terms of discussion: [thread=10277645]A Prime Example[/thread]. Were you around back then? Doesn't look like you got in on the contest, though. 🙁

[Derokorian]

While I don't have a specific solution to what you're asking. I can share something I wrote before to determine the number of primes in a range. Taking that, you could determine a range for an integer specified by N-bits and pass that range in:

function rangePrimes($A, $B) {
    $intPrimes = 0;
    for( $i = $A; $i <= $B; $i++ ) {
        if( isPrime($i) )
            $intPrimes++;
    }
    return $intPrimes;
}

function isPrime($num) {
    //1 is not prime
    if($num == 1)
        return false;

//2 is prime 
if($num == 2)
    return true;

// even numbers other than 2 are not prime
if($num % 2 == 0) {
    return false;
}

/**
 * Checks the odd numbers. If any of them is a factor, then it returns false.
 * The sqrt can be an aproximation, hence just for the sake of
 * security, one rounds it to the next highest integer value.
 */
for($i = 3; $i <= ceil(sqrt($num)); $i = $i + 2) {
    if($num % $i == 0)
        return false;
}

return true;
}

[NogDog]

The GMP extension has some prime number functions that might help.

[johanafm]

sneakyimp;11036459 wrote:

That said, anyone have suggestions for a good algorithm for locating all the Primes for an N-bit integer?

.
In case you want to try finding a good algorithm yourself you only really need to know that
- any number N can have no more than one primefactor > floor(sqrt(N))
(which sets the upper limit for what numbers you need to test)
- any prime above 3 can be expressed as 6k + 1 or 6k - 1
(since you are trying for prime factors, only numbers expresseable on this form needs be tried)

Doing it this way is faster than Derokorians algorithm since it tries a lot fewer numbers. However, I didn't come up with that algorithm on my own, but rather read about it at projecteuler.net after solving their 7th problem - finding the 10001st prime - using an algorithm similar to that of Derokorian's. The solution presented at projecteuler after you solve the problem (in any way) goes through the first 200k numbers looking for primes in 0.85 s on my laptop, compared to 6.4 s for the algorithm presented here.

But they don't want people spreading solutions. So in case you'd grow tired of trying or don't care about trying, you could creatre an account at projecteuler, use Derokorian's algorithm as a basis for solving their 7th problem, post the solution and then get your hands on their suggested algorithm

[Derokorian]

johanafm - I really like that site! Great puzzles, got 10 of them done out of boredom tonight. FWIW, my solution in js runs on my little ultrabook in ~360ms whereas the provided solution runs in ~315ms not sure why you are seeing such a significant difference with PHP. I will have to test them head to head later. My ultrabook doesn't have PHP or any type of web server installed.

[johanafm]

Derokorian;11036493 wrote:

johanafm - I really like that site! Great puzzles, got 10 of them done out of boredom tonight.

Oh yes - it's fun! Just recently found it myself and I've only done the first 8. And I really like the PDFs you get access to after each problem. I'm usually good enough to understand maths while I lack the brilliance to come up with those "smart solutions". But I do enjoy reading about them.

Perhaps I should mention that I found projecteuler in a list of sites with programming problems. There are lots of other (non-mathematical type) problem types listed too should anyone ever need something to do 🙂

Derokorian;11036493 wrote:

~360ms whereas the provided solution runs in ~315ms not sure why you are seeing such a significant difference with PHP

Ah, there was one major difference in my implementation of the 6k1 and mod2 functions, which are my name for their and your algorithm respectively. In mod2 I had the calculation of ceil(sqrt($num)) in the for loop comparison, while this was done once for the 6k1 implementation. I.e. the difference between

# This
$done = doItOnce();
for ($i = 0; $i < $done;)

# And this
for ($i = 0; $i < doItOverAndOver();)

Edit: And I had retained the same difference using copy-paste in js too…

Testing 00 .. 200k - php through webserver & CLI - same results
Mod2: 1.3012518882751 s
6k1: 0.8339319229126 s
ratio mod2:6k1: 1.56

Somewhat surprised to find that processing is now by far quicker in js
Testing 00 .. 200k - js implementations of the same two functions - firefox
mod2: 0.166 s
6k +- 1: 0.155 s
ratio mod2k:6k1: 1.07

Testing 00 .. 200k - js implementations of the same two functions - chrome
mod2: 0.023 s
result: 17984
6k +- 1: 0.015 s
ratio mod2:6k1: 1.53

Surprised that 6k1 does not have the same relative relation to mod2 in ff and chrome. But now Chrome is in line with the php results.

[sneakyimp]

johanafm;11036489 wrote:

.
- any number N can have no more than one primefactor > floor(sqrt(N))
(which sets the upper limit for what numbers you need to test)
- any prime above 3 can be expressed as 6k + 1 or 6k - 1
(since you are trying for prime factors, only numbers expresseable on this form needs be tried)

Are we certain of these assertions? I'm fairly certain that there are infinitely many primes of the form 6k -1 or 6k + 1, but are you certain that ALL primes fit this form? I searched around and have seen some attempts at a proof that seem a bit fast and loose.

johanafm;11036489 wrote:

Doing it this way is faster than Derokorians algorithm since it tries a lot fewer numbers. However, I didn't come up with that algorithm on my own, but rather read about it at projecteuler.net after solving their 7th problem - finding the 10001st prime - using an algorithm similar to that of Derokorian's. The solution presented at projecteuler after you solve the problem (in any way) goes through the first 200k numbers looking for primes in 0.85 s on my laptop, compared to 6.4 s for the algorithm presented here.

That's definitely an impressive improvement. Yes I'm certain we can improve on Derkorian's algorithm -- and your suggestions are excellent if we can verify them. So far they seem spot on to me.

johanafm;11036489 wrote:

But they don't want people spreading solutions. So in case you'd grow tired of trying or don't care about trying, you could creatre an account at projecteuler, use Derokorian's algorithm as a basis for solving their 7th problem, post the solution and then get your hands on their suggested algorithm

I'll have to go over there and look. FYI, you can get huge prime numbers from wolframalpha.com just by typing in "1000000th prime number"

[sneakyimp]

NogDog;11036477 wrote:

The GMP extension has some prime number functions that might help.

Thanks for that. I believe the built-in PHP functions are pretty usually much faster than any we might write ourselves.

[sneakyimp]

BTW if anyone is interested about how many primes exist for LARGE key values, the second post here is pretty cool:
http://math.stackexchange.com/questions/263588/how-many-all-prime-numbers-p-with-length-of-bits-of-p-1024-bits

Can't verify it's veracity, but the values are easily calculated at http://wolframalpha.com.

There are apparently a LOT of prime numbers with a bit length of 4096: http://www.wolframalpha.com/input/?i=pi++%282^{4096%29+-+pi+}+%282^4095%29

[johanafm]

sneakyimp;11036507 wrote:

That's definitely an impressive improvement.

See my last post above. Unfortunately that was simply due to moving an expression out of the for conditional. It is still an improvement, but they are not that far apart. One evaluates every third number (2 out of 6). The other evaluates every other number (odd numbers).

sneakyimp;11036507 wrote:

Are we certain of these assertions?

The first part is that there is at most one primefactor > sqrt(N). If sqrt(N) is R and R is an integer, then RR = N and there is no factor > R. If sqrt(N) is r and r is not an integer, let ceil(r) be R, the smallest integer > than r. Then R r is > N and therefor any other factor(s) of N has to be < r. Therefor, if you have found no divisor in the range 2 .. floor(sqrt(N)), there can be no divisor in ceil(sqrt(N)) .. N.

sneakyimp;11036507 wrote:

are you certain that ALL primes fit this form?

For the second part, all integers can be written as 6k + i, for some integer k and i = -1, 0, 1, 2, 3, 4

1. 6k + 0

2. 6k + 1

3. 6k + 2

4. 6k + 3

5. 6k + 4

6. 6k + 5 => 6 (k+1) -1

Every time you reach 6k + 5 (every 6th integer), you increase k by 1 and reset i to -1

Looking at these expressions, you will find that
1. is divisible by 2 and 3
3. is divisible by 2
4. is divisible by 3
5. is divisible by 2

This only leaves cases 2. and 6. that we are not certain about. Thus, any prime will be expressible as 6k ± 1. But not all numbers expressible as 6k ± 1 are prime though.

[Weedpacket]

The above analysis can be taken further (e.g., all primes greater than 5 are of the form (2×3×5)n±{1,7,11,13}) - what it amounts to is precomputing the first iterations of Eratosthenes's sieve.