Algorithm Implementation/Mathematics/Prime number generation

From Wikibooks, open books for an open world
Jump to: navigation, search

C[edit]

#define MAX 10000
 
int i, j, sieve[MAX], primecount=0, prime[MAX];
 
for(i=0; i<MAX; i++){
	sieve[i]=1;
}
sieve[0]=sieve[1]=0;
for(i=2; i<MAX; i++){
	while(sieve[i]==0 && i<MAX){i++;}
	prime[primecount]=i;
	for(j=i*i; j<MAX; j+=i){sieve[j]=0;}
	primecount++;
}

D[edit]

int[] sieve(int n) {
    int[] sieve = new int[](n + 1); 
    int m = cast(int) sqrt(n);
 
    sieve[1..3] = 1;
 
    for (int i = 3; i <= m; i+=2) {
        for ( size_t j = 2; j < cast(int)sqrt(i); ++j) {
            if (sieve[i] == 1 || i % j != 0) {
                 continue;
            }
        }
        sieve[i] = 1;
    }
 
    return sieve;
}

Haskell[edit]

For each number, check whether it is divisible by any of the primes not above its square root. This is an optimal trial division algorithm:

isPrime primes n = foldr (\p r -> p*p > n || (rem n p /= 0 && r))
                         True primes
primes = 2 : filter (isPrime primes) [3..]

Note that this code utilizes Haskell's lazy evaluation, to allow primes and isPrime to be defined in terms of each other.

The Sieve of Erastosthenes has better theoretical complexity than trial division. Its intent demonstrated by:

primesTo m = 2 : sieve [3,5..m] where
  sieve s@(p:xs) | p*p>m = s
                 | True  = p : sieve (xs `minus` [p*p, p*p+2*p..m])

Were a (minus a b) operation to have a complexity of O(|b|), as it indeed is so on mutable arrays in imperative languages, this code would achieve the theoretical complexity of the S. of E. Unfortunately it is O(|a|) here. So this code ends up having same complexity as optimal trial division, although it still runs faster.

Recast as unbounded definition, yet faster and with much better complexity:

primesU () = 2 : [3,5..] `minus` unionAll [[p*p, p*p+2*p..] | p <- primes'] 
  where
    primes' = 3 : [5,7..] `minus` unionAll [[p*p, p*p+2*p..] | p <- primes']
    unionAll ((x:xs):t) = x : union xs (unionAll (pairs t))
    pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t

This replaces the linear nested chain of minus operations in primesTo function, with an unbounded right-deepening tree of unions. A sizable constant factor speed-up can be further achieved with wheel factorization. The utility functions used above are:

minus (x:xs) (y:ys) = case compare x y of LT -> x : minus  xs (y:ys)
                                          EQ ->     minus  xs    ys
                                          GT ->     minus (x:xs) ys
minus  a      b     = a
union (x:xs) (y:ys) = case compare x y of LT -> x : union  xs (y:ys)
                                          EQ -> x : union  xs    ys
                                          GT -> y : union (x:xs) ys
union  a      b     = a ++ b

Java[edit]

simple naive algorithm

import java.util.*;
 
// finds all prime numbers up to max
public static List<Integer> generatePrimes(int max)
{
    List<Integer> primes = new ArrayList<Integer>();
 
    // start from 2
    OUTERLOOP:
    for (int i = 2; i <= max; i++) {
        // try to divide it by all known primes
        for (Integer p : primes)
            if (i % p == 0)
                continue OUTERLOOP; // i is not prime
 
        // i is prime
        primes.add(i);
    }
    return primes;
}

Sieve of Erastosthenes, from here

import java.util.*;
 
public class Sieve
{
    private BitSet sieve;
 
    private Sieve() {}
 
    private Sieve(int size) {
        sieve = new BitSet((size+1)/2);
    }
 
    private boolean is_composite(int k)
    {
        assert k >= 3 && (k % 2) == 1;
        return sieve.get((k-3)/2);
    }
 
    private void set_composite(int k)
    {
        assert k >= 3 && (k % 2) == 1;
        sieve.set((k-3)/2);
    }
 
    public static List<Integer> sieve_of_eratosthenes(int max)
    {
        Sieve sieve = new Sieve(max + 1); // +1 to include max itself
        for (int i = 3; i*i <= max; i += 2) {
            if (sieve.is_composite(i))
                continue;
 
            // We increment by 2*i to skip even multiples of i
            for (int multiple_i = i*i; multiple_i <= max; multiple_i += 2*i)
                sieve.set_composite(multiple_i);
 
        }
        List<Integer> primes = new ArrayList<Integer>();
        primes.add(2);
        for (int i = 3; i <= max; i += 2)
            if (!sieve.is_composite(i))
                primes.add(i);
 
        return primes;
    }
}

Javascript[edit]

<!DOCTYPE html>
<html>
<body>
 
<button onclick="myFunction()">Sieve</button>
 
<script>
var sieve = [false, false];
function myFunction(){
  var n = sieve.length, nn = n * n;
  for(var i = n; i < nn; i++){ sieve.push(true); }
  for(var p = 0; p < n; p++){
    if(sieve[p]){
      for(var m = p * p; m < nn; m += p){
        sieve[m] = false;
      }
    }
  }
  var txt = document.body.
   appendChild(document.createElement('p')).
   appendChild(document.createTextNode(''));
  for(var i = n; i < nn; i++){
    if(sieve[i]){ txt.nodeValue += ' ' + i; }
  }
}
</script>
 
</body>
</html>

Perl[edit]

$MAX = 10000;
 
for($i=2; $i<$MAX; $i++){
	$sieve[$i]=1;
}
$sieve[0]=0;
$sieve[1]=0;
$primecount=0;
for($i=2; $i<$MAX; $i++){
	while($sieve[$i]==0 && $i<$MAX){$i++;}
	$prime[$primecount]=$i;
	for($j=2*$i; $j<$MAX; $j+=$i){$sieve[$j]=0;}
	$primecount++;
}

Python[edit]

def prime_sieve(n):
    """Python implementation of the Sieve of Eratosthenes"""
    n = int(n)
    primelist=[]
    markedlist=[False]*(n+1)
    for i in range(2,n+1):
        if not markedlist[i]:
            primelist.append(i)
            for m in range(i**2, n+1, i):
                markedlist[m] = True
    return primelist

Faster variation:

def sieve(n):
    s = xrange(3, n + 1, 2)
    r = set(s)
    [r.difference_update(xrange(n << 1, s[-1] + 1, n)) for n in s if n in r]
    return r.union([2])

Relatively slow but elegant .. .. so use it for educational and scientific research purposes .. it includes a time measurement of the calculation

def pgenl(end=100):
  '''
    generates the sequence of all primes (lower end): 2,3,5 ... < end
    ... based on Wilson' theorem
  '''
  l=[]
  p=2 
  f=1 
 
  try:
 
    ts2=time.time()
 
    while(p<end or end==True):
 
      if f%p%2: l.append(p)
      p+=1
      f*=(p-2)
 
    ts2-=time.time()
    sys.stderr.write("time: "+str(round(-ts2,3))+" sec\n")
    return l
 
  except:
 
  #   raise KeyboardInterrupt .. Ctrl-C
 
    ts2-=time.time()
    sys.stderr.write("time: "+str(round(-ts2,3))+" sec\n")
    return l

Ruby[edit]

def sieve_of_eratosthenes(max)
  arr=(2..max).to_a
  (2..Math::sqrt(max)).each do |i|
     arr.delete_if {|a|a % i == 0 && a!=i}
   end
   arr
end
 
# Example Call
puts sieve_of_eratosthenes(20)

Note that Ruby has a built-in Sieve of Erastosthenes prime number generator:

require 'mathn'
for p in Prime.new
  puts p
end

F# / OCaml[edit]

let prime_series_sieve (limit:bigint) = 
    let series = List.to_array [0I..limit]
    series.SetValue(0I, 1)
 
    let rec eliminate_multiples (n:bigint) (i:bigint) = 
        let index = (i * n)
        if index < bigint.Parse(series.Length.ToString()) then 
            series.SetValue(0I, (bigint.ToInt64(index)))
            eliminate_multiples n (i + 1I)
 
    for n in [2I..(limit/2I)] do
        eliminate_multiples n 2I
 
    series;;

C (bitwise)[edit]

#include <math.h>
#define MAX 1000000
const int S=(int)sqrt((double)MAX);
 
#define rP(n) (sieve[n>>6]|=(1<<((n>>1)&31)))
#define gP(n) sieve[n>>6]&(1<<((n>>1)&31))
 
unsigned sieve[(MAX>>6)+1]={0};
int i, j,k,l=0 , prime[MAX];
prime[0]=2;
for(i=3;i<=S;i+=2) if(!(gP(i))) {
	k=(i<<1);
	prime[l++]=i;
	for(j=i*i;j<=MAX;j+=k) rP(j);
}
 
for(;i<=MAX;i++) if(!(gP(i))) prime[l++]=i;

C (Fast bitwise)[edit]

Fast implementation of bitwise Sieve of Eratosthenes. Little bit cheating, because it's counts only odd primes (read as: you need to add first prime - "2" manually).

There is also OpenMP #pragma's along the way, so you can try compile it with -fopenmp

#include <stdio.h>
#include <stdlib.h>
#include <string.h> 
 
#define MAX_N 1000000000
#define MAX MAX_N/2 
 
/*
Good way to compile is:
gcc --std=c99 -O4 -o sieve_fast sieve_fast.c  
*/
 
int main(int argc, char * args[])
{
  size_t mem_sz = MAX / 8 + 1;
  size_t superpage_size = 2 * 1024 * 1024; // 2MB
  unsigned long cnt = 0;
  unsigned char * arr;
 
  if (posix_memalign((void **)&arr, superpage_size, mem_sz)) {
    perror("memalign");
    exit(1);
  }
 
  bzero(arr, mem_sz);
 
  #pragma omp parallel 
  for(register unsigned long n=1; ((n*(n+1))<<1)<MAX; n++)
  {
    register unsigned long t = n >> 3;
    if((arr[t] & (1 << (n & 7)))){
      //printf("Skipped: %lu\n", n*2+1);
      continue;
    }
    //printf("prime: %lu; n: %lu\n", n<<1+1, n);
    #pragma omp for schedule(static, 65536)
    for(register unsigned long c=(n*(n+1))<<1; c<MAX; c+=(n<<1)+1){
      t = c >> 3;
      arr[t] |= (1 << (c & 7));
      //printf("Sieving %lu with step %lu\n", c*2+1, n*2+1);
    }
  }
 
  #pragma omp parallel for reduction(+:cnt)
  for(register unsigned long n=1; n<MAX; n++)
  {
    register unsigned long t = n >> 3;
    if(!(arr[t] & (1 << (n & 7))))
    {
      //printf("%lu, ", n*2+1);
      cnt++;
    }  
  }
 
  printf("Found total: %lu\n", cnt);
}

See also[edit]