Sieve of Eratosthenes
February 19, 2009
Over two millenia ago, Eratosthenes, who calculated the circumference of the earth, the distance to the Sun and the tilt of the Earth’s axis, developed a system of latitude and longitude, and invented the leap day, created a systematic method to enumerate the prime numbers that is still in use today. Eratosthenes was born in Cyrene (present-day Libya), lived from 276 B.C. to 194 B.C., and spent most of his life in Alexandria, Egypt, where he was the second Chief Librarian of the Great Library, succeeding Apollonius of Rhodes; he was a good friend of Archimedes.
The Sieve of Eratosthenes starts by making a list of all the numbers up to a desired maximum; we’ll illustrate the method by calculating the prime numbers through thirty:
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
Now take the first number on the list, 2, and cross off every second number:
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
(Although it may not be obvious, the number 4 is crossed off the list; in some fonts, the cross-bar of the 4 coincides with the strike-through bar.) Next, take the next number on the list that isn’t crossed off, 3, and cross off every third number; some of them have already been crossed off:
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
Repeat that last step for the next un-crossed number on the list, 5:
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
And so on, each time crossing off all multiples of the next un-crossed number on the list. The list of prime numbers are all those that haven’t been crossed off:
2 3 5 7 11 13 17 19 23 29
This method is called a sieve because it sweeps through a range of numbers, with each prime number, as it is discovered, blocking all its multiples from falling through as prime numbers. The sieve admits several optimizations. First, only odd numbers are considered, since the initial sifting crosses off all the even numbers except 2, which is handled separately. Second, crossing off starts at the square of the number being sifted, since all smaller primes have already been crossed off by previous steps of the sieve; for instance, sifting by 3 starts at 9, since 6 was already crossed off when sifting by 2. Third, sifting stops at the square root of the maximum number in the sieve, since any non-primes larger than the square root must have already been crossed off at previous levels of the sieve; thus, in the above example there is no need to sieve on the prime number 7, or any larger prime number, since the square of 7 is greater than 30, which is the largest number in the list.
Write a function that takes a single argument n and returns a list of prime numbers less than or equal to n using the optimized sieving algorithm described above. Apply the function to the argument 15485863 and count the number of primes returned.
Programming Praxis 
My haskell solution : http://codepad.org/OMywU3lJ
In C. http://programmingpraxis.pastebin.com/f6f889fa4
[…] this is the second Programming Praxis exercise I’ve done, and this one was slightly more difficult to get right, or at least test. […]
I’m confused by one aspect of the above solution. I don’t understand why startj is defined as (+ (* 2 i i) (* 6 i) 3).
I understand why p is (+ i i 3) – that’s the sequence of odd numbers, starting at 3. But what relationship does that have with startj?
Can anyone enlighten me?
That’s the implementation of the second optimization.
Suppose you have already sieved 2, 3, and 5 and are now beginning to sieve 7. Start by adding 7 to the list of primes. Then 14, but that has already been sifted out by 2; likewise, 21 has been sifted out by 3, 28 has been sifted out by 2, 35 has been sifted out by 5, and 42 has been sifted out by 2 (and also 3, but that doesn’t matter). The first multiple of 7 that gets sifted out is 49, which is 7 times 7. In general, when sifting by n, sifting starts at n-squared because all the previous multiples of n have already been sieved.
The rest of the expression has to do with the cross-reference between numbers and sieve indexes. There’s a 2 in the expression because we eliminated all the even numbers before we ever started. There’s a 3 in the expression because Scheme vectors are zero-based, and the numbers 0, 1 and 2 aren’t part of the sieve. I think the 6 is actually a combination of the 2 and the 3, but it’s been a while since I looked at the code, so I’ll leave it to you to figure out.
Good question!
The number 2 is handled specially, so the vector v looks this:
index 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23number 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49
When you want to begin sieving 7, the first number you sift out is 7 × 7 = 49. The number 7 is at index i = 2 in the vector, and the number 49 is at index startj = (+ (* 2 i i) (* 6 i) 3) = 23 in the vector.
So, it seems we can alternatively define startj as:
(quotient (- (expt p 2) 3) 2)That is,
– p-squared, to give us the number to start sifting at
– shifted by 3 because we ignore the first 3 numbers (0, 1 and 2)
– divided by 2 to get the index, because we ignore every second number.
Is that correct?
The straightforward implementation in C, storing all the odd numbers in an array and then zeroing all non-primes. Not being much of a programmer, I’ll admit I’m amazed my laptop can compute 1,000,000 primes in under half a second… And I suppose the code could even be made to run twice as fast on my c2duo by having two threads sifting through the array.
#include <stdio.h> #include <stdlib.h> #include <math.h> #define index(A) (A - 3) / 2 //index of A in the array of odd numbers #define vector(A) 2 * A + 3 //inverse operation of index() long *fill_num(long n); //Generates array of odd numbers > 2 int seek_primes(long *numbers, long n); //Zeroes all non-primes int main() { long ulimit; long *nombres; int premiers; char c; printf("Upper limit for the calculation of primes?\n"); scanf("%ld", &ulimit); ulimit = index(ulimit) + 1; //+ 1: the array is indexed from 0 nombres = fill_num(ulimit); premiers = seek_primes(nombres, ulimit) + 1; //2 is also prime printf("Found %d primes\n", premiers); free(nombres); return 0; } int seek_primes(long *numbers, long n) { long i, j, primes = 0, max; //Third optimisation: sift until sqrt of upper limit printf("sift until %ld\n", max = (long)sqrt(vector(n))); for (i = 0; i < n; i++) if (numbers[i] != 0) { primes++; if (i < max) { j = numbers[i] * numbers[i]; while (j < vector(n)) { numbers[index(j)] = 0; j = j + 2 * numbers[i]; } } } return primes; } long *fill_num(long n) //n is the size of the array, vector(n - 1) the highest odd number { long i; long *tmp; tmp = (long *)malloc(n * sizeof(long)); for (i = 0; i < n; i++) tmp[i] = vector(i); return tmp; }…turns out the multi-threaded version does not work, at least the way I’ve implemented it. I’ve effectively split the seek_primes in two threads, one sifting for (i = 0; i < n; i+=2), and the other one for (i = 1; i < n; i+=2).
It appears my first thread finds non-zero elements in the numbers[] array in seek_primes, but there is no guarantee these same elements are not going to be zeroed by the second thread later.
In fact up to 100 I find 3 primes too many, and up to 1000 about 38 (and of course the actual number found might even be different from one run to the next).
My python solution is at http://pastebin.com/Qup0TuFu.
I generated all primes from the integer list of 0-100000. Since the last prime in that list is 99991, using the square of the primes method, it is relatively easy and fast deriving the prime sequence of most numbers. I think I began having problems around the number of primes for the integer 9990000000.
Here’s my Clojure version:
(require ‘[clojure.contrib.math :as math])
;; faster for dropping a few leading nils than (into [] (drop-while nil? sieve))
(defn drop-while-nil? [sieve] (if (first sieve) sieve (recur (subvec sieve 1))))
;; returns sieve with all multiples of (first sieve) set to nil,
;; leading nils removed
(defn remove-every-nth [sieve]
(let [count-sieve (count sieve), n (first sieve)]
(loop [s sieve, i 1, ni n]
(if (> ni count-sieve)
(drop-while-nil? (subvec s 1))
(recur (assoc s ni nil), (inc i), (* n (inc i)))))))
(defn eratosthenes [r]
(loop [s (into [] (range 2 r)), primes [1] ]
(if (> (peek primes) (math/sqrt r))
(concat primes (remove nil? s))
(recur (remove-every-nth s), (conj primes (first s))))))
(map #(count (time (eratosthenes %))) [10 100 1000 10000 100000 ])
Simple and simple tail-recursive solutions, but missing the “start eliminating at (* p p)” optimization:
#lang racket ;; .. : num num -> list-of-nums ;; given first and last, return a list of the range first .. last, inclusive (define (.. first last) (if (= first last) (list first) (cons first (.. (+ first 1) last)))) (define (erat l n) (let ([p (car l)]) (if (> (sqr p) n) l (cons p (erat (filter (lambda (n) (not (= (modulo n p) 0))) (cdr l)) n))))) (define (tail-erat l n accum) (let ([p (car l)]) (if (> (sqr p) n) (append (reverse accum) l) (tail-erat (filter (lambda (n) (not (= (modulo n p) 0))) (cdr l)) n (cons p accum))))) ; main program (let* ([args (vector->list (current-command-line-arguments))] [n (string->number (car args))] [l (.. 2 n)]) ; (display (erat l n)) (display (tail-erat l n '())) (newline))Did you calculate the number of primes less than 15485863? I think it will take a while.
Yes, it takes about fifty times longer than the vector-based solution above, mostly in GC.
Using Petite Chez Scheme for profiling:
(time (display (length (…))))
6 collections
3021 ms elapsed cpu time, including 248 ms collecting
3022 ms elapsed real time, including 247 ms collecting
110325408 bytes allocated, including 40188320 bytes reclaimed
(time (length (tail-erat (…) …)))
2767 collections
167705 ms elapsed cpu time, including 103201 ms collecting
167817 ms elapsed real time, including 103268 ms collecting
23602922896 bytes allocated, including 24373515936 bytes reclaimed
Moving to an in-place filtering of the list (see below) brought this time down by a factor of ten, though it’s still slower than the vector version:
(time (length (tail-erat (…) …)))
105 collections
11498 ms elapsed cpu time, including 5577 ms collecting
11503 ms elapsed real time, including 5590 ms collecting
1135547536 bytes allocated, including 3928896 bytes reclaimed
This was achieved by replacing filter with the following in-place filter! routine (and staying in Chez Scheme, where set-cdr! is available).
A more general filter!, implemented as syntax! so it could consistently alter the list in place but not depend on returning a list would be interesting, but isn’t needed for this — this implementation returns the result of modifying the list in place, but the original value passed in is only modified from the first value not matching pred? on…
(define (filter! pred? l) (cond [(null? l) '()] [(pred? (car l)) (filter! pred? (cdr l))] [else (let loop ([follower l] [leader (cdr l)]) (cond [(null? leader) l] [(pred? (car leader)) (set-cdr! follower (cdr leader)) (loop follower (cdr leader))] [else (loop leader (cdr leader))]))]))You’re missing the point. The Sieve is fast because it is based on addition. You are using division (the modulo function). Your algorithm is not the Sieve of Eratosthenes.
Oh, I see… let me revisit…
Naive attempt at moving to vector and using addition over vector indices gives a slight speedup, looking further:
(time (length (primes 15485863)))
4 collections
10296 ms elapsed cpu time, including 366 ms collecting
10300 ms elapsed real time, including 366 ms collecting
155888768 bytes allocated, including 124098560 bytes reclaimed
(define (primes n) (define (zeroing-pass! vec start step) (cond [(> start n) vec] [(not (vector-ref vec (sub1 step))) vec] [else (vector-set! vec (sub1 start) #f) (zeroing-pass! vec (+ start step) step)])) (define (erat! vec) (let* ([stop (div n 2)]) (let loop ([p 2]) (if (= p stop) vec (begin (zeroing-pass! vec (+ p p) p) (loop (+ p 1))))))) (let ([vec (make-vector n #t)]) (erat! vec) (reverse (let loop ([i 1] [l '()]) (if (= i n) l (loop (add1 i) (if (vector-ref vec i) (cons (add1 i) l) l)))))))I got curious as to how big a difference the additon-vs-module distinction is on modern hardware (historically, it has been very large, as you point out).
I used the following code to run 10^8 integer additions, modulo operations, and modulo operations with a power-of-two modulus (another traditional distinction in performance).
I may well be missing something big, but the output shows about a 25% difference in speed between additions and modulo operations. This obviously adds up, but looking at the above code, it seems to be small compared to other differences (particularly, the various operations which control how many of whichever operation are carried out, and which reduce the infrastructural cost of the algorithm).
(In the code below, I use three vectors in an attempt to minimize the effect of cache read-ahead on subsequent vector-maps during earlier testing with a smaller n.)
(let* ([n (expt 10 6)] [times 100] [time-reps (lambda (f v) (time (do ([i 0 (+ i 1)]) ([= i times]) (vector-map f v))))] [v1 (make-vector n 1203006)] [v2 (make-vector n 1203006)] [v3 (make-vector n 1203006)]) (display "Timing addition\n") (time-reps (lambda (n) (+ n 57365)) v1) (newline) (display "Timing modulo (power of 2)\n") (time-reps (lambda (n) (mod n 65536)) v2) (newline) (display "Timing modulo (non-power of 2)\n") (time-reps (lambda (n) (mod n 57365)) v3) (newline))Timing addition (time (do ((...)) ...)) 299 collections 6471 ms elapsed cpu time, including 2642 ms collecting 6474 ms elapsed real time, including 2641 ms collecting 2833256224 bytes allocated, including 2762557728 bytes reclaimed Timing modulo (power of 2) (time (do ((...)) ...)) 299 collections 8536 ms elapsed cpu time, including 2586 ms collecting 8545 ms elapsed real time, including 2597 ms collecting 2833185520 bytes allocated, including 2838706016 bytes reclaimed Timing modulo (non-power of 2) (time (do ((...)) ...)) 299 collections 8544 ms elapsed cpu time, including 2601 ms collecting 8548 ms elapsed real time, including 2596 ms collecting 2833185520 bytes allocated, including 2821785344 bytes reclaimedmy OCaml version. Because 15000000 is bigger than the bound for basic array, i had to use Bigarray.
open Bigarray open Sys open Format let timeit f = let tb = time () in let res = Lazy.force_val f in printf "time: %f s\n" ((time()) -. tb); res let ba_fold_left f init arr = let adim = Array1.dim arr in let rec ba_fold_left' i res = match i >= 0 with | true -> ba_fold_left' (i - 1) (f res i (Array1.get arr i)) | _ -> res in ba_fold_left' (adim - 1) init ;; let sieve n = let adim = (n - 1)/ 2 in let arr = Array1.create int8_unsigned c_layout adim in Array1.fill arr 0; let rec erase i mo = if i < adim then begin if (2 * i + 3) mod mo = 0 then Array1.set arr i 1; erase (i + mo) mo end in let rec sieve' i = let e = 2 * i + 3 in match e*e with | p when p > n -> () | p when (Array1.get arr i) = 0 -> erase (p / 2 - 1) e; sieve' (i + 1) | _ -> sieve' (i + 1) in sieve' 0; 2 :: (ba_fold_left (fun x k v -> if v = 0 then (k * 2 + 3) :: x else x) [] arr) ;; let _ = Printf.printf " %d " (List.length (timeit (Lazy.lazy_from_fun (fun () -> (sieve 15485863)))))[…] ancient Sieve of Eratosthenes that computes the list of prime numbers is inefficient in the sense that some composite numbers are […]
This was my first Factor program I wrote a few months ago, to learn the language. Updated it with the optimization to start sieving at the square of the prime being sieved. (Didn’t know about that trick before :)
USING: kernel math sequences math.ranges bit-sets locals formatting io ; FROM: sets => adjoin delete in? members ; IN: sieve CONSTANT: items/line 10 : .row ( seq -- ) [ "%5d " printf ] each nl ; : .vec ( seq -- ) [ dup length items/line > ] [ items/line cut-slice swap .row ] while .row ; : init-sieve ( max-prime -- bitset ) dup <bit-set> 2 rot [a,b) [ over adjoin ] each ; : filter-composite ( bitset start end -- ) [ dup sq ] dip rot <range> [ over delete ] each drop ; :: primes ( maxprime -- vec ) [let maxprime 1 + init-sieve :> sieve 2 [ dup sq maxprime > ] [ dup sieve in? [ dup maxprime sieve -rot filter-composite ] when 1 + ] until drop sieve members ] ; : .primes ( maxprime -- ) primes .vec ;Session:
( scratchpad ) 15485863 primes length . 1000000 ( scratchpad ) 127 .primes 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 ( scratchpad )[…] As mentioned in my last article, I started doing some challenges from the Programming Praxis website. And here comes my PHP solution to the second challenge. […]
I’m working my way through old exercises; here’s my submission.
Python was a bit slow for me, and I’m trying to learn Common Lisp, so here’s
another version. I was only able to figure out how to get 2/3 of the optimizations
working, but it’s still absurdly fast compared to my last submission.
#!/usr/bin/env python import itertools import sys # note: the 'modulo' version with the same cut-offs uses much less memory in Python # but is also way slower def sieve(n): # optimization 1: only odd numbers primes = range(3, n + 1, 2) upper_bound = n ** 0.5 for base in xrange(len(primes)): if not primes[base]: continue # optimization 3: stop at the square root of n if primes[base] >= upper_bound: break # optimization 2: start with the square for i in xrange(base + (base + 1) * primes[base], len(primes), primes[base]): primes[i] = None primes.insert(0, 2) return filter(None, primes) def main(): print 'Enter an upper bound for the sieve.' n = int(sys.stdin.readline()) result = sieve(n) print 'Number of primes: {0}'.format(len(result)) if __name__ == '__main__': main()I just started learning ruby. Here is my first try.
def prime_numbers(n)
return [] if n < 2 # no prime numbers before 2
max = (n**0.5).to_i; # define until when we need to loop
primes = []
# optimization 1: checking only odd numbers
primes = 3.step( n, 2 ).to_a
primes.each do |p|
next unless p
# optimization 3: checking until the root of the maximum number
break unless p <= max
# optimization 2: checking starts at the square
((p**2)/2-1).step( primes.length-1, p) do |i|
primes[i] = nil
end
end
primes.insert(0, 2).delete(nil) # add the prime number 2 and remove false numbers
primes
end
puts prime_numbers(15485863).length
[\sourcecode]
Another clojure version,
because the first one die on me with OutOfMemoryError to compute (eratosthenes 15485863)
; we work only with odd numbers till N (defn odds [n] (vec (range 3 (inc n) 2))) ; it's easy to calculate index of a value in our vector of odds (defn odds-idx [v] (/ (- v 3) 2)) ; just sieve the vector from start=(index of v²) (defn sieve [max-idx nums v] (let [start (odds-idx (* v v))] (reduce #(assoc %1 %2 nil) nums (range start max-idx v)))) ; get primes equal or less than N (defn primes [n] (let [max-i (Math/floor (odds-idx n))] (loop [nums (odds n) i 0] (if-let [v (nums i)] (if (> (* v v) n) (cons 2 (remove nil? nums)) (recur (sieve max-i nums v) (inc i))) (recur nums (inc i)))))) (count (primes 15485863)) 1000000An imperative OCaml Batteries
BitSet-based version, without multiplications, only additions, subtractions and shifts:let sieve_of_eratosthenes (proc : int -> unit) n = let limit = (n - 3) asr 1 in let sieve = BitSet.create_full (limit + 1) in let i = ref 0 and p = ref 3 and q = ref 9 in while !q < n do if BitSet.is_set sieve !i then begin let j = ref ((!q - 3) asr 1) in while !j <= limit do BitSet.unset sieve !j; j := !j + !p done end; incr i; p := !p + 2; q := !q + (!i + 1) lsl 3 done; proc 2; foreach (BitSet.enum sieve) (fun i -> if i <= limit then proc (i lsl 1 + 3))My C++ solution where a “crossing out” means “made 0”
http://ideone.com/fdmux
Here’s a Matlab solution:
function primes=sieve(maxnum) %Sieve of Erasthenes for finding primes less than or equal to the input %argument mask=ones(ceil((maxnum-1)/2),1); for i=3:2:floor(sqrt(maxnum)) if mask((i-1)/2) mask((i^2-1)/2:i:end)=0; end end primes=[2;2*find(mask)+1];It’s a tiny bit faster than the nearly-identical Matlab implementation, called primes.
A rather simple version of a Python solution.
def primes_sieve(limit):
limit += 1
not_prime = [False] * limit
primes = []
for i in xrange(2, limit):
if not_prime[i]:
continue
for j in xrange(i*2, limit, i):
not_prime[j] = True
primes.append(i)
return primes
if __name__ == "__main__":
print len(primes_sieve(15485863))
I tried another solution in C++ using 2 threads, one running upwards on the numbers, the other running downwards. They stop if the reach another.
The solution is correct, but the performance decreases severly.
Here is my code:
#include
#include
#include
#include
#include
using namespace std;
const unsigned long MAX_NUMBER = 15485863;
vector numbers(MAX_NUMBER+1);
pthread_mutex_t mutex1 = PTHREAD_MUTEX_INITIALIZER;
unsigned long l1,l2;
// thread specific arguments
struct thread_data
{
unsigned int thread_id;
unsigned long n; // data to sieve
bool bw; // backward or forward search
};
struct thread_data thread_data_array[2];
void init_numbers(vector * numbers);
void *prime_sieve(void *arg);
void init_numbers(vector * numbers)
{
unsigned long i;
(*numbers)[2]=true;
for (i=3;ibw))
{
while (l1 <= l2)
{
if (numbers[l1] != false)
{
c = l1*l1;
pos_move = (l1 << 1);
while (c n)
{
pthread_mutex_lock(&mutex1);
numbers[c]=false;
pthread_mutex_unlock(&mutex1);
c += pos_move;
}
}
pthread_mutex_lock(&mutex1);
l1++;
pthread_mutex_unlock(&mutex1);
}
}
else
{
while(l2 >= l1)
{
if (numbers[l2] != false)
{
c = l2*l2;
pos_move = (l2 << 1);
while (c n)
{
pthread_mutex_lock(&mutex1);
numbers[c]=false;
pthread_mutex_unlock(&mutex1);
c += pos_move;
}
}
pthread_mutex_lock(&mutex1);
l2–;
pthread_mutex_unlock(&mutex1);
}
}
pthread_exit(NULL);
}
unsigned long get_primes(vector numbers, unsigned long n)
{
unsigned long primes=0;
for(unsigned int i=0;i<=n;i++)
{
primes+=numbers[i];
}
return primes;
}
int main(int argc, char *argv[])
{
pthread_t threads[2];
int rc;
unsigned long n = MAX_NUMBER;
init_numbers(&numbers);
l1 = 2;
l2 = (int)(sqrt(n))+1;
thread_data_array[0].thread_id = 0;
thread_data_array[0].n = MAX_NUMBER;
thread_data_array[0].bw = false;
rc = pthread_create(&threads[0], NULL, prime_sieve, (void *) &thread_data_array[0]);
thread_data_array[1].thread_id = 1;
thread_data_array[1].n = MAX_NUMBER;
thread_data_array[1].bw = true;
rc = pthread_create(&threads[1], NULL, prime_sieve, (void *) &thread_data_array[1]);
void* retval;
pthread_join(threads[0], &retval);
pthread_join(threads[1], &retval);
cout << get_primes(numbers, n) << endl;
}
Using python:
from math import *
def sieve2(n):
index = 0
numberList = range(3,n+1,2)
while numberList[index] = sqrt(n):
break
return [2] + numberList
^
|
|
Didn’t post correctly
https://github.com/gcapell/ProgrammingPraxis/blob/master/02_sieve_of_eratosthenes/sieve.go
package main import "fmt" func PrimesLessEqual(n int) [] int { c, p, w, l := 0, 2, make([] bool, n + 1), make([] int, n) for ; p <= n; p++ { for p <= n && w[p] { p++ } if p <= n { l[c] = p; c++ for k := p; k <= n; k += p { w[k] = true } } } return l[0:c] } func main() { fmt.Println(len(PrimesLessEqual(15485863))) }compact ruby implementation:
list = (3..ARGV[0].to_i).step(2).to_a (0..(Math.sqrt(list.last)/2)-1).each do |i| (((list[i]**2)/2) - 1..list.size).step(list[i]) { |j| list[j] = 0 } if (list[i] > 0) end puts ([2] + list.select {|x| x > 0}).sizeA Scala solution:
import scala.math._ import scala.collection.mutable._ object Main extends Application { val limit = 15485863 //20 val limitSquareRoot = pow(limit,0.5) toInt def sieve(factor: Int, xs: Array[Int]): Array[Int] = { for (i <- (factor * factor) until (limit + 1) by factor) { xs.update((i-2), -1) } xs } def sieveOfEtratosthenes: Array[Int] = { def numbersArray(start: Int, end: Int): Array[Int] = start until end toArray val xs = numbersArray(2, limit + 1) for (i x != -1} } println(sieveOfEtratosthenes.length)// foreach println) } MainPrevious post did not copy properly:
import scala.math._ import scala.collection.mutable._ object Main extends Application { val limit = 15485863 //20 val limitSquareRoot = pow(limit,0.5) toInt def sieve(factor: Int, xs: Array[Int]): Array[Int] = { for (i <- (factor * factor) until (limit + 1) by factor) { xs.update((i-2), -1) } xs } def sieveOfEtratosthenes: Array[Int] = { def numbersArray(start: Int, end: Int): Array[Int] = start until end toArray val xs = numbersArray(2, limit + 1) for (i x != -1} } println(sieveOfEtratosthenes.length)// foreach println) } Main#include
#include
#include
#include
#include
unsigned LIMIT 1000000000;
int main()
{
int all[LIMIT],i,j=2,k=0,lim;
float start,stop;
printf(“Enter the limit(excluding 0)”);
scanf(“%d”,&lim);
for(i=0;i<lim;i++)
{
all[i]=i;
}
all[0]=0;
all[1]=0;
stop=(unsigned)sqrt(lim);
while(j<stop)
{
start=0;
if(all[j]!=0)
{
start=all[j]*all[j];
if(all[start]%all[j]==0)
{
all[start]=0;
start++;
}
}
j++;
}
for(i=2;i<lim-2;i++)
{
if(all[i]!=0)
{
printf(" %d",all[i]);
}
}
getch();
return 0;
}
Corrected solution :) very simple works like a charm
#include “stdafx.h”
#include
#include
#include
#include
#include
int main()
{
int all[100],i,j=2,k=0,lim,start;
float stop;
printf(“Enter the limit(excluding 0)”);
scanf(“%d”,&lim);
lim=lim+1;
for(i=0;i<=lim;i++)
{
all[i]=i;
}
all[0]=0;
all[1]=0;
stop=sqrt((float)lim);
while(j<stop)
{
start=0;
if(all[j]!=0)
{
start=all[j]*all[j];
while(start<=lim)
{
all[start]=0;
start=start+j;
}
}
j++;
}
for(i=2;i<=lim-2;i++)
{
if(all[i]!=0)
{
printf(" %d",all[i]);
}
}
getch();
return 0;
}
now if any of you guys take a look here how can you make this work for a billionth using the constant did not give me result,are linked list only solution???
A memory-hungry Python solution: https://gist.github.com/gmcclure/5475130
C# solution:
public class Primer
{
public int[] GetPrimes(int max)
{
if (max < 0) max = max * -1;
if (max 0)
{
for (int i = start * start; i < primeIndexes.Length; i += start)
{
primeIndexes[i] = true;
}
start++;
}
List primes = new List(max / 10);
for (int i = 1; i primeIndexes.Length)
return 0;
for (int i = minIndex; i < primeIndexes.Length; i++)
{
if (primeIndexes[i] == false)
return i;
}
return 0;
}
}
Just over a second on Core2 Duo.
Hm, it didn’t format properly, here’s a gist:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Primer.cs
hosted with ❤ by GitHub
Hi, just learn python and still beginner level, so I don’t know
def sieve(number):
array = [False,False] + [True] * (number – 1)
result = []
prime_found = 0
index = 2
while index <= number:
if array[index]:
prime_found = index
result.append(index)
index *= index
while index <= number:
array[index] = False
index += prime_found
index = prime_found + 1
else:
index += 1
return result
But I remember in java, an array of boolean is really cheap, only 1 bit per element. and if I switch so False is Boolean while True is not or not examined yet, this will probably safe some space and array creation time.
Sorry, Just disappointed with the way my code shown here, so I just try this to see how to format my code appropriately:
def just_a_test(‘my apology’):
print ‘Again, really sorry’
Java solution:
import java.util.*; class Sieve { public static void main(String[] args) { for (Integer prime : sieveFor(Integer.parseInt(args[0]))) { System.out.println(prime.toString()); } } public static Integer[] sieveFor(int maximum) { return (new Sieve(maximum)).primes(); } private int maximum; public Sieve(int maximum) { this.maximum = maximum; } public Integer[] primes() { TreeSet<Integer> primes = allNumbers(); Integer prime = primes.first(); while (prime <= largestSievableNumber()) { prime = primes.higher(prime); Integer potential_composite = primes.higher(firstUnsievedComposite(prime) - 1); while (potential_composite != null) { if (potential_composite % prime == 0) { primes.remove(potential_composite); } potential_composite = primes.higher(potential_composite); } } return primes.toArray(new Integer[0]); } private TreeSet<Integer> allNumbers() { TreeSet<Integer> numbers = new TreeSet<Integer>(); numbers.add(2); for (int i = 3;i <= this.maximum;i += 2) { numbers.add(i); } return numbers; } private Integer largestSievableNumber() { return new Integer((int) Math.pow(this.maximum, 0.5)); } private ArrayList<Integer> compositesFor(int prime, int maximum) { ArrayList<Integer> composites = new ArrayList<Integer>(); for (int i = firstUnsievedComposite(prime);i <= this.maximum;i += prime) { composites.add(i); } return composites; } private Integer firstUnsievedComposite(int prime) { return new Integer((int) Math.pow(prime, 2)); } }jUnit Tests:
import static org.junit.Assert.*; import org.junit.Test; public class SieveTest { @Test public void even() { assertArrayEquals(new Integer[]{2}, sieveFor(2)); } @Test public void firstOdd() { assertArrayEquals(new Integer[]{2, 3}, sieveFor(3)); assertArrayEquals(new Integer[]{2, 3}, sieveFor(4)); } @Test public void compositeOdd() { assertArrayEquals(new Integer[]{2, 3, 5, 7}, sieveFor(9)); } @Test public void largeList() { assertArrayEquals(new Integer[]{2, 3, 5, 7, 11, 13, 17, 19, 23, 29}, sieveFor(30)); } private Integer[] sieveFor(int maximum) { return Sieve.sieveFor(maximum); } }Output on 15485863 took a while, although I didn’t time it. Does anybody have some performance suggestions, my java’s a bit rusty.
My solution in C using all optimizations. It prints primes as it finds them, and then once it hits the square root mark prints any remaining untouched variables in the array. Chews through N=15485863 in about 3.3 seconds.
In lua, not efficient, but short:
function sieve(n,s)
for i=2,math.ceil(math.sqrt(n)) do
for j=2*i,n,i do s[j]=false end
end
return s
end
s=sieve(15485863,{false,false})
for i=1,15485863 do if s[i]==nil then print(i) end end
Well, it “only” takes only 5.5s for 15485863 (and 262MB RAM)…
function sieve(n,s)
for i=2,math.ceil(math.sqrt(n)) do
if (i==2) or (i%2==1) then for j=i*i,n,(i%2+1)*i do s[j]=false end end
end
return s
end
s=sieve(15485863,{false,false})
for i=1,15485863 do if s[i]==nil then print(i) end end
(a few improvements, now it only takes 2.3s [and 1.5s for printing…])
Just stumbled across this page. I have written a multi-threaded nth-prime algorithm in C# (it finds the millionth prime in 7 milliseconds using 16 threads on a 4 core – 8 hyper-threaded core – i7 3770s running at 3.1 GHz while using only 405 KB for the sieve) and was looking for information on how to convert it to racket. I’m new to racket and still trying to wrap my head around the functional approach to the sieve.
The various examples in functional languages should give me what I was looking for. Thanks to everyone for sharing.
If anyone is interested in the multi-threaded algorithm, let me know (I will be notified of follow up comments). Alternatively, as it appears that this page is all about coding challenges, how about a multi-threaded sieve challenge? I also want to convert my algo to c and compile it with clang to see what kind of boost clang’s optimers might provide vs. compiling with gcc (and vs. my visual studio compiled c#).
Note that my algorithm is an nth-prime algorithm – the inverse of what was requested on this page. In other words, given the input 1000000, my algo outputs 15485863.
There are three reasons I want to convert it to racket:
1. To use racket’s native big integers (the C# version begins to return clipped answers somewhere in the neighborhood of the 240 millionth prime).
2. To see what kind of performance penalty racket’s native big integers incurr (a simple tail recursive factorial in wicked fast in racket – faster than any big integer factorial I can find in c#).
3. It’s fun to learn new languages.
This is my PHP version that calculates half a million before the terminal crashes, don’t have a local copy of PHP to test with at the minute so using a remote low powered box, had to chunk the number sequences as a result of the low performing box.
https://gist.github.com/71138a75c285709df950.git
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
sieve-php
hosted with ❤ by GitHub
In Rust: http://xojoc.pw/programming-praxis/sieve-eratosthenes.html
Python / Numpy solution
import itertools as it import numpy as np def sieve(n): if n > 2: yield 2 if n > 3: limit = int(n ** 0.5) size = (n - 1) // 2 flags = np.ones([size], dtype=bool) primes = it.compress(it.count(3, 2), flags) for p in primes: yield p if p > limit: yield from primes break flags[(p*p - 3) // 2::p] = False len(list(sieve(15485863)))[…] start at p^2 instead of p and the outer loop can stop at the square root of n. I’ll leave the optimized version for you to work […]