Algorithm Implementation/Mathematics/Prime number generation

C

```#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

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

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 ${\displaystyle 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 ${\displaystyle O(|a\cup b|)}$ 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 `union`s. 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

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
}
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>();
for (int i = 3; i <= max; i += 2)
if (!sieve.is_composite(i))

return primes;
}
}
```

Javascript

```<!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

```\$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

```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

```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

```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)

```#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)

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