COME ON CODE ON

A blog about programming and more programming.

Prime Numbers

with 7 comments

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 :

  1. Create a contiguous list of numbers from two to some highest number n.
  2. Strike out from the list all multiples of two (4, 6, 8 etc.).
  3. The list’s next number that has not been struck out is a prime number.
  4. Strike out from the list all multiples of the number you identified in the previous step.
  5. 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).
  6. 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 p2 = 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

Advertisements

Written by fR0DDY

March 10, 2009 at 3:34 PM

Posted in Maths, Programming

Tagged with , , , ,

7 Responses

Subscribe to comments with RSS.

  1. 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?

    taruniiita

    October 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.

      fR0D

      October 29, 2009 at 11:40 AM

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

        taruniiita

        October 29, 2009 at 1:34 PM

  2. 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 :

    #include<iostream>
    using namespace std;
    #include<vector>
    
    int main()
    {
        vector<bool> B;
        cout<<B.max_size()<<endl;
        vector<int> I;
        cout<<I.max_size()<<endl;
    }
    

    fR0D

    October 29, 2009 at 8:09 PM

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

    #include<iostream>
    using namespace std;
    
    int main()
    {
    	int a[10000000];
    	for (int i = 0; i < 10000000; i++) {
    		a[i] = i+1;
    	}
    	cout << a[0] << endl;
    }
    

    taruniiita

    October 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 :

      #include<iostream>
      using namespace std;
      
      int a[10000000];
      int main()
      {
      for (int i = 0; i < 10000000; i++) {
      a[i] = i+1;
      }
      cout << a[0] << endl;
      }

      fR0D

      October 29, 2009 at 8:27 PM

      • okay, got that. thanks 🙂

        taruniiita

        October 29, 2009 at 8:30 PM


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: