LTD Puzzle 4: Finding Prime Numbers

[damber]

This weeks puzzle is quite straight forward.

All we want you to do is find every prime number between 0 and 1,000,000. That's it.

A prime number is a number that can only be divided by itself and 1 to result in a whole number. E.g. 2, 3 and 5 are all prime, but 4 can be divided by 2 aswell as 4 and 1, so it is not.

The faster your code is, the better - so if you were thinking of looping through each number, that may not be the most efficient route...

For extra points, find the Mersenne primes too.

For even more extra points we challenge you to find all the primes between 0 and 10,000,000. And for those hardcore geeks - how about 0 to 2,147,483,647 (the maximum for a 32 bit integer).

Good luck...

[MrChaos]

Ugle C code but, it should work.

cLine number Off | Show | Select all

1. #include <math.h>
2. #include <stdio.h>
3.  
4. /* How many numbers to search */
5. #define MAX_SEARCH 10000000
6.  
7. /* room for all primes 0 -> 2^32-1 (I think)*/
8. int primes[1024*1024*128];
9.  
10. void main(int argc, char** argv) {
11.   unsigned int current = 3;
12.   int maxSearch;
13.   unsigned int currentPrime;
14.   unsigned int primeCount = 1;
15.   primes[0] = 2;
16.   int isPrime = 0;
17.   while(current<MAX_SEARCH) {
18.     isPrime = 1;
19.     /* only search to the square root of the value */
20.     maxSearch = (int)sqrt((double)current);
21. #ifdef _DEBUG
22.     printf("MaxSearch for %d is %d\n", current, maxSearch);
23. #endif
24.     /* search as long as the number might still be a prime and as long as the current prime we use for division is lower that the max we need to check*/
25.     for(currentPrime = 0 ; isPrime==1 && primes[currentPrime]<=maxSearch; currentPrime++) {
26. #ifdef _DEBUG
27.       printf("Checking if %d division with %d is 0\n", current, primes[currentPrime]);
28. #endif
29.       isPrime = current%primes[currentPrime]!=0;
30.     }
31.     if(isPrime) {
32. #ifdef _DEBUG
33.       printf("Adding prime %d at index %d\n", current, primeCount);
34. #endif
35.       primes[primeCount++] = current;
36.     }
37.     current+=2;
38.   }
39.   for(current=0;current<primeCount;current++) {
40.     printf("%d = %d\n", current, primes[current]);
41.   }
42. }

Code is hidden, SHOW

[damber]

Well done to MrChaos for getting the first (and quick) answer in !

Just in case you were all wondering what the highest prime number was, AlexCuse has already posted that answer here: Topic1983: The Highest Prime Number in the World

[AlexCuse]

Yes, thank you Mr. Chaos, was not expecting to see much C on here. A most pleasant surprise.

damber, the veracity of that claim has come into dispute. Maybe those crazy mathematicians are right and we'll never know the highest prime number

[Emtucifor]

The highest prime number is infinitely high!

There is only the highest prime number man has identified/proven so far.

[traingamer]

I found a spiffy, easily modified online JavaScript example based on the Sieve of Erastoshenes, but I haven't gotten around to modifying it in VBA (or any other language). Maybe tonight.

[SQLDenis]

I know at least 2 people who have a Sieve of Erastoshenes SQL version

[damber]

The highest prime number is infinitely high!

Erik,
You should probably check out the link, and I think you'll agree that that's the highest prime number in the world... Denis tried to offer another possible highest prime, but lost out on a technicality...

traingamer, SQLDenis, - only original code please - otherwise where's the fun/challenge ?

[Emtucifor]

Should we get more/less points for not using/using someone else's code as a starting place?

Programming the puzzle from scratch is a different task than porting another's code, though if you get stuck I'd say that looking at other code is legitimate...

[SQLDenis]

traingamer, SQLDenis, - only original code please - otherwise where's the fun/challenge ?

I will not post the code, hopefully they will

[Emtucifor]

Here's an HTML/javascript version (yes, I know, funny). It doesn't perform too well. Even in Firefox you have to tell it to continue 3 times to calculate just 1 million of them.

I like this code because I didn't give it any seed numbers at all, not even the sieve list. The sieve size is not hardcoded, though this code can only handle sieve sizes of 1-3 because 4 comes out to 48 members in the sieve list and that's too many bits for storing the value as a long integer. I wrote this all from scratch.

I'd like to rewrite this code using C++ or even assembler, with some libraries for doing math on very large numbers. It would be cool to do a sort of sliding sieve: storage becomes a real problem when dealing with large quantities of primes, and with a sieve you get a fairly decent storage algorithm with 1 bit per candidate, though of course it becomes very sparse when you get to high numbers. So the sliding sieve would stairstep into the next size as it went along, conserving storage space at the expense of some CPU, but with a real tradeoff of checking fewer numbers along the way.

There are other optimizations besides the sieve that could be done. I'd love to have the time and energy to learn more about them and write code that actually uses them.

This program has a bug in that it may display a few more numbers than you asked for (though it calculates the primes correctly). I don't have time to fix it right now...

javascriptLine number Off | Show | Select all

1. <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01//EN" "http://www.w3.org/TR/html4/strict.dtd">
2. <html>
3. <head>
4. <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
5. <title>Primes</title>
6. </head>
7. <body>
8. <script type="text/javascript">
9. var
10.     sieveprimecount = 3,
11.     calclimit = 1000000,
12.     sieveprimes, sievesize,  sieve = [],  sievecount, result = [];
13.  
14. function arrmult(arr, start, finish) {
15.     start = start === undefined ? 0 : start;
16.     finish = finish === undefined ? arr.length - 1 : finish;
17.     var result = 1;
18.     while (start<=finish){result*=arr[start++];}
19.     return result;
20. }
21.  
22. function isDiv(num, testarr, startInd, endInd) {
23.     while(startInd<=endInd){
24.         if(num%testarr[startInd]===0){return true;}
25.         startInd++;
26.     }
27.     return false;
28. }
29.  
30. function buildArr(num, limit, testArr, startInd, endInd, add, seed) {
31.     seed=seed===undefined ? [] : seed;
32.     while((num <= limit||limit===undefined)&&(testArr.length<=endInd||endInd===undefined)) {
33.         if (!isDiv(num, testArr, startInd, endInd===undefined?testArr.length-1:endInd)){
34.             seed[seed.length]=num;
35.             if(add){testArr[testArr.length]=num;}
36.         }
37.         num++;
38.     }
39.     return seed;
40. }
41.  
42. function arrIndexOf(testArr, num){
43.     var i = testArr.length;
44.     while(i--){
45.         if(testArr[i]===num){return i;}
46.     }
47.     return -1;
48. }
49.  
50. function arrToBits(inputArr, compArr){
51.     var bits = 0,
52.         i=inputArr.length,
53.         pos;
54.     while(i--){
55.         pos = arrIndexOf(compArr, inputArr[i]);
56.         if(pos>=0){bits |= (1 << pos);}
57.     }
58.     return bits;
59. }
60.  
61. function testPrime(num, limit) {
62.     var chunk=-1,
63.         i;
64.     while(++chunk<=limit){
65.         for(i=0;i<sievecount;i++){
66.             if (((result[chunk]>>i) & 1===1) && (num%(chunk*sievesize+sieve[i])===0)) {return false;}
67.         }
68.     }
69.     return true;
70. }
71.  
72. function displayPrimes() {
73.     var i, limit, M = 32;
74.     for(var pos=sievesize;pos <= calclimit; pos+=sievesize){
75.         result[pos / sievesize] = 0;
76.         limit = Math.floor(Math.floor(Math.pow(pos+sievesize-1, 0.5)) / sievesize);
77.         for(i=0;i<sievecount && (pos + i <= calclimit);i++){
78.             if(testPrime(pos+sieve[i], limit)){
79.                 result[pos / sievesize] |= (1 << i);
80.                 document.write(', ' + (pos+sieve[i]));
81.                 if(pos+sieve[i]+1===M) {document.write(' (M)');}
82.                 if(pos+sieve[i]>M){M *= 2;}
83.             }
84.         }
85.     }
86. }
87.  
88. function displayResult(){
89.     displayArr(sieveprimes)
90.     var arr=[];
91.     for(var i=1;i<sievesize;i++){
92.         if((result[0]>>i) & 1 ===1) {
93.             arr[arr.length] = sieve[i];
94.         }
95.     }
96.     document.write(', ');
97.     displayArr(arr);
98.     displayPrimes();
99. }
100.  
101. function displayArr(arr) {
102.     var M;
103.     for(var i=0;i<arr.length;i++){
104.         M = arr[i] + 1;
105.         while (M > 2) {
106.             if (M%2===0) {M=M/2} else {M = 0}
107.         }
108.         if(M===2){arr[i]+=' (M)'};
109.     }
110.     document.write(arr.join(', '));
111. }
112.  
113. sieveprimes = buildArr(2, undefined, [], 0, sieveprimecount - 1, true);
114. sievesize = arrmult(sieveprimes);
115. sieve = [1].concat(buildArr(sieveprimes[sieveprimes.length - 1] + 1, sievesize, sieveprimes, 0, undefined, false));
116. sievecount = sieve.length;
117. result = [arrToBits(buildArr(sieveprimes[sieveprimes.length - 1] + 1, sievesize, sieveprimes.slice, 0, undefined, true), sieve)];
118.  
119. displayResult();
120. </script>
121. </body>
122. </html>

Code is hidden, SHOW

P.S. Damber, I have an idea for a puzzle some week: create your own functions to add, subtract, divide, multiply, and calculate power and log for numbers up to 1000 digits (or some other huge value).

[spoulson]

Ugle C code but, it should work.

I like the code, but unfortunately it's not completely copy/paste. Visual C++ 2008 complains of a couple points. It didn't like the huge static size on primes[] and threw a stack overflow. I replaced it with a pointer and a malloc() call. Also, in C, you can't intermix declaration lines with code, so the "primes[0] = 2" had to be moved down below the decls.

Overall, it works great! Thanks for the submission.

[traingamer]

When I said modifying it, I meant the sieve, not the JavaScript code.

[Emtucifor]

Oops... I hardcoded the Mersenne prime seed number. A fix:

javascriptLine number Off | Show | Select all

1. function displayPrimes() {
2.    var i, limit, M = 2;
3.    while (M<=sievesize){M*=2};
4.    ...

Code is hidden, SHOW

[AlexCuse]

I had WAY too much fun with this one. I must have tried at least 10 different ways of doing it in F#. Here's what I came up with in the end. While its slower than what I had in C#, I find it aesthetically pleasing. Finds all primes up to 2147483616 (the max capacity of the data structure I chose, even though it is initialized with an int32!) in about ten minutes:

fsharpLine number Off | Show | Select all

1. #light
2. open System
3. open System.Collections
4. open System.Collections.Generic
5. let start = DateTime.Now
6. let capacity = 2147483616
7. let isComposite = new BitArray(capacity)
8. isComposite.Set(0, true)
9. isComposite.Set(1, true)
10.  
11. let mutable i = 2
12.  
13. let rec next n =
14.     if not (isComposite.Get n) or (n=capacity) then
15.         n
16.     else next (n+1)
17.    
18. while i<= int (Math.Sqrt((double)capacity)) do
19.     for q in (i + i) .. i .. capacity do
20.         if q<capacity then
21.             isComposite.Set(q, true)
22.            
23.     i <- next (i + 1)
24.  
25. let primes = new List<int>()
26.  
27. for p in 2 .. 1 .. (capacity - 1) do
28.     if not (isComposite.Get p) then
29.         primes.Add p
30.  
31. printfn "Primes found is %d"primes.Count
32. Console.WriteLine (DateTime.Now - start)
33. Console.ReadKey(true)

Code is hidden, SHOW

However, I'm not sure if it is because of my lack of skill in F#, or because the F# compiler doesn't have as many years of wisdom as C#'s, but I was able to do much better in C# using the same method (minus the recursion, which seemed to actually make it slower). I had a really cool solution that returned an IEnumerable, but this seemed to get bogged down in the higher numbers. Using a similar method to the one used in F# could find all the primes to 2147483616 in 4:27 on my laptop (with everything coded in main), but I wanted to find as high as I could. This led me to some guy's own implementation of BitArray (http://www.jaggersoft.com/csharp_standard/17.8.htm), which does offer the size advertised. Its' a bit slower, but I thought it beat writing a ton of extra code to work with an "overflow" BitArray.

So, here's the BitArray2 class:

csharpLine number Off | Show | Select all

1. using System;
2. using System.Collections.Generic;
3. using System.Linq;
4. using System.Text;
5.  
6. namespace primes
7. {
8.     public class BitArray2
9.     {
10.         int[] bits;
11.         int length;
12.  
13.         public BitArray2(int length)
14.         {
15.             if (length < 0) throw new ArgumentException();
16.             bits = new int[((length - 1) >> 5) + 1];
17.             this.length = length;
18.         }
19.  
20.         public int Length
21.         {
22.             get { return length; }
23.         }
24.  
25.         public bool this[int index]
26.         {
27.             get
28.             {
29.                 if (index < 0 || index >= length)
30.                 {
31.                     throw new IndexOutOfRangeException();
32.                 }
33.                 return (bits[index >> 5] & 1 << index) != 0;
34.             }
35.             set
36.             {
37.                 if (index < 0 || index >= length)
38.                 {
39.                     throw new IndexOutOfRangeException();
40.                 }
41.                 if (value)
42.                 {
43.                     bits[index >> 5] |= 1 << index;
44.                 }
45.                 else
46.                 {
47.                     bits[index >> 5] &= ~(1 << index);
48.                 }
49.             }
50.         }
51.     }
52. }

Code is hidden, SHOW

And the "sledgehammer" method of getting the primes (I added int32.maxvalue at the end only if that's the capacity given, because we know that's prime and it is JUST out of reach). I'm not sure if the try/catch works any faster than checking the values prior to assignment, but I would imagine it does (i'm out of time for testing). This saves us from having to check for values <0 in the outer loop as well (since if you try to add past int32.maxvalue you get a negative #).

Time for this (with checking the values twice on each loop iteration) was 5:46. Hoping to have a test with the try/catch done running by the time I type the rest of this.

csharpLine number Off | Show | Select all

1. using System;
2. using System.Collections.Generic;
3. using System.Linq;
4. using System.Text;
5.  
6. namespace primes
7. {
8.     class PrimeFinderNE
9.     {
10.         private static BitArray2 _isComposite;
11.         public static List<int> FindPrimes(int capacity)
12.         {
13.             _isComposite = new BitArray2(capacity);
14.  
15.             _isComposite[0] = true;
16.  
17.             _isComposite[1] = true;
18.  
19.             int i = 2;
20.  
21.             while (i <= Math.Sqrt(_isComposite.Length))
22.             {
23.                 for (int z = i + i; z < _isComposite.Length; z += i)
24.                 {
25.                     try
26.                     {
27.                         _isComposite[z] = true;
28.                     }
29.                     catch (IndexOutOfRangeException ex)
30.                     {
31.                         break;
32.                     }
33.                 }
34.  
35.                 for (int q = i + 1; q < _isComposite.Length; q++)
36.                 {
37.                     if (!_isComposite[q])
38.                     {
39.                         i = q;
40.                         break;
41.                     }
42.                 }
43.             }
44.  
45.             List<int> primes = new List<int>();
46.  
47.             for (int pi = 0; pi < _isComposite.Length; pi++)
48.                 if (!_isComposite[pi])
49.                     primes.Add(pi);
50.  
51.             if (capacity == int.MaxValue)
52.                 primes.Add(int.MaxValue);
53.  
54.             return primes;
55.         }
56.     }
57. }

Code is hidden, SHOW

Times for the C# methods include querying the list to get Min and Max values as well as count (which is no trivial task when you're working with such a long list). And, handling the exceptions was slightly faster than preventing them at 5:29

[Emtucifor]

And a fix for my "sometimes it displays more numbers than you asked for":

javascriptLine number Off | Show | Select all

1. function displayPrimes() {
2. ...
3.    for(i=0;i<sievecount && (pos + sieve[i] < calclimit);i++){
4. ...

Code is hidden, SHOW

I needed to use sieve[i] instead of i.

[damber]

P.S. Damber, I have an idea for a puzzle some week: create your own functions to add, subtract, divide, multiply, and calculate power and log for numbers up to 1000 digits (or some other huge value).

Thanks Erirk - I've added it to our list of puzzle questions - so look out for it over the next few weeks

[chrissie1]

There are 2 types of winner - a "People's Champion" (the one with the most rating points across all their posts) and the LTD Admin's Champion (selected with our special formula....honest...well, maybe it's just the ones that the admins thought were particularly good )

So...

People's Champion: MrChaos scoring 11 points

LTD Admin's Champion: AlexCuse for usinf F#.


Congratulations to all.. and good luck on next weeks puzzle

[tisodotsk]

My solution in PHP:

phpLine number Off | Show | Select all

1. <?php
2. //ltd4.php - LTD Puzzle 4: Finding Prime Numbers
3. //via Topic1938: LTD Puzzle 4: Finding Prime Numbers
4. //------------------------------------------------------------------------------
5. $primes[]=2;
6. echo 'Prime numbers less than milion: '.$primes[0];
7. for($i=3; $i<=1000000; $i+=2){
8.     $is_prime=TRUE;
9.     foreach($primes as $val){
10.         if(ceil($ratio=$i/$val) === floor($ratio)){
11.             $is_prime=FALSE;
12.             break(1);
13.         }
14.         if($val > sqrt($i)){break(1);}
15.     }
16.     if($is_prime){
17.         $primes[]=$i;
18.         echo ', '.$i;
19.         //is Mersenne prime number?
20.         if(preg_match('/^0*1+$/', decbin($i), $match)){$m[]=$i;}
21.     }
22. }
23. echo '<br>End, '.count($primes).' prime numbers founded<br>';
24. echo 'Mersenne prime numbers: '.implode(', ', $m);
25. //------------------------------------------------------------------------------
26. ?>

Code is hidden, SHOW

I use some optimizations.
Output:

textLine number Off | Hide | Select all

1. Prime numbers less than milion: 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, ...
2. ...
3. End, 78498 prime numbers founded
4. Mersenne prime numbers: 3, 7, 31, 127, 8191, 131071, 524287

[Hazar]

my simple C++ solution:

textLine number Off | Show | Select all

1. #include <iostream>
2. using namespace std;
3.  
4. int main(void){
5.     //feel free to calculate to your hearts content
6.     int limit = 10000001;
7.  
8.     int* numberArray = (int*)calloc(sizeof(int)*limit, sizeof(int));
9.  
10.     for(int i = 2; i < limit; i++){
11.         if( numberArray[i] == 0){
12.             cout << i << "\n";
13.  
14.             //cross multiples off the list
15.             for(int j = i; j <= limit; j += i){
16.                 numberArray[j] = 1;
17.             }
18.         }
19.     }  
20.  
21.     return 0;
22. }
23.  

Code is hidden, SHOW

[AlexCuse]

Excellent Hazar. In the future could you please use the "hidecode" tag on your answers (only in the puzzle forum of course!)

You've just given me an idea too.

[GothGargoyle]

Here's a version in Java. On my machine it determined the primes < 2,147,483,647 in 2 minutes 16 seconds.

It is based on the Sieve of Eratosthenes algorithm, and uses a BitSet for minimizing the amount of memory required.

javaLine number Off | Show | Select all

1. package lessthandot;
2.  
3. import java.util.BitSet;
4.  
5. public class PrimeSieveOfEratosthenes {
6.  
7.     public static void main(final String[] args) {
8.         PrimeSieveOfEratosthenes sieve = new PrimeSieveOfEratosthenes();
9.         if (args.length < 1) {
10.             System.out.println("Please supply a maximum number");
11.             return;
12.         }
13.         long start = System.nanoTime();
14.         int max = Integer.parseInt(args[0]);
15.         BitSet primes = sieve.getPrimes(max);
16.         System.out.println("Found " + (primes.length() - primes.cardinality()) + " primes < " + max);
17.         System.out.println("Calculation took " + (System.nanoTime() - start) + " nano seconds");
18.         sieve.printPrimes(primes);
19.     }
20.  
21.     public void printPrimes(final BitSet primes) {
22.         for (int i = 0; i < primes.length(); i++) {
23.             if (!primes.get(i)) {
24.                 System.out.println(i + 1);
25.             }
26.         }
27.     }
28.    
29.     public BitSet getPrimes(final int max) {
30.         BitSet primes = new BitSet(max);
31.        
32.         primes.set(0);
33.         int candidate = 2;
34.         while (candidate < Math.sqrt(max)) {
35.             for (int i = 2; i <= (max / candidate); i++) {
36.                 primes.set(candidate * i - 1);
37.             }
38.             candidate = primes.nextClearBit(candidate) + 1;
39.         }
40.        
41.         return primes;
42.     }
43.  
44. }
45.  

Code is hidden, SHOW

[AlexCuse]

That's a pretty slick one GothGargoyle. Very clean and easy to read.