# Algorithm Implementation/Mathematics/Prime number generation

## C

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

```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. 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 `union`s. 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

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):
"""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

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