Algorithm Implementation/Mathematics/Prime number generation

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

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  , 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   with linked lists.

Recast as unbounded definition, yet faster and with better complexity due to the tree-shaped folding:

primesU = 2 : _Y ((3 :) . minus [5,7..] . _U . map (\p-> [p*p, p*p+2*p..]))

_Y g = g (_Y g)                                      -- g . g . g . g . .... 
_U ((x:xs):t) = x : (union xs . _U . pairs) t        -- big_U, or unionAll 
   where  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 (also found in data-ordlist package's Data.List.Ordered module) 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, tails)
import Data.Array
 
primesSA = _Y (\ps -> 
                2 : [n | (r:q:_, px) <- (zip . tails . (2:) . map (^2)) ps (inits ps),
                         (n,True)    <- assocs (
                                          accumArray (\_ _ -> False) True (r+1,q-1)
                                    [(m,()) | p <- px, 
                                              let s=(r+p)`div`p*p, m <- [s,s+p..q-1]] )])

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

Phix edit

constant limit = 1000
sequence primes = {}
sequence flags = repeat(1, limit)
for i=2 to floor(sqrt(limit)) do
    if flags[i] then
        for k=i*i to limit by i do
            flags[k] = 0
        end for
    end if
end for
for i=2 to limit do
    if flags[i] then
        primes &= i
    end if
end for
? primes

You can get exactly the same output (the first 168 primes) from

?get_primes(-168)

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