Algorithm Implementation/Mathematics/Prime number generation

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

C[edit]

#include <stdbool.h>
#include <stdio.h>

#define MAX 10000

/*
 * Prints all primes less than MAX using the Sieve of Eratosthenes.
 */
int
main(void)
{
	bool sieve[MAX];
	int i, j, primecount = 0, prime[MAX];

	for (i = 0; i < MAX; i++)
		sieve[i] = true;
	sieve[0] = sieve[1] = false;
	for (i = 2; i < MAX; i++) {
		if (!sieve[i])
			continue;
		prime[primecount++] = i;
		for (j = i * i; j < MAX; j += i)
			sieve[j] = false;
	}
	for (i = 0; i < primecount; i++)
		printf("%d\n", prime[i]);
}

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 \cup b|) here.

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

primesU = 2 : _Y ((3 :) . minus [5,7..] . _U . map (\p-> [p*p, p*p+2*p..]))
  where
    _Y g = g (_Y g)                                      -- g . g . g . g . .... 
    _U ((x:xs):t) = x : (union xs . _U . pairs) t        -- big_U, or unionAll 
    pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t

This replaces the linear left-leaning nested chain of minus operations in primesTo function, with an unbounded right-deepening tree of progressively growing balanced trees of unions. A sizable constant-factor speedup 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      _     = 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      []    = a
union  []     b     = b

The removal of composites is easy with arrays. Processing is naturally broken into segments between squares of consecutive primes:

import Data.List (inits)
import Data.Array.Unboxed
 
primesSAE = 2 : sieve 3 4 (tail primesSAE) (inits primesSAE)
  where
  sieve x q ps (fs:ft) = [i | (i,True) <- assocs (
         accumArray (\ _ _ -> False)
                    True  (x,q-1)
                    [(i,()) | p <- fs, let c = p * div (x+p-1) p,
                              i <- [c, c+p..q-1]] :: UArray Int Bool )]
      ++ sieve q (head ps^2) (tail ps) ft

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):
    """Generate the primes less than ``n`` using the Sieve of Eratosthenes."""
    a = [True] * n
    a[0] = a[1] = False
    for i, isprime in enumerate(a):
        if isprime:
            yield i
            for j in range(i * i, n, i):
                a[j] = False


def primes_wilson(n):
    """Generate the primes less than ``n`` using Wilson's theorem."""
    fac = 1
    for i in range(2, n):
        fac *= i - 1
        if (fac+1) % i == 0:
            yield i

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);
}