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]);
}
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 , 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
editsimple 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
editconstant 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
editdef 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
editdef 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
editlet 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)
editFast 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);
}