24
Are you crazy about primes?
Algorithm Design
Primes

Sieve of Eratosthenes is a well known prime generation algorithm: link : http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

Using the property that all prime numbers essentially fall in the "6k+1" or "6k-1" series, we can reduce the space and time complexity a little bit.

We need to distinguish the composite numbers that occur in "6k+1" and "6k-1" series using an Algorithm.

Concept

To list all prime numbers less than "10" we use the following representation

  • A bitset of size 10 = 1111111111 (initially marking all numbers as prime)
  • After processing the bitset = 0001010100 (later we set those bits to zero whose index are composite numbers and index starting from "1")
  • To illustrate the method we use an example, which is to generate all prime numbers less than "31700"

Approach 1:

  • We take a bitset of size "31700" and initialize all bits to "1" .
  • Mark "0" and "1" as composite numbers.
  • Mark all multiples of "2" as composite numbers i.e. the index values 4,6,8 ... are set to "0".
  • Find the next index after "2" that happens to be set to "1". It is the index "3". We can mark out its multiples 6, 9, 12, 15... to "0".
  • Similarly the next index would be "5" and we can mark out its multiples.
  • Repeat this for until the last index less than square_root(31700) and is
    marked as prime (set to "1").

Approach 2:

We change our representation here by mapping each index of the bitset with "k" in "6k+1" or "6k-1" series, which is as follows:

We take two bitsets, one mapping to "6k+1" series and one more mapping to "6k-1" series. If I called these bitsets "u" and "v" respectively:

  • In a bitset u = 1111 of size "4" index starting from "0":
  • u[1] maps to 6*(1)+1 = 7
  • u[2] maps to 6*(2)+1 = 13
  • u[3] maps to 6*(3)+1 = 19 and so... on

and

  • In a bitset l = 1111 of size "4" index starting from "0":
  • l[1] maps to 6*(1)-1 = 5
  • l[2] maps to 6*(2)-1 = 11
  • l[3] maps to 6*(3)-1 = 17 and so.. on

Computing size of the bitset needed !!

This is the maximum value that 'k' can take when I need to map "6k-1" to the number close to 31700 . we compute this as follows:

  • 6*k - 1 = 31700 => k = (31700+1)/6 => 5283
  • 6k + 1 = 6(5283) + 1 =31699 is the last probable prime number less than 31700
  • The size of the bitset would be '5184' since indexing starts from '0' [A size of 10 can hold '0' to '9' when '0' is included]

Case 1: In "6*k+1" series we observe the follows:

series: *7 13 19 25 31 37 43 49 55 61 67 73 79 85 91 97 103 109* ...

Multiples of '7' are distanced '7' units apart from itself, Multiples of '13' are distanced '13' units apart from itself .(units being the numbers in the series) This way we can set multiples of '7' , '13' etc to '0' in the "6k+1" series.

2:We need to mark-out the multiples of '7', '13' etc in the "6k-1" series as well

series: 5 11 17 23 29 35 41 47 53 59 65 71 77 83 89 95 101 107 113 119 ...

we observe the first multiple of '7' occurs at a position 'k1' where k1=7-k, [k-position where '7' occurs i.e in "6k+1" series. ]

  • 6*(1)+1 = 7 => k=1;
  • k1 = 7- k => 7-1 = 6; [ k1 is the position of the first multiple of '7' ( '7' belongs to "6k+ 1" series), in the "6k-1" series].
  • 6(k1) - 1 = 6(6) - 1 =35 ;
  • '35' is the first multiple of '7' in the "6k-1" series.

similarly,

  • 6*(2)+1 = 13 => k=2;
  • k1= 13- k =>13-2 = 11; [ k1 is the position of the first multiple of '13' ( '13' belongs to "6k+ 1" series), in the "6k-1" series].
  • 6(k1) -1 = 6(11)-1 = 65;
  • '65' is the first multiple of '13' in "6k-1" series.

Once we find 'k1' we can successively mark out the intervals at distance of '7' or '13' respectively from itself in the "6k-1" series to remove their multiples.

Case 2:

In "6k-1" series we need to perform similar operations:

series: *5 11 17 23 29 35 41 47 53 59 65 71 77 83 89 95 101 107 113 119 ...*

Mark-out multiples of '5', '11' etc in "6k-1" series which follow the same principle as above.Multiples of '5' lie at a distance of '5' units from itself as shown in the above series.

We need to mark-out the multiple of '5' in "6k+1" series also.

series: 7 13 19 25 31 37 43 49 55 61 67 73 79 85 91 97 103 109 ...

Mark-out multiples of '5' , '11' etc in "6k-1" series similarly.The first multiple of 5 in "6k+1" series lies at a position k1 where k1=5-k .

  • 6*(1)- 1 = 5 => k=1;
  • k1 = 5- k => 5-1 = 4; [ k1 is the position of the first multiple of '5' ( '5' belongs to "6k- 1" series), in the "6k+1" series].
  • 6(k1) - 1 = 6(4) + 1 =25 ;
  • '25' is the first multiple of '5' in the "6k+1" series.

[ you can reason this logic easily with a little bit of work around the formula 6k+1 and 6k-1 :)]

Last Step :

We know how to mark-out all the composites , but how much do we need to iterate?

For n=31700 we iterate up to square_root(31700) = 178 Here we need to iterate with k i.e. 6k-1 = 178 => k = (178+1)/6 = 29

Within approximately 30 basic iterations we can generate all prime numbers less than 31700, or infact with fewer basic iterations generate prime numbers less than a fairly large value of 'n'.

Code

#include<iostream>
#include<bitset>
#include<fstream>
#include<ctime>
using namespace std;
int main()
{
bitset<5284> u,l;
u.set();
l.set();
int start= clock();
for(int k=1;k<=29;k++)
{

      if(u[k]==1)
      {
            int m=6*k+1;int n=m-k;
            l.set(n,0);
            int i=1;
            while(k+ i*m <= 5283 && n+ i*m <= 5283 )
            {
                  u.set(k+i*m, 0);
                  l.set(n+i*m , 0);
                  i++;
            }
            if(k+i*m <= 5283)
                  u.set(k+i*m,0);
            if(n+i*m <= 5283)
                  l.set(n+i*m,0);
      }
      if(l[k]==1)
      {
            int m=6*k-1; int n=m-k;
            int i=1;
            u.set(n,0);
            while(k+i*m<= 5283 && n+i*m <= 5283)
            {
                  l.set(k+i*m,0);
                  u.set(n+i*m,0);
                  i++;
            }
            if(k+i*m <= 5283)
                  l.set(k+i*m,0);
            if(n+i*m <= 5283)
                  u.set(n+i*m,0);
      }
}
int stop=clock();
cout<<"time "<<(stop-start)/double(CLOCKS_PER_SEC);
fstream f("p1");
f<<"2\n3\n";
for(int i=1;i<=5283;i++)
{
if(l[i])
      f<<6*i-1<<endl;
if(u[i])
      f<<6*i+1<<endl;
}
}
//open the file named "p1" to view the list of prime numbers
//hard-coding is used for illustration.
Author

?