## Prime Numbers

Suppose that you have to find all prime numbers till say some value N or let’s say sum of all prime numbers till N. How would you proceed. Let’s play around a bit. A grade 5 student would tell you that a prime has only two divisors 1 and itself and would code it like this :

void basic1(int N) { int i,j,k,c; long long s=0; for (i=1;i<=N;i++) { c=0; for (j=1;j<=i && c<3;j++) if (i%j==0) c++; if (c==2) s+=i; } cout<<s<<endl; }

A little smarter kid would tell you that a prime number has no divisor between 2 and N/2 and would code it like this :

void basic2(int N) { int i,j,k,F; long long s=0; for (i=2;i<=N;i++) { F=1; for (j=2;j<=i/2 && F==1;j++) if (i%j==0) F=0; if (F) s+=i; } cout<<s<<endl; }

An even smarter kid will tell you that a prime number has no divisor between 2 and its root and would code it like this :

void moderate1(int N) { vector<int> p(N/2,0); int i,j,k,F,c=0; long long s=0; for (i=2;i<=N;i++) { F=1; for (j=2;j*j<=i && F==1;j++) if (i%j==0) F=0; if (F) { p[c++]=i; s+=i; } } cout<<s<<endl; }

A good maths student can tell you that a prime number has no prime divisors between 2 and its root.

void moderate2(int N) { vector<int> p(N/2,0); p[0]=2; int i,j,F,c=1; long long s=2; for (i=3;i<=N;i+=2) { F=1; for (j=0; p[j]*p[j]<=i && F==1; j++) if (i%p[j] == 0) F=0; if (F) { p[c++]=i; s+=i; } } cout<<s<<endl; }

A good programmer will tell you that Sieve of Eratosthenes is the best way to find the list of prime numbers.

Its algorithm is as follows :

- Create a contiguous list of numbers from two to some highest number n.
- Strike out from the list all multiples of two (4, 6, 8 etc.).
- The list’s next number that has not been struck out is a prime number.
- Strike out from the list all multiples of the number you identified in the previous step.
- Repeat steps 3 and 4 until you reach a number that is greater than the square root of n (the highest number in the list).
- All the remaining numbers in the list are prime.

And the C++ implementation would be :

void sieve(int N) { int x=sqrt(N),i,j; vector<bool> S(N+1,0); for (i=4;i<=N;i+=2) S[i]=1; for (i=3;i<=x;i+=2) if (!S[i]) for (j=i*i;j<=N;j+=2*i) S[j]=1; long long s=0; for (i=2;i<=N;i++) if (!S[i]) s+=i; printf("%lld\n",s); }

But then we can optimize even the Sieve of Eratosthenes. If you look closer you will realise that apart from 2 all the other even numbers are of no use and just waste our time and memory. So we should get rid of them. We can do that by sieving only odd numbers. Let an element with index i correspond to number 2*i+1. So in the sieve we only need to go upto (N-1)/2. Also if p = 2*i+1 *p*^{2} = 4*i^2 + 4*i + 1 which is represented by index 2*i*(i+1). Also the step will be of order 2*i+1. So the code would look something like this:

void improved_sieve(int N) { int M=(N-1)/2; int x=(floor(sqrt(N))-1)/2,i,j; vector<bool> S(M+1,0); for (i=1;i<=x;i++) if (!S[i]) for (j=2*i*(i+1);j<=M;j+=(2*i+1)) S[j]=1; long long s=2; for (i=1;i<=M;i++) if (!S[i]) s+=(2*i+1); printf("%lld\n",s); }

and since i am learning Python. Here’s the python code as well :

import math def improved_sieve(N): M=(N-1)/2 x=(int(math.sqrt(N))-1)/2 S=[] for i in range(M+1): S.append(False); for i in range (1,x+1): if (S[i]==False): for j in range (2*i*(i+1),M+1,2*i+1): S[j]=True s=2; for i in range (1,M+1): if (S[i]==False): s+=(2*i+1) print s N=input() improved_sieve(N)

Lets compare all these codes on runtime :

Limit10^N where N= | basic1 | basic2 | moderate1 | moderate2 | sieve | improved_sieve |

4 | 0.079 | 0.026 | 0.003 | 0.007 | 0.002 | 0.001 |

5 | 7.562 | 2.611 | 0.050 | 0.117 | 0.020 | 0.011 |

6 | – | – | 1.160 | 2.218 | 0.208 | 0.121 |

7 | – | – | 26.350 | 44.946 | 2.183 | 1.269 |

All time are in seconds. If you notice that though algorithm 4 seems to be better than 3 it takes more time because of the array refrences that are made several times. Choose the best.

Happy Coding!

-fR0D

hey, in seive of eratosthenes function, suppose N is of order of 10^8 or 10^9, then vector S(N+1,0); this line will cause error, how to deal with such large numbers?

taruniiitaOctober 28, 2009 at 9:07 PM

It doesn’t give any error on my machine. If it gives error on your’s try declaring the S array globally as bool S[1000000000]; . May be that will help.

fR0DOctober 29, 2009 at 11:40 AM

maximum array size is 5*(10^5). so how can you declare array of size 10^9 ?

taruniiitaOctober 29, 2009 at 1:34 PM

Tarun i don’t know how did you get that number, but to get the maximum size of the vector array’s you can do this :

fR0DOctober 29, 2009 at 8:09 PM

hey, just run this, this gives segmentation fault. when we take array size > 5*(10^5) this gives segmentation fault.

taruniiitaOctober 29, 2009 at 8:23 PM

You cannot declare that big a static array locally. If you remember i wrote in my earlier post, declare the array globally. Have you tried this :

fR0DOctober 29, 2009 at 8:27 PM

okay, got that. thanks 🙂

taruniiitaOctober 29, 2009 at 8:30 PM