Sieve of Eratosthenes
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
The Sieve of Eratosthenes is a simple algorithm that finds the prime numbers up to a given integer. Implement this algorithm, with the only allowed optimization that the outer loop can stop at the square root of the limit, and the inner loop may start at the square of the prime just found. That means especially that you shouldn't optimize by using pre-computed wheels, i.e. don't assume you need only to cross out odd numbers (wheel based on 2), numbers equal to 1 or 5 modulo 6 (wheel based on 2 and 3), or similar wheels based on low primes.
If there's an easy way to add such a wheel based optimization, implement this as an alternative version.
- Note
- It is important that the sieve algorithm be the actual algorithm used to find prime numbers for the task.
- Cf
360 Assembly
For maximum compatibility, this program uses only the basic instruction set. <lang 360_Assembly>* Sieve of Eratosthenes ERATOST CSECT
USING ERATOST,R12
SAVEAREA B STM-SAVEAREA(R15)
DC 17F'0' DC CL8'ERATOST'
STM STM R14,R12,12(R13) save calling context
ST R13,4(R15) ST R15,8(R13) LR R12,R15 set addessability
- ---- CODE
LA R4,1 I=1 LA R6,1 increment L R7,N limit
LOOPI BXH R4,R6,ENDLOOPI do I=2 to N
LR R1,R4 R1=I BCTR R1,0 LA R14,CRIBLE(R1) CLI 0(R14),X'01' BNE ENDIF if not CRIBLE(I) LR R5,R4 J=I LR R8,R4 LR R9,R7
LOOPJ BXH R5,R8,ENDLOOPJ do J=I*2 to N by I
LR R1,R5 R1=J BCTR R1,0 LA R14,CRIBLE(R1) MVI 0(R14),X'00' CRIBLE(J)='0'B B LOOPJ
ENDLOOPJ EQU * ENDIF EQU *
B LOOPI
ENDLOOPI EQU *
LA R4,1 I=1 LA R6,1 L R7,N
LOOP BXH R4,R6,ENDLOOP do I=1 to N
LR R1,R4 R1=I BCTR R1,0 LA R14,CRIBLE(R1) CLI 0(R14),X'01' BNE NOTPRIME if not CRIBLE(I) CVD R4,P P=I UNPK Z,P Z=P MVC C,Z C=Z OI C+L'C-1,X'F0' zap sign MVC WTOBUF(8),C+8 WTO MF=(E,WTOMSG)
NOTPRIME EQU *
B LOOP
ENDLOOP EQU * RETURN EQU *
LM R14,R12,12(R13) restore context XR R15,R15 set return code to 0 BR R14 return to caller
- ---- DATA
I DS F J DS F
DS 0F
P DS PL8 packed Z DS ZL16 zoned C DS CL16 character WTOMSG CNOP 0,4
DC H'80' length of WTO buffer DC H'0' must be binary zeroes
WTOBUF DC 80C' '
LTORG
N DC F'100000' CRIBLE DC 100000X'01'
YREGS END ERATOST</lang>
- Output:
00000002 00000003 00000005 00000007 00000011 00000013 00000017 00000019 00000023 00000029 00000031 00000037 00000041 00000043 00000047 00000053 00000059 00000061 00000067 ... 00099767 00099787 00099793 00099809 00099817 00099823 00099829 00099833 00099839 00099859 00099871 00099877 00099881 00099901 00099907 00099923 00099929 00099961 00099971 00099989 00099991
68000 Assembly
Algorithm somewhat optimized: array omits 1, 2, all higher odd numbers. Optimized for storage: uses bit array for prime/composite flags.
Some of the macro code is derived from the examples included with EASy68K. See 68000 "100 Doors" listing for additional information. <lang 68000devpac>*-----------------------------------------------------------
- Title : BitSieve
- Written by : G. A. Tippery
- Date : 2014-Feb-24, 2013-Dec-22
- Description: Prime number sieve
- -----------------------------------------------------------
ORG $1000
- ---- Generic macros ---- **
PUSH MACRO MOVE.L \1,-(SP) ENDM
POP MACRO MOVE.L (SP)+,\1 ENDM
DROP MACRO ADDQ #4,SP ENDM
PUTS MACRO ** Print a null-terminated string w/o CRLF ** ** Usage: PUTS stringaddress ** Returns with D0, A1 modified MOVEQ #14,D0 ; task number 14 (display null string) LEA \1,A1 ; address of string TRAP #15 ; display it ENDM
GETN MACRO MOVEQ #4,D0 ; Read a number from the keyboard into D1.L. TRAP #15 ENDM
- ---- Application-specific macros ---- **
val MACRO ; Used by bit sieve. Converts bit address to the number it represents. ADD.L \1,\1 ; double it because odd numbers are omitted ADDQ #3,\1 ; add offset because initial primes (1, 2) are omitted ENDM
- ** ================================================================================ **
- ** Integer square root routine, bisection method **
- ** IN: D0, should be 0<D0<$1000 (65536) -- higher values MAY work, no guarantee
- ** OUT: D1
SquareRoot:
MOVEM.L D2-D4,-(SP) ; save registers needed for local variables
- DO == n
- D1 == a
- D2 == b
- D3 == guess
- D4 == temp
- a = 1;
- b = n;
MOVEQ #1,D1 MOVE.L D0,D2
- do {
REPEAT
- guess = (a+b)/2;
MOVE.L D1,D3 ADD.L D2,D3 LSR.L #1,D3
- if (guess*guess > n) { // inverse function of sqrt is square
MOVE.L D3,D4 MULU D4,D4 ; guess^2 CMP.L D0,D4 BLS .else
- b = guess;
MOVE.L D3,D2 BRA .endif
- } else {
.else:
- a = guess;
MOVE.L D3,D1
- } //if
.endif:
- } while ((b-a) > 1); ; Same as until (b-a)<=1 or until (a-b)>=1
MOVE.L D2,D4 SUB.L D1,D4 ; b-a UNTIL.L D4 <LE> #1 DO.S
- return (a) ; Result is in D1
- } //LongSqrt()
MOVEM.L (SP)+,D2-D4 ; restore saved registers RTS
- ** ================================================================================ **
- ======================================================================= **
-
- Prime-number Sieve of Eratosthenes routine using a big bit field for flags **
- Enter with D0 = size of sieve (bit array)
- Prints found primes 10 per line
- Returns # prime found in D6
- Register usage:
- D0 == n
- D1 == prime
- D2 == sqroot
- D3 == PIndex
- D4 == CIndex
- D5 == MaxIndex
- D6 == PCount
- A0 == PMtx[0]
- On return, all registers above except D0 are modified. Could add MOVEMs to save and restore D2-D6/A0.
- ------------------------ **
GetBit: ** sub-part of Sieve subroutine ** ** Entry: bit # is on TOS ** Exit: A6 holds the byte number, D7 holds the bit number within the byte ** Note: Input param is still on TOS after return. Could have passed via a register, but
** wanted to practice with stack. :)
MOVE.L (4,SP),D7 ; get value from (pre-call) TOS ASR.L #3,D7 ; /8 MOVEA D7,A6 ; byte # MOVE.L (4,SP),D7 ; get value from (pre-call) TOS AND.L #$7,D7 ; bit # RTS
- ------------------------ **
Sieve: MOVE D0,D5 SUBQ #1,D5 JSR SquareRoot ; sqrt D0 => D1 MOVE.L D1,D2 LEA PArray,A0 CLR.L D3
PrimeLoop: MOVE.L D3,D1 val D1 MOVE.L D3,D4 ADD.L D1,D4
CxLoop: ; Goes through array marking multiples of d1 as composite numbers CMP.L D5,D4 BHI ExitCx PUSH D4 ; set D7 as bit # and A6 as byte pointer for D4'th bit of array JSR GetBit DROP BSET D7,0(A0,A6.L) ; set bit to mark as composite number ADD.L D1,D4 ; next number to mark BRA CxLoop ExitCx: CLR.L D1 ; Clear new-prime-found flag ADDQ #1,D3 ; Start just past last prime found PxLoop: ; Searches for next unmarked (not composite) number CMP.L D2,D3 ; no point searching past where first unmarked multiple would be past end of array BHI ExitPx ; if past end of array TST.L D1 BNE ExitPx ; if flag set, new prime found PUSH D3 ; check D3'th bit flag JSR GetBit ; sets D7 as bit # and A6 as byte pointer DROP ; drop TOS BTST D7,0(A0,A6.L) ; read bit flag BNE IsSet ; If already tagged as composite MOVEQ #-1,D1 ; Set flag that we've found a new prime IsSet: ADDQ #1,D3 ; next PIndex BRA PxLoop ExitPx: SUBQ #1,D3 ; back up PIndex TST.L D1 ; Did we find a new prime #? BNE PrimeLoop ; If another prime # found, go process it
; fall through to print routine
- ------------------------ **
- Print primes found
- D4 == Column count
- Print header and assumed primes (#1, #2)
PUTS Header ; Print string @ Header, no CR/LF
MOVEQ #2,D6 ; Start counter at 2 because #1 and #2 are assumed primes MOVEQ #2,D4
MOVEQ #0,D3 PrintLoop: CMP.L D5,D3 BHS ExitPL PUSH D3 JSR GetBit ; sets D7 as bit # and A6 as byte pointer DROP ; drop TOS BTST D7,0(A0,A6.L) BNE NotPrime
- printf(" %6d", val(PIndex)
MOVE.L D3,D1 val D1 AND.L #$0000FFFF,D1 MOVEQ #6,D2 MOVEQ #20,D0 ; display signed RJ TRAP #15 ADDQ #1,D4 ADDQ #1,D6
- *** Display formatting ***
- if((PCount % 10) == 0) printf("\n");
CMP #10,D4 BLO NoLF PUTS CRLF MOVEQ #0,D4 NoLF: NotPrime: ADDQ #1,D3 BRA PrintLoop ExitPL: RTS
- ======================================================================= **
N EQU 5000 ; *** Size of boolean (bit) array *** SizeInBytes EQU (N+7)/8
START: ; first instruction of program MOVE.L #N,D0 ; # to test JSR Sieve
- printf("\n %d prime numbers found.\n", D6); ***
PUTS Summary1,A1 MOVE #3,D0 ; Display signed number in D1.L in decimal in smallest field. MOVE.W D6,D1 TRAP #15 PUTS Summary2,A1
SIMHALT ; halt simulator
- ======================================================================= **
- Variables and constants here
ORG $2000 CR EQU 13 LF EQU 10 CRLF DC.B CR,LF,$00
PArray: DCB.B SizeInBytes,0
Header: DC.B CR,LF,LF,' Primes',CR,LF,' ======',CR,LF DC.B ' 1 2',$00
Summary1: DC.B CR,LF,' ',$00 Summary2: DC.B ' prime numbers found.',CR,LF,$00
END START ; last line of source</lang>
ABAP
<lang Lisp> PARAMETERS: p_limit TYPE i OBLIGATORY DEFAULT 100.
AT SELECTION-SCREEN ON p_limit.
IF p_limit LE 1. MESSAGE 'Limit must be higher then 1.' TYPE 'E'. ENDIF.
START-OF-SELECTION.
FIELD-SYMBOLS: <fs_prime> TYPE flag. DATA: gt_prime TYPE TABLE OF flag, gv_prime TYPE flag, gv_i TYPE i, gv_j TYPE i.
DO p_limit TIMES. IF sy-index > 1. gv_prime = abap_true. ELSE. gv_prime = abap_false. ENDIF.
APPEND gv_prime TO gt_prime. ENDDO.
gv_i = 2. WHILE ( gv_i <= trunc( sqrt( p_limit ) ) ). IF ( gt_prime[ gv_i ] EQ abap_true ). gv_j = gv_i ** 2. WHILE ( gv_j <= p_limit ). gt_prime[ gv_j ] = abap_false. gv_j = ( gv_i ** 2 ) + ( sy-index * gv_i ). ENDWHILE. ENDIF. gv_i = gv_i + 1. ENDWHILE.
LOOP AT gt_prime INTO gv_prime. IF gv_prime = abap_true. WRITE: / sy-tabix. ENDIF. ENDLOOP.
</lang>
ACL2
<lang Lisp>(defun nats-to-from (n i)
(declare (xargs :measure (nfix (- n i)))) (if (zp (- n i)) nil (cons i (nats-to-from n (+ i 1)))))
(defun remove-multiples-up-to-r (factor limit xs i)
(declare (xargs :measure (nfix (- limit i)))) (if (or (> i limit) (zp (- limit i)) (zp factor)) xs (remove-multiples-up-to-r factor limit (remove i xs) (+ i factor))))
(defun remove-multiples-up-to (factor limit xs)
(remove-multiples-up-to-r factor limit xs (* factor 2)))
(defun sieve-r (factor limit)
(declare (xargs :measure (nfix (- limit factor)))) (if (zp (- limit factor)) (nats-to-from limit 2) (remove-multiples-up-to factor (+ limit 1) (sieve-r (1+ factor) limit))))
(defun sieve (limit)
(sieve-r 2 limit))</lang>
Ada
<lang Ada>with Ada.Text_IO, Ada.Command_Line;
procedure Eratos is
Last: Positive := Positive'Value(Ada.Command_Line.Argument(1)); Prime: array(1 .. Last) of Boolean := (1 => False, others => True); Base: Positive := 2; Cnt: Positive;
begin
loop exit when Base * Base > Last; if Prime(Base) then Cnt := Base + Base; loop exit when Cnt > Last; Prime(Cnt) := False; Cnt := Cnt + Base; end loop; end if; Base := Base + 1; end loop; Ada.Text_IO.Put("Primes less or equal" & Positive'Image(Last) &" are:"); for Number in Prime'Range loop if Prime(Number) then Ada.Text_IO.Put(Positive'Image(Number)); end if; end loop;
end Eratos;</lang>
- Output:
> ./eratos 31 Primes less or equal 31 are : 2 3 5 7 11 13 17 19 23 29 31
ALGOL 60
Works with: ALGOL 60 for OS/360 <lang algol60>'BEGIN'
'INTEGER' 'ARRAY' CANDIDATES(/0..1000/); 'INTEGER' I,J,K; 'COMMENT' SET LINE-LENGTH=120,SET LINES-PER-PAGE=62,OPEN; SYSACT(1,6,120); SYSACT(1,8,62); SYSACT(1,12,1); 'FOR' I := 0 'STEP' 1 'UNTIL' 1000 'DO' 'BEGIN' CANDIDATES(/I/) := 1; 'END'; CANDIDATES(/0/) := 0; CANDIDATES(/1/) := 0; I := 0; 'FOR' I := I 'WHILE' I 'LESS' 1000 'DO' 'BEGIN' 'FOR' I := I 'WHILE' I 'LESS' 1000 'AND' CANDIDATES(/I/) 'EQUAL' 0 'DO' I := I+1; 'IF' I 'LESS' 1000 'THEN' 'BEGIN' J := 2; K := J*I; 'FOR' K := K 'WHILE' K 'LESS' 1000 'DO' 'BEGIN' CANDIDATES(/K/) := 0; J := J + 1; K := J*I; 'END'; I := I+1; 'END' 'END'; 'FOR' I := 0 'STEP' 1 'UNTIL' 999 'DO' 'IF' CANDIDATES(/I/) 'NOTEQUAL' 0 'THEN' 'BEGIN' OUTINTEGER(1,I); OUTSTRING(1,'(' IS PRIME')'); 'COMMENT' NEW LINE; SYSACT(1,14,1) 'END' 'END'
'END'</lang>
ALGOL 68
<lang algol68>BOOL prime = TRUE, non prime = FALSE; PROC eratosthenes = (INT n)[]BOOL: (
[n]BOOL sieve; FOR i TO UPB sieve DO sieve[i] := prime OD; INT m = ENTIER sqrt(n); sieve[1] := non prime; FOR i FROM 2 TO m DO IF sieve[i] = prime THEN FOR j FROM i*i BY i TO n DO sieve[j] := non prime OD FI OD; sieve
);
print((eratosthenes(80),new line))</lang>
- Output:
FTTFTFTFFFTFTFFFTFTFFFTFFFFFTFTFFFFFTFFFTFTFFFTFFFFFTFFFFFTFTFFFFFTFFFTFTFFFFFTF
ALGOL W
<lang algolw>begin
% implements the sieve of Eratosthenes % % s(i) is set to true if i is prime, false otherwise % % algol W doesn't have a upb operator, so we pass the size of the % % array in n % procedure sieve( logical array s ( * ); integer value n ) ; begin
% start with everything flagged as prime % for i := 1 until n do s( i ) := true;
% sieve out the non-primes % s( 1 ) := false; for i := 2 until truncate( sqrt( n ) ) do begin if s( i ) then begin for p := i * i step i until n do s( p ) := false end if_s_i end for_i ;
end sieve ;
% test the sieve procedure %
integer sieveMax;
sieveMax := 100; begin
logical array s ( 1 :: sieveMax );
i_w := 2; % set output field width % s_w := 1; % and output separator width %
% find and display the primes % sieve( s, sieveMax ); for i := 1 until sieveMax do if s( i ) then writeon( i );
end
end.</lang>
- Output:
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
AutoHotkey
Search autohotkey.com: of Eratosthenes
Source: AutoHotkey forum by Laszlo
<lang autohotkey>MsgBox % "12345678901234567890`n" Sieve(20)
Sieve(n) { ; Sieve of Eratosthenes => string of 0|1 chars, 1 at position k: k is prime
Static zero := 48, one := 49 ; Asc("0"), Asc("1") VarSetCapacity(S,n,one) NumPut(zero,S,0,"char") i := 2 Loop % sqrt(n)-1 { If (NumGet(S,i-1,"char") = one) Loop % n//i If (A_Index > 1) NumPut(zero,S,A_Index*i-1,"char") i += 1+(i>2) } Return S
}</lang>
AutoIt
<lang autoit>#include <Array.au3> $M = InputBox("Integer", "Enter biggest Integer") Global $a[$M], $r[$M], $c = 1 For $i = 2 To $M -1 If Not $a[$i] Then $r[$c] = $i $c += 1 For $k = $i To $M -1 Step $i $a[$k] = True Next EndIf Next $r[0] = $c - 1 ReDim $r[$c] _ArrayDisplay($r)</lang>
AWK
An initial array holds all numbers 2..max (which is entered on stdin); then all products of integers are deleted from it; the remaining are displayed in the unsorted appearance of a hash table. Here, the script is entered directly on the commandline, and input entered on stdin:
$ awk '{for(i=2;i<=$1;i++) a[i]=1; > for(i=2;i<=sqrt($1);i++) for(j=2;j<=$1;j++) delete a[i*j]; > for(i in a) printf i" "}' 100 71 53 17 5 73 37 19 83 47 29 7 67 59 11 97 79 89 31 13 41 23 2 61 43 3
The following variant does not unset non-primes, but sets them to 0, to preserve order in output:
$ awk '{for(i=2;i<=$1;i++) a[i]=1; > for(i=2;i<=sqrt($1);i++) for(j=2;j<=$1;j++) a[i*j]=0; > for(i=2;i<=$1;i++) if(a[i])printf i" "}' 100 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
Now with the script from a file, input from commandline as well as stdin, and input is checked for valid numbers: <lang awk>
- usage: gawk -v n=101 -f sieve.awk
function sieve(n) { # print n,":" for(i=2; i<=n; i++) a[i]=1; for(i=2; i<=sqrt(n);i++) for(j=2;j<=n;j++) a[i*j]=0; for(i=2; i<=n; i++) if(a[i]) printf i" " print "" }
BEGIN { print "Sieve of Eratosthenes:" if(n>1) sieve(n) }
{ n=$1+0 } n<2 { exit } { sieve(n) }
END { print "Bye!" } </lang>
BASIC
<lang freebasic>DIM n AS Integer, k AS Integer, limit AS Integer
INPUT "Enter number to search to: "; limit DIM flags(limit) AS Integer
FOR n = 2 TO SQR(limit)
IF flags(n) = 0 THEN FOR k = n*n TO limit STEP n flags(k) = 1 NEXT k END IF
NEXT n
' Display the primes FOR n = 2 TO limit
IF flags(n) = 0 THEN PRINT n; ", ";
NEXT n</lang>
Applesoft BASIC
<lang basic>10 INPUT "ENTER NUMBER TO SEARCH TO: ";LIMIT 20 DIM FLAGS(LIMIT) 30 FOR N = 2 TO SQR (LIMIT) 40 IF FLAGS(N) < > 0 GOTO 80 50 FOR K = N * N TO LIMIT STEP N 60 FLAGS(K) = 1 70 NEXT K 80 NEXT N 90 REM DISPLAY THE PRIMES 100 FOR N = 2 TO LIMIT 110 IF FLAGS(N) = 0 THEN PRINT N;", "; 120 NEXT N</lang>
Locomotive Basic
<lang locobasic>10 DEFINT a-z 20 INPUT "Limit";limit 30 DIM f(limit) 40 FOR n=2 TO SQR(limit) 50 IF f(n)=1 THEN 90 60 FOR k=n*n TO limit STEP n 70 f(k)=1 80 NEXT k 90 NEXT n 100 FOR n=2 TO limit 110 IF f(n)=0 THEN PRINT n;","; 120 NEXT</lang>
ZX Spectrum Basic
<lang basic>10 INPUT "Enter number to search to: ";l 20 DIM p(l) 30 FOR n=2 TO SQR l 40 IF p(n)<>0 THEN NEXT n 50 FOR k=n*n TO l STEP n 60 LET p(k)=1 70 NEXT k 80 NEXT n 90 REM Display the primes 100 FOR n=2 TO l 110 IF p(n)=0 THEN PRINT n;", "; 120 NEXT n</lang>
BBC BASIC
<lang bbcbasic> limit% = 100000
DIM sieve% limit% prime% = 2 WHILE prime%^2 < limit% FOR I% = prime%*2 TO limit% STEP prime% sieve%?I% = 1 NEXT REPEAT prime% += 1 : UNTIL sieve%?prime%=0 ENDWHILE REM Display the primes: FOR I% = 1 TO limit% IF sieve%?I% = 0 PRINT I%; NEXT</lang>
bash
See solutions at UNIX Shell.
Batch File
<lang dos>:: Sieve of Eratosthenes for Rosetta Code - PG @echo off setlocal ENABLEDELAYEDEXPANSION setlocal ENABLEEXTENSIONS rem echo on set /p n=limit: rem set n=100 for /L %%i in (1,1,%n%) do set crible.%%i=1 for /L %%i in (2,1,%n%) do (
if !crible.%%i! EQU 1 ( set /A w = %%i * 2 for /L %%j in (!w!,%%i,%n%) do (
set crible.%%j=0 )
)
) for /L %%i in (2,1,%n%) do (
if !crible.%%i! EQU 1 echo %%i
) pause</lang>
- Output:
limit: 100 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
Befunge
2>:3g" "-!v\ g30 < |!`"O":+1_:.:03p>03g+:"O"`| @ ^ p3\" ":< 2 234567890123456789012345678901234567890123456789012345678901234567890123456789
Bracmat
This solution does not use an array. Instead, numbers themselves are used as variables. The numbers that are not prime are set (to the silly value "nonprime"). Finally all numbers up to the limit are tested for being initialised. The uninitialised (unset) ones must be the primes. <lang bracmat>( ( eratosthenes
= n j i . !arg:?n & 1:?i & whl ' ( (1+!i:?i)^2:~>!n:?j & ( !!i | whl ' ( !j:~>!n & nonprime:?!j & !j+!i:?j ) ) ) & 1:?i & whl ' ( 1+!i:~>!n:?i & (!!i|put$(!i " ")) ) )
& eratosthenes$100 )</lang>
- Output:
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
C
Plain sieve, without any optimizations:<lang c>#include <stdlib.h>
- include <math.h>
char* eratosthenes(int n, int *c) { char* sieve; int i, j, m;
if(n < 2) return NULL;
*c = n-1; /* primes count */ m = (int) sqrt((double) n);
/* calloc initializes to zero */ sieve = calloc(n+1,sizeof(char)); sieve[0] = 1; sieve[1] = 1; for(i = 2; i <= m; i++) if(!sieve[i]) for (j = i*i; j <= n; j += i) if(!sieve[j]){ sieve[j] = 1; --(*c); }
return sieve;
}</lang>Possible optimizations include sieving only odd numbers (or more complex wheels), packing the sieve into bits to improve locality (and allow larger sieves), etc.
Another example:
We first fill ones into an array and assume all numbers are prime. Then, in a loop, fill zeroes into those places where i * j is less than or equal to n (number of primes requested), which means they have multiples! To understand this better, look at the output of the following example. To print this back, we look for ones in the array and only print those spots. <lang C>#include <stdio.h>
- include <malloc.h>
void sieve(int *, int);
int main(int argc, char *argv) {
int *array, n=10; array =(int *)malloc(sizeof(int)); sieve(array,n); return 0;
}
void sieve(int *a, int n) {
int i=0, j=0;
for(i=2; i<=n; i++) { a[i] = 1; }
for(i=2; i<=n; i++) { printf("\ni:%d", i); if(a[i] == 1) { for(j=i; (i*j)<=n; j++) { printf ("\nj:%d", j); printf("\nBefore a[%d*%d]: %d", i, j, a[i*j]); a[(i*j)] = 0; printf("\nAfter a[%d*%d]: %d", i, j, a[i*j]); } } }
printf("\nPrimes numbers from 1 to %d are : ", n); for(i=2; i<=n; i++) { if(a[i] == 1) printf("%d, ", i); } printf("\n\n");
}</lang>
- Output:
<lang Shell>i
j:2 Before a[2*2]: 1 After a[2*2]: 0 j:3 Before a[2*3]: 1 After a[2*3]: 0 j:4 Before a[2*4]: 1 After a[2*4]: 0 j:5 Before a[2*5]: 1 After a[2*5]: 0 i:3 j:3 Before a[3*3]: 1 After a[3*3]: 0 i:4 i:5 i:6 i:7 i:8 i:9 i:10 Primes numbers from 1 to 10 are : 2, 3, 5, 7, </lang>
C++
<lang cpp>// yield all prime numbers less than limit. template<class UnaryFunction> void primesupto(int limit, UnaryFunction yield) {
std::vector<bool> is_prime(limit, true); const int sqrt_limit = static_cast<int>(std::sqrt(limit)); for (int n = 2; n <= sqrt_limit; ++n) if (is_prime[n]) {
yield(n);
for (unsigned k = n*n, ulim = static_cast<unsigned>(limit); k < ulim; k += n)
//NOTE: "unsigned" is used to avoid an overflow in `k+=n` for `limit` near INT_MAX
is_prime[k] = false;
}
for (int n = sqrt_limit + 1; n < limit; ++n) if (is_prime[n])
yield(n); }</lang>
Full program:
<lang cpp>/**
$ g++ -I/path/to/boost sieve.cpp -o sieve && sieve 10000000 */
- include <inttypes.h> // uintmax_t
- include <limits>
- include <cmath>
- include <iostream>
- include <sstream>
- include <vector>
- include <boost/lambda/lambda.hpp>
int main(int argc, char *argv[])
{
using namespace std; using namespace boost::lambda;
int limit = 10000; if (argc == 2) { stringstream ss(argv[--argc]); ss >> limit;
if (limit < 1 or ss.fail()) { cerr << "USAGE:\n sieve LIMIT\n\nwhere LIMIT in the range [1, "
<< numeric_limits<int>::max() << ")" << endl;
return 2; } }
// print primes less then 100 primesupto(100, cout << _1 << " "); cout << endl;
// find number of primes less then limit and their sum int count = 0; uintmax_t sum = 0; primesupto(limit, (var(sum) += _1, var(count) += 1));
cout << "limit sum pi(n)\n" << limit << " " << sum << " " << count << endl;
}</lang>
C#
<lang csharp>using System; using System.Collections; using System.Collections.Generic;
namespace SieveOfEratosthenes {
class Program { static void Main(string[] args) { int maxprime = int.Parse(args[0]); var primelist = GetAllPrimesLessThan(maxprime); foreach (int prime in primelist) { Console.WriteLine(prime); } Console.WriteLine("Count = " + primelist.Count); Console.ReadLine(); }
private static List<int> GetAllPrimesLessThan(int maxPrime) { var primes = new List<int>() { 2 }; var maxSquareRoot = Math.Sqrt(maxPrime); var eliminated = new BitArray(maxPrime + 1);
for (int i = 3; i <= maxPrime; i += 2) { if (!eliminated[i]) { primes.Add(i); if (i < maxSquareRoot) { for (int j = i * i; j <= maxPrime; j += 2 * i) { eliminated[j] = true; } } } } return primes; } }
}</lang>
Unbounded
Richard Bird Sieve
To show that C# code can be written in somewhat functional paradigms, the following in an implementation of the Richard Bird sieve from the Epilogue of Melissa E. O'Neill's definitive article in Haskell: <lang csharp>using System; using System.Collections; using System.Collections.Generic; using System.Linq; using PrimeT = UInt32;
class PrimesBird : IEnumerable<PrimeT> { private struct CIS<T> { public T v; public Func<CIS<T>> cont; public CIS(T v, Func<CIS<T>> cont) { this.v = v; this.cont = cont; } } private CIS<PrimeT> pmlts(PrimeT p) { Func<PrimeT, CIS<PrimeT>> fn = null; fn = (c) => new CIS<PrimeT>(c, () => fn(c + p)); return fn(p * p); } private CIS<CIS<PrimeT>> allmlts(CIS<PrimeT> ps) { return new CIS<CIS<PrimeT>>(pmlts(ps.v), () => allmlts(ps.cont())); } private CIS<PrimeT> merge(CIS<PrimeT> xs, CIS<PrimeT> ys) { var x = xs.v; var y = ys.v; if (x < y) return new CIS<PrimeT>(x, () => merge(xs.cont(), ys)); else if (y < x) return new CIS<PrimeT>(y, () => merge(xs, ys.cont())); else return new CIS<PrimeT>(x, () => merge(xs.cont(), ys.cont())); } private CIS<PrimeT> cmpsts(CIS<CIS<PrimeT>> css) { return new CIS<PrimeT>(css.v.v, () => merge(css.v.cont(), cmpsts(css.cont()))); } private CIS<PrimeT> minusat(PrimeT n, CIS<PrimeT> cs) { var nn = n; var ncs = cs; for (; ; ++nn) { if (nn >= ncs.v) ncs = ncs.cont(); else return new CIS<PrimeT>(nn, () => minusat(++nn, ncs)); } } private CIS<PrimeT> prms() { return new CIS<PrimeT>(2, () => minusat(3, cmpsts(allmlts(prms())))); } public IEnumerator<PrimeT> GetEnumerator() { for (var ps = prms(); ; ps = ps.cont()) yield return ps.v; } IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); } }</lang>
Tree Folding Sieve
The above code can easily be converted to "odds-only" and a infinite tree-like folding scheme with the following minor changes: <lang csharp>using System; using System.Collections; using System.Collections.Generic; using System.Linq; using PrimeT = UInt32;
class PrimesTreeFold : IEnumerable<PrimeT> { private struct CIS<T> { public T v; public Func<CIS<T>> cont; public CIS(T v, Func<CIS<T>> cont) { this.v = v; this.cont = cont; } } private CIS<PrimeT> pmlts(PrimeT p) { var adv = p + p; Func<PrimeT, CIS<PrimeT>> fn = null; fn = (c) => new CIS<PrimeT>(c, () => fn(c + adv)); return fn(p * p); } private CIS<CIS<PrimeT>> allmlts(CIS<PrimeT> ps) { return new CIS<CIS<PrimeT>>(pmlts(ps.v), () => allmlts(ps.cont())); } private CIS<PrimeT> merge(CIS<PrimeT> xs, CIS<PrimeT> ys) { var x = xs.v; var y = ys.v; if (x < y) return new CIS<PrimeT>(x, () => merge(xs.cont(), ys)); else if (y < x) return new CIS<PrimeT>(y, () => merge(xs, ys.cont())); else return new CIS<PrimeT>(x, () => merge(xs.cont(), ys.cont())); } private CIS<CIS<PrimeT>> pairs(CIS<CIS<PrimeT>> css) { var nxtcss = css.cont(); return new CIS<CIS<PrimeT>>(merge(css.v, nxtcss.v), () => pairs(nxtcss.cont())); } private CIS<PrimeT> cmpsts(CIS<CIS<PrimeT>> css) { return new CIS<PrimeT>(css.v.v, () => merge(css.v.cont(), cmpsts(pairs(css.cont())))); } private CIS<PrimeT> minusat(PrimeT n, CIS<PrimeT> cs) { var nn = n; var ncs = cs; for (; ; nn += 2) { if (nn >= ncs.v) ncs = ncs.cont(); else return new CIS<PrimeT>(nn, () => minusat(nn + 2, ncs)); } } private CIS<PrimeT> oddprms() { return new CIS<PrimeT>(3, () => minusat(5, cmpsts(allmlts(oddprms())))); } public IEnumerator<PrimeT> GetEnumerator() { yield return 2; for (var ps = oddprms(); ; ps = ps.cont()) yield return ps.v; } IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); } }</lang>
The above code runs over ten times faster than the original Richard Bird algorithm.
Priority Queue Sieve
First, an implementation of a Min Heap Priority Queue is provided as extracted from the entry at RosettaCode, with only the necessary methods duplicated here: <lang csharp>namespace PriorityQ {
using KeyT = UInt32; using System; using System.Collections.Generic; using System.Linq; class Tuple<K, V> { // for DotNet 3.5 without Tuple's public K Item1; public V Item2; public Tuple(K k, V v) { Item1 = k; Item2 = v; } public override string ToString() { return "(" + Item1.ToString() + ", " + Item2.ToString() + ")"; } } class MinHeapPQ<V> { private struct HeapEntry { public KeyT k; public V v; public HeapEntry(KeyT k, V v) { this.k = k; this.v = v; } } private List<HeapEntry> pq; private MinHeapPQ() { this.pq = new List<HeapEntry>(); } private bool mt { get { return pq.Count == 0; } } private Tuple<KeyT, V> pkmn { get { if (pq.Count == 0) return null; else { var mn = pq[0]; return new Tuple<KeyT, V>(mn.k, mn.v); } } } private void psh(KeyT k, V v) { // add extra very high item if none if (pq.Count == 0) pq.Add(new HeapEntry(UInt32.MaxValue, v)); var i = pq.Count; pq.Add(pq[i - 1]); // copy bottom item... for (var ni = i >> 1; ni > 0; i >>= 1, ni >>= 1) { var t = pq[ni - 1]; if (t.k > k) pq[i - 1] = t; else break; } pq[i - 1] = new HeapEntry(k, v); } private void siftdown(KeyT k, V v, int ndx) { var cnt = pq.Count - 1; var i = ndx; for (var ni = i + i + 1; ni < cnt; ni = ni + ni + 1) { var oi = i; var lk = pq[ni].k; var rk = pq[ni + 1].k; var nk = k; if (k > lk) { i = ni; nk = lk; } if (nk > rk) { ni += 1; i = ni; } if (i != oi) pq[oi] = pq[i]; else break; } pq[i] = new HeapEntry(k, v); } private void rplcmin(KeyT k, V v) { if (pq.Count > 1) siftdown(k, v, 0); } public static MinHeapPQ<V> empty { get { return new MinHeapPQ<V>(); } } public static Tuple<KeyT, V> peekMin(MinHeapPQ<V> pq) { return pq.pkmn; } public static MinHeapPQ<V> push(KeyT k, V v, MinHeapPQ<V> pq) { pq.psh(k, v); return pq; } public static MinHeapPQ<V> replaceMin(KeyT k, V v, MinHeapPQ<V> pq) { pq.rplcmin(k, v); return pq; }
}</lang>
The following code implements an improved version of the odds-only O'Neil algorithm, which provides the improvements of only adding base prime composite number streams to the queue when the sieved number reaches the square of the base prime (saving a huge amount of memory and considerable execution time, including not needing as wide a range of a type for the internal prime numbers) as well as minimizing stream processing using fusion: <lang csharp>using System; using System.Collections; using System.Collections.Generic; using System.Linq; using PrimeT = UInt32;
class PrimesPQ : IEnumerable<PrimeT> { private IEnumerator<PrimeT> nmrtr() { MinHeapPQ<PrimeT> pq = MinHeapPQ<PrimeT>.empty; PrimeT bp = 3; PrimeT q = 9; IEnumerator<PrimeT> bps = null; yield return 2; yield return 3; for (var n = (PrimeT)5; ; n += 2) { if (n >= q) { // always equal or less... if (q <= 9) { bps = nmrtr(); bps.MoveNext(); bps.MoveNext(); } // move to 3... bps.MoveNext(); var nbp = bps.Current; q = nbp * nbp; var adv = bp + bp; bp = nbp; pq = MinHeapPQ<PrimeT>.push(n + adv, adv, pq); } else { var pk = MinHeapPQ<PrimeT>.peekMin(pq); var ck = (pk == null) ? q : pk.Item1; if (n >= ck) { do { var adv = pk.Item2; pq = MinHeapPQ<PrimeT>.replaceMin(ck + adv, adv, pq); pk = MinHeapPQ<PrimeT>.peekMin(pq); ck = pk.Item1; } while (n >= ck); } else yield return n; } } } public IEnumerator<PrimeT> GetEnumerator() { return nmrtr(); } IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); } }</lang>
The above code is at least about 2.5 times faster than the Tree Folding version.
Dictionary (Hash table) Sieve
The above code adds quite a bit of overhead in having to provide a version of a Priority Queue for little advantage over a Dictionary (hash table based) version as per the code below: <lang csharp>using System; using System.Collections; using System.Collections.Generic; using System.Linq; using PrimeT = UInt32;
class PrimesDict : IEnumerable<PrimeT> { private IEnumerator<PrimeT> nmrtr() { Dictionary<PrimeT, PrimeT> dct = new Dictionary<PrimeT, PrimeT>(); PrimeT bp = 3; PrimeT q = 9; IEnumerator<PrimeT> bps = null; yield return 2; yield return 3; for (var n = (PrimeT)5; ; n += 2) { if (n >= q) { // always equal or less... if (q <= 9) { bps = nmrtr(); bps.MoveNext(); bps.MoveNext(); } // move to 3... bps.MoveNext(); var nbp = bps.Current; q = nbp * nbp; var adv = bp + bp; bp = nbp; dct.Add(n + adv, adv); } else { if (dct.ContainsKey(n)) { PrimeT nadv; dct.TryGetValue(n, out nadv); dct.Remove(n); var nc = n + nadv; while (dct.ContainsKey(nc)) nc += nadv; dct.Add(nc, nadv); } else yield return n; } } } public IEnumerator<PrimeT> GetEnumerator() { return nmrtr(); } IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); } }</lang>
The above code runs in about three quarters of the time as the above Priority Queue based version for a range of a million primes which will fall even further behind for increasing ranges due to the Dictionary providing O(1) access times as compared to the O(log n) access times for the Priority Queue; the only slight advantage of the PQ based version is at very small ranges where the constant factor overhead of computing the table hashes becomes greater than the "log n" factor for small "n".
Page Segmented Array Sieve
All of the above unbounded versions are really just an intellectual exercise as with very little extra lines of code above the fastest Dictionary based version, one can have an bit-packed page-segmented array based version as follows: <lang csharp>using System; using System.Collections; using System.Collections.Generic; using System.Linq; using PrimeT = UInt32;
class PrimesPgd : IEnumerable<PrimeT> { private const int PGSZ = 1 << 14; // L1 CPU cache size in bytes private const int BFBTS = PGSZ * 8; // in bits private const int BFRNG = BFBTS * 2; public IEnumerator<PrimeT> nmrtr() { IEnumerator<PrimeT> bps = null; List<uint> bpa = new List<uint>(); uint[] cbuf = new uint[PGSZ / 4]; // 4 byte words yield return 2; for (var lowi = (PrimeT)0; ; lowi += BFBTS) { for (var bi = 0; ; ++bi) { if (bi < 1) { if (bi < 0) { bi = 0; yield return 2; } PrimeT nxt = 3 + lowi + lowi + BFRNG; if (lowi <= 0) { // cull very first page for (int i = 0, p = 3, sqr = 9; sqr < nxt; i++, p += 2, sqr = p * p) if ((cbuf[i >> 5] & (1 << (i & 31))) == 0) for (int j = (sqr - 3) >> 1; j < BFBTS; j += p) cbuf[j >> 5] |= 1u << j; } else { // cull for the rest of the pages Array.Clear(cbuf, 0, cbuf.Length); if (bpa.Count == 0) { // inite secondar base primes stream bps = nmrtr(); bps.MoveNext(); bps.MoveNext(); bpa.Add((uint)bps.Current); bps.MoveNext(); } // add 3 to base primes array // make sure bpa contains enough base primes... for (PrimeT p = bpa[bpa.Count - 1], sqr = p * p; sqr < nxt; ) { p = bps.Current; bps.MoveNext(); sqr = p * p; bpa.Add((uint)p); } for (int i = 0, lmt = bpa.Count - 1; i < lmt; i++) { var p = (PrimeT)bpa[i]; var s = (p * p - 3) >> 1; // adjust start index based on page lower limit... if (s >= lowi) s -= lowi; else { var r = (lowi - s) % p; s = (r != 0) ? p - r : 0; } for (var j = (uint)s; j < BFBTS; j += p) cbuf[j >> 5] |= 1u << ((int)j); } } } while (bi < BFBTS && (cbuf[bi >> 5] & (1 << (bi & 31))) != 0) ++bi; if (bi < BFBTS) yield return 3 + (((PrimeT)bi + lowi) << 1); else break; // outer loop for next page segment... } } } public IEnumerator<PrimeT> GetEnumerator() { return nmrtr(); } IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); } }</lang>
The above code is about 25 times faster than the Dictionary version at computing the first about 50 million primes (up to a range of one billion), with the actual enumeration of the result sequence now taking longer than the time it takes to cull the composite number representation bits from the arrays, meaning that it is over 50 times faster at actually sieving the primes. The code owes its speed as compared to a naive "one huge memory array" algorithm to using an array size that is the size of the CPU L1 or L2 caches and using bit-packing to fit more number representations into this limited capacity; in this way RAM memory access times are reduced by a factor of from about four to about 10 (depending on CPU and RAM speed) as compared to those naive implementations, and the minor computational cost of the bit manipulations is compensated by a large factor in total execution time.
The time to enumerate the result primes sequence can be reduced somewhat (about a second) by removing the automatic iterator "yield return" statements and converting them into a "rull-your-own" IEnumerable<PrimeT> implementation, but for page segmentation of odds-only, this iteration of the results will still take longer than the time to actually cull the composite numbers from the page arrays.
In order to make further gains in speed, custom methods must be used to avoid using iterator sequences. If this is done, then further gains can be made by extreme wheel factorization (up to about another about four times gain in speed) and multi-processing (with another gain in speed proportional to the actual independent CPU cores used).
Note that all of these gains in speed are not due to C# other than it compiles to reasonably efficient machine code, but rather to proper use of the Sieve of Eratosthenes algorithm.
All of the above unbounded code can be tested by the following "main" method (replace the name "PrimesXXX" with the name of the class to be tested): <lang csharp> static void Main(string[] args) {
Console.WriteLine(PrimesXXX().ElementAt(1000000 - 1)); // zero based indexing... }</lang>
To produce the following output for all tested versions (although some are considerably faster than others):
- Output:
15485863
Chapel
This solution uses nested iterators to create new wheels at run time: <lang chapel>// yield prime and remove all multiples of it from children sieves iter sieve(prime):int {
yield prime;
var last = prime; label candidates for candidate in sieve(prime+1) do { for composite in last..candidate by prime do {
// candidate is a multiple of this prime if composite == candidate then { // remember size of last composite last = composite; // and try the next candidate continue candidates; } }
// candidate cannot need to be removed by this sieve // yield to parent sieve for checking yield candidate; }
}</lang>The topmost sieve needs to be started with 2 (the smallest prime): <lang chapel>config const N = 30; for p in sieve(2) {
write(" ", p); if p > N then { writeln(); break; }
}</lang>
Clojure
Calculates primes up to and including n using a mutable boolean array but otherwise entirely functional code. <lang clojure> (defn primes-to
"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes" [n] (let [root (-> n Math/sqrt long), cmpsts (boolean-array (inc n)), cullp (fn [p] (loop [i (* p p)] (if (<= i n) (do (aset cmpsts i true) (recur (+ i p))))))] (do (dorun (map #(cullp %) (filter #(not (aget cmpsts %)) (range 2 (inc root))))) (filter #(not (aget cmpsts %)) (range 2 (inc n))))))
</lang>
Alternative implementation using Clojure's side-effect oriented list comprehension.
<lang clojure> (defn primes-to
"Returns a lazy sequence of prime numbers less than lim" [lim] (let [refs (boolean-array (+ lim 1) true) root (int (Math/sqrt lim))] (do (doseq [i (range 2 lim) :while (<= i root) :when (aget refs i)] (doseq [j (range (* i i) lim i)] (aset refs j false))) (filter #(aget refs %) (range 2 lim)))))
</lang>
Alternative very slow entirely functional implementation using lazy sequences <lang clojure> (defn primes-to
"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes" [n] (letfn [(nxtprm [cs] ; current candidates (let [p (first cs)] (if (> p (Math/sqrt n)) cs (cons p (lazy-seq (nxtprm (-> (range (* p p) (inc n) p) set (remove cs) rest)))))))] (nxtprm (range 2 (inc n)))))
</lang>
The reason that the above code is so slow is that it has has a high constant factor overhead due to using a (hash) set to remove the composites from the future composites stream, each prime composite stream removal requires a scan across all remaining composites (compared to using an array or vector where only the culled values are referenced, and due to the slowness of Clojure sequence operations as compared to iterator/sequence operations in other languages.
Version based on immutable Vector's
Here is an immutable boolean vector based non-lazy sequence version other than for the lazy sequence operations to output the result: <lang clojure> (defn primes-to
"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes" [max-prime] (let [sieve (fn [s n] (if (<= (* n n) max-prime) (recur (if (s n) (reduce #(assoc %1 %2 false) s (range (* n n) (inc max-prime) n)) s) (inc n)) s))] (->> (-> (reduce conj (vector-of :boolean) (map #(= % %) (range (inc max-prime)))) (assoc 0 false) (assoc 1 false) (sieve 2)) (map-indexed #(vector %2 %1)) (filter first) (map second))))
</lang>
The above code is still quite slow due to the cost of the immutable copy-on-modify operations.
Odds only bit packed mutable array based version
The following code implements an odds-only sieve using a mutable bit packed long array, only using a lazy sequence for the output of the resulting primes: <lang clojure> (set! *unchecked-math* true)
(defn primes-to
"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes" [n] (let [root (-> n Math/sqrt long), rootndx (long (/ (- root 3) 2)), ndx (long (/ (- n 3) 2)), cmpsts (long-array (inc (/ ndx 64))), isprm #(zero? (bit-and (aget cmpsts (bit-shift-right % 6)) (bit-shift-left 1 (bit-and % 63)))), cullp (fn [i] (let [p (long (+ i i 3))]
(loop [i (bit-shift-right (- (* p p) 3) 1)] (if (<= i ndx) (do (let [w (bit-shift-right i 6)] (aset cmpsts w (bit-or (aget cmpsts w) (bit-shift-left 1 (bit-and i 63))))) (recur (+ i p))))))),
cull (fn [] (loop [i 0] (if (<= i rootndx) (do (if (isprm i) (cullp i)) (recur (inc i))))))] (letfn [(nxtprm [i] (if (<= i ndx) (cons (+ i i 3) (lazy-seq (nxtprm (loop [i (inc i)] (if (or (> i ndx) (isprm i)) i (recur (inc i)))))))))] (if (< n 2) nil (cons 3 (if (< n 3) nil (do (cull) (lazy-seq (nxtprm 0)))))))))
</lang>
The above code is about as fast as any "one large sieving array" type of program in any computer language with this level of wheel factorization other than the lazy sequence operations are quite slow: it takes about ten times as long to enumerate the results as it does to do the actual sieving work of culling the composites from the sieving buffer array. The slowness of sequence operations is due to nested function calls, but primarily due to the way Clojure implements closures by "boxing" all arguments (and perhaps return values) as objects in the heap space, which then need to be "un-boxed" as primitives as necessary for integer operations. Some of the facilities provided by lazy sequences are not needed for this algorithm, such as the automatic memoization which means that each element of the sequence is calculated only once; it is not necessary for the sequence values to be retraced for this algorithm.
If further levels of wheel factorization were used, the time to enumerate the resulting primes would be an even higher overhead as compared to the actual composite number culling time, would get even higher if page segmentation were used to limit the buffer size to the size of the CPU L1 cache for many times better memory access times, most important in the culling operations, and yet higher again if multi-processing were used to share to page segment processing across CPU cores.
The following code overcomes many of those limitations by using an internal (OPSeq) "deftype" which implements the ISeq interface as well as the Counted interface to provide immediate count returns (based on a pre-computed total), as well as the IReduce interface which can greatly speed come computations based on the primes sequence (eased greatly using facilities provided by Clojure 1.7.0 and up): <lang clojure> (defn primes-tox
"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes" [n] (let [root (-> n Math/sqrt long), rootndx (long (/ (- root 3) 2)), ndx (max (long (/ (- n 3) 2)) 0), lmt (quot ndx 64), cmpsts (long-array (inc lmt)), cullp (fn [i] (let [p (long (+ i i 3))]
(loop [i (bit-shift-right (- (* p p) 3) 1)] (if (<= i ndx) (do (let [w (bit-shift-right i 6)]
(aset cmpsts w (bit-or (aget cmpsts w) (bit-shift-left 1 (bit-and i 63))))) (recur (+ i p))))))), cull (fn [] (do (aset cmpsts lmt (bit-or (aget cmpsts lmt) (bit-shift-left -2 (bit-and ndx 63)))) (loop [i 0] (when (<= i rootndx) (when (zero? (bit-and (aget cmpsts (bit-shift-right i 6)) (bit-shift-left 1 (bit-and i 63)))) (cullp i)) (recur (inc i)))))) numprms (fn [] (let [w (dec (alength cmpsts))] ;; fast results count bit counter (loop [i 0, cnt (bit-shift-left (alength cmpsts) 6)] (if (> i w) cnt (recur (inc i) (- cnt (java.lang.Long/bitCount (aget cmpsts i))))))))] (if (< n 2) nil (cons 2 (if (< n 3) nil (do (cull) (deftype OPSeq [^long i ^longs cmpsa ^long cnt ^long tcnt] ;; for arrays maybe need to embed the array so that it doesn't get garbage collected??? clojure.lang.ISeq (first [_] (if (nil? cmpsa) nil (+ i i 3))) (next [_] (let [ncnt (inc cnt)] (if (>= ncnt tcnt) nil (OPSeq. (loop [j (inc i)] (let [p? (zero? (bit-and (aget cmpsa (bit-shift-right j 6)) (bit-shift-left 1 (bit-and j 63))))] (if p? j (recur (inc j))))) cmpsa ncnt tcnt)))) (more [this] (let [ncnt (inc cnt)] (if (>= ncnt tcnt) (OPSeq. 0 nil tcnt tcnt) (.next this)))) (cons [this o] (clojure.core/cons o this)) (empty [_] (if (= cnt tcnt) nil (OPSeq. 0 nil tcnt tcnt))) (equiv [this o] (if (or (not= (type this) (type o)) (not= cnt (.cnt ^OPSeq o)) (not= tcnt (.tcnt ^OPSeq o)) (not= i (.i ^OPSeq o))) false true)) clojure.lang.Counted (count [_] (- tcnt cnt)) clojure.lang.Seqable (clojure.lang.Seqable/seq [this] (if (= cnt tcnt) nil this)) clojure.lang.IReduce (reduce [_ f v] (let [c (- tcnt cnt)] (if (<= c 0) nil (loop [ci i, n c, rslt v] (if (zero? (bit-and (aget cmpsa (bit-shift-right ci 6)) (bit-shift-left 1 (bit-and ci 63)))) (let [rrslt (f rslt (+ ci ci 3)), rdcd (reduced? rrslt), nrslt (if rdcd @rrslt rrslt)] (if (or (<= n 1) rdcd) nrslt (recur (inc ci) (dec n) nrslt))) (recur (inc ci) n rslt)))))) (reduce [this f] (if (nil? i) (f) (if (= (.count this) 1) (+ i i 3) (.reduce ^clojure.lang.IReduce (.next this) f (+ i i 3))))) clojure.lang.Sequential Object (toString [this] (if (= cnt tcnt) "()" (.toString (seq (map identity this)))))) (->OPSeq 0 cmpsts 0 (numprms))))))))
</lang>
'(time (count (primes-tox 10000000)))' takes about 40 milliseconds (compiled) to produce 664579.
Due to the better efficiency of the custom CIS type, the primes to the above range can be enumerated in about the same 40 milliseconds that it takes to cull and count the sieve buffer array.
Under Clojure 1.7.0, one can use '(time (reduce (fn [] (+ (long sum) (long v))) 0 (primes-tox 2000000)))' to find "142913828922" as the sum of the primes to two million as per Euler Problem 10 in about 40 milliseconds total with about half the time used for sieving the array and half for computing the sum.
To show how sensitive Clojure is to forms of expression of functions, the simple form '(time (reduce + (primes-tox 2000000)))' takes about twice as long even though it is using the same internal routine for most of the calculation due to the function not having the type coercion's.
Before one considers that this code is suitable for larger ranges, it is still lacks the improvements of page segmentation with pages about the size of the CPU L1/L2 caches (produces about a four times speed up), maximal wheel factorization (to make it another about four times faster), and the use of multi-processing (for a further gain of about 4 times for a multi-core desktop CPU such as an Intel i7), will make the sieving/counting code about 50 times faster than this, although there will only be a moderate improvement in the time to enumerate/process the resulting primes. Using these techniques, the number of primes to one billion can be counted in a small fraction of a second.
Unbounded Versions
For some types of problems such as finding the nth prime (rather than the sequence of primes up to m), a prime sieve with no upper bound is a better tool.
The following variations on an incremental Sieve of Eratosthenes are based on or derived from the Richard Bird sieve as described in the Epilogue of Melissa E. O'Neill's definitive paper:
A Clojure version of Richard Bird's Sieve using Lazy Sequences (sieves odds only) <lang clojure> (defn primes-Bird
"Computes the unbounded sequence of primes using a Sieve of Eratosthenes algorithm by Richard Bird." [] (letfn [(mltpls [p] (let [p2 (* 2 p)] (letfn [(nxtmltpl [c] (cons c (lazy-seq (nxtmltpl (+ c p2)))))] (nxtmltpl (* p p))))), (allmtpls [ps] (cons (mltpls (first ps)) (lazy-seq (allmtpls (next ps))))), (union [xs ys] (let [xv (first xs), yv (first ys)] (if (< xv yv) (cons xv (lazy-seq (union (next xs) ys))) (if (< yv xv) (cons yv (lazy-seq (union xs (next ys)))) (cons xv (lazy-seq (union (next xs) (next ys)))))))), (mrgmltpls [mltplss] (cons (first (first mltplss)) (lazy-seq (union (next (first mltplss)) (mrgmltpls (next mltplss)))))), (minusStrtAt [n cmpsts] (loop [n n, cmpsts cmpsts] (if (< n (first cmpsts)) (cons n (lazy-seq (minusStrtAt (+ n 2) cmpsts))) (recur (+ n 2) (next cmpsts)))))] (do (def oddprms (cons 3 (lazy-seq (let [cmpsts (-> oddprms (allmtpls) (mrgmltpls))] (minusStrtAt 5 cmpsts))))) (cons 2 (lazy-seq oddprms)))))
</lang>
The above code is quite slow due to both that the data structure is a linear merging of prime multiples and due to the slowness of the Clojure sequence operations.
A Clojure version of the tree folding sieve using Lazy Sequences
The following code speeds up the above code by merging the linear sequence of sequences as above by pairs into a right-leaning tree structure: <lang clojure> (defn primes-treeFolding
"Computes the unbounded sequence of primes using a Sieve of Eratosthenes algorithm modified from Bird." [] (letfn [(mltpls [p] (let [p2 (* 2 p)] (letfn [(nxtmltpl [c] (cons c (lazy-seq (nxtmltpl (+ c p2)))))] (nxtmltpl (* p p))))), (allmtpls [ps] (cons (mltpls (first ps)) (lazy-seq (allmtpls (next ps))))), (union [xs ys] (let [xv (first xs), yv (first ys)] (if (< xv yv) (cons xv (lazy-seq (union (next xs) ys))) (if (< yv xv) (cons yv (lazy-seq (union xs (next ys)))) (cons xv (lazy-seq (union (next xs) (next ys)))))))), (pairs [mltplss] (let [tl (next mltplss)] (cons (union (first mltplss) (first tl)) (lazy-seq (pairs (next tl)))))), (mrgmltpls [mltplss] (cons (first (first mltplss)) (lazy-seq (union (next (first mltplss)) (mrgmltpls (pairs (next mltplss))))))), (minusStrtAt [n cmpsts] (loop [n n, cmpsts cmpsts] (if (< n (first cmpsts)) (cons n (lazy-seq (minusStrtAt (+ n 2) cmpsts))) (recur (+ n 2) (next cmpsts)))))] (do (def oddprms (cons 3 (lazy-seq (let [cmpsts (-> oddprms (allmtpls) (mrgmltpls))] (minusStrtAt 5 cmpsts))))) (cons 2 (lazy-seq oddprms)))))
</lang>
The above code is still slower than it should be due to the slowness of Clojure's sequence operations.
A Clojure version of the above tree folding sieve using a custom Co Inductive Sequence
The following code uses a custom "deftype" non-memoizing Co Inductive Stream/Sequence (CIS) implementing the ISeq interface to make the sequence operations more efficient and is about four times faster than the above code: <lang clojure> (defn primes-treeFoldingx
"Computes the unbounded sequence of primes using a Sieve of Eratosthenes algorithm modified from Bird." [] (do (deftype CIS [v cont] clojure.lang.ISeq (first [_] v) (next [_] (if (nil? cont) nil (cont))) (more [this] (let [nv (.next this)] (if (nil? nv) (CIS. nil nil) nv))) (cons [this o] (clojure.core/cons o this)) (empty [_] (if (and (nil? v) (nil? cont)) nil (CIS. nil nil))) (equiv [this o] (loop [cis1 this, cis2 o] (if (nil? cis1) (if (nil? cis2) true false) (if (or (not= (type cis1) (type cis2)) (not= (.v cis1) (.v ^CIS cis2)) (and (nil? (.cont cis1)) (not (nil? (.cont ^CIS cis2)))) (and (nil? (.cont ^CIS cis2)) (not (nil? (.cont cis1))))) false (if (nil? (.cont cis1)) true (recur ((.cont cis1)) ((.cont ^CIS cis2)))))))) (count [this] (loop [cis this, cnt 0] (if (or (nil? cis) (nil? (.cont cis))) cnt (recur ((.cont cis)) (inc cnt))))) clojure.lang.Seqable (seq [this] (if (and (nil? v) (nil? cont)) nil this)) clojure.lang.Sequential Object (toString [this] (if (and (nil? v) (nil? cont)) "()" (.toString (seq (map identity this)))))) (letfn [(mltpls [p] (let [p2 (* 2 p)] (letfn [(nxtmltpl [c] (->CIS c (fn [] (nxtmltpl (+ c p2)))))] (nxtmltpl (* p p))))), (allmtpls [^CIS ps] (->CIS (mltpls (.v ps)) (fn [] (allmtpls ((.cont ps)))))), (union [^CIS xs ^CIS ys] (let [xv (.v xs), yv (.v ys)] (if (< xv yv) (->CIS xv (fn [] (union ((.cont xs)) ys))) (if (< yv xv) (->CIS yv (fn [] (union xs ((.cont ys))))) (->CIS xv (fn [] (union (next xs) ((.cont ys))))))))), (pairs [^CIS mltplss] (let [^CIS tl ((.cont mltplss))] (->CIS (union (.v mltplss) (.v tl)) (fn [] (pairs ((.cont tl))))))), (mrgmltpls [^CIS mltplss] (->CIS (.v ^CIS (.v mltplss)) (fn [] (union ((.cont ^CIS (.v mltplss))) (mrgmltpls (pairs ((.cont mltplss)))))))), (minusStrtAt [n ^CIS cmpsts] (loop [n n, cmpsts cmpsts] (if (< n (.v cmpsts)) (->CIS n (fn [] (minusStrtAt (+ n 2) cmpsts))) (recur (+ n 2) ((.cont cmpsts))))))] (do (def oddprms (->CIS 3 (fn [] (let [cmpsts (-> oddprms (allmtpls) (mrgmltpls))] (minusStrtAt 5 cmpsts))))) (->CIS 2 (fn [] oddprms))))))
</lang>
'(time (count (take-while #(<= (long %) 10000000) (primes-treeFoldingx))))' takes about 3.4 seconds for a range of 10 million.
The above code is useful for ranges up to about fifteen million primes, which is about the first million primes; it is comparable in speed to all of the bounded versions except for the fastest bit packed version which can reasonably be used for ranges about 100 times as large.
Incremental Hash Map based unbounded "odds-only" version
The following code is a version of the O'Neill Haskell code but does not use wheel factorization other than for sieving odds only (although it could be easily added) and uses a Hash Map (constant amortized access time) rather than a Priority Queue (log n access time for combined remove-and-insert-anew operations, which are the majority used for this algorithm) with a lazy sequence for output of the resulting primes; the code has the added feature that it uses a secondary base primes sequence generator and only adds prime culling sequences to the composites map when they are necessary, thus saving time and limiting storage to only that required for the map entries for primes up to the square root of the currently sieved number: <lang clojure> (defn primes-hashmap
"Infinite sequence of primes using an incremental Sieve or Eratosthenes with a Hashmap" [] (letfn [(nxtoddprm [c q bsprms cmpsts] (if (>= c q) ;; only ever equal (let [p2 (* (first bsprms) 2), nbps (next bsprms), nbp (first nbps)] (recur (+ c 2) (* nbp nbp) nbps (assoc cmpsts (+ q p2) p2))) (if (contains? cmpsts c) (recur (+ c 2) q bsprms (let [adv (cmpsts c), ncmps (dissoc cmpsts c)] (assoc ncmps (loop [try (+ c adv)] ;; ensure map entry is unique (if (contains? ncmps try) (recur (+ try adv)) try)) adv))) (cons c (lazy-seq (nxtoddprm (+ c 2) q bsprms cmpsts))))))] (do (def baseoddprms (cons 3 (lazy-seq (nxtoddprm 5 9 baseoddprms {})))) (cons 2 (lazy-seq (nxtoddprm 3 9 baseoddprms {}))))))
</lang>
The above code is slower than the best tree folding version due to the added constant factor overhead of computing the hash functions for every hash map operation even though it has computational complexity of (n log log n) rather than the worse (n log n log log n) for the previous incremental tree folding sieve. It is still about 100 times slower than the sieve based on the bit-packed mutable array due to these constant factor hashing overheads.
There is almost no benefit of converting the above code to use a CIS as most of the time is expended in the hash map functions.
Incremental Priority Queue based unbounded "odds-only" version
In order to implement the O'Neill Priority Queue incremental Sieve of Eratosthenes algorithm, one requires an efficient implementation of a Priority Queue, which is not part of standard Clojure. For this purpose, the most suitable Priority Queue is a binary tree heap based MinHeap algorithm. The following code implements a purely functional (using entirely immutable state) MinHeap Priority Queue providing the required functions of (emtpy-pq) initialization, (getMin-pq pq) to examinte the minimum key/value pair in the queue, (insert-pq pq k v) to add entries to the queue, and (replaceMinAs-pq pq k v) to replaace the minimum entry with a key/value pair as given (it is more efficient that if functions were provided to delete and then re-insert entries in the queue; there is therefore no "delete" or other queue functions supplied as the algorithm does not requrie them: <lang clojure> (deftype PQEntry [k, v]
Object (toString [_] (str "<" k "," v ">")))
(deftype PQNode [^PQEntry ntry, lft, rght, lvl]
Object (toString [_] (str "<" lvl ntry " left: " (str lft) " right: " (str rght) ">")))
(defn empty-pq [] nil)
(defn getMin-pq ^PQEntry [pq] (condp instance? pq
PQEntry pq, PQNode (.ntry ^PQNode pq) nil))
(defn insert-pq [opq k v]
(loop [kv (->PQEntry k v), msk 0, pq opq, cont identity] (condp instance? pq PQEntry (if (< k (.k ^PQEntry pq)) (cont (->PQNode kv pq nil 2)) (cont (->PQNode pq kv nil 2))), PQNode (let [^PQNode pqn pq, kvn (.ntry pqn), l (.lft pqn), r (.rght pqn), nlvl (+ (.lvl pqn) 1), nmsk (if (zero? msk) ;; never ever 0 again with the bit or'ed 1 (bit-or (bit-shift-left nlvl (- 64 (long (quot (Math/log (double nlvl)) (Math/log (double 2)))))) 1) (bit-shift-left msk 1))] (if (<= k (.k ^PQEntry kvn)) (if (neg? nmsk) (recur kvn nmsk r (fn [npq] (cont (->PQNode kv l npq nlvl)))) (recur kvn nmsk l (fn [npq] (cont (->PQNode kv npq r nlvl))))) (if (neg? nmsk) (recur kv nmsk r (fn [npq] (cont (->PQNode kvn l npq nlvl)))) (recur kv nmsk l (fn [npq] (cont (->PQNode kvn npq r nlvl))))))), (cont kv))))
(defn replaceMinAs-pq [opq k v]
(let [kv (->PQEntry k v)] (loop [pq opq, cont identity] (if (instance? PQNode pq) (let [^PQNode pqn pq, l (.lft pqn), r (.rght pqn), lvl (.lvl pqn)] (cond (and (instance? PQEntry r) (> k (.k ^PQEntry r))) (cond ;; right not empty so left is never empty (and (instance? PQEntry l) (> k (.k ^PQEntry l))) ;; both qualify; choose least (if (> (.k ^PQEntry l) (.k ^PQEntry r)) (cont (->PQNode r l kv lvl)) (cont (->PQNode l kv r lvl))), (and (instance? PQNode l) (> k (.k ^PQEntry (.ntry ^PQNode l)))) (let [^PQEntry kvl (.ntry ^PQNode l)] (if (> (.k kvl) (.k ^PQEntry r)) ;; both qualify; choose least (cont (->PQNode r l kv lvl)) (recur l (fn [npq] (cont (->PQNode kvl npq r lvl)))))), :else (cont (->PQNode r l kv lvl))), ;; only right qualifies; no recursion (and (instance? PQNode r) (> k (.k ^PQEntry (.ntry ^PQNode r)))) (let [^PQEntry kvr (.ntry ^PQNode r)] (if (and (instance? PQNode l) (> k (.k ^PQEntry (.ntry ^PQNode l)))) (let [^PQEntry kvl (.ntry ^PQNode l)] (if (> (.k kvl) (.k kvr)) ;; both qualify; choose least (recur r (fn [npq] (cont (->PQNode kvr l npq lvl)))) (recur l (fn [npq] (cont (->PQNode kvl npq r lvl)))))) (recur r (fn [npq] (cont (->PQNode kvr l npq lvl)))))), ;; only right qualifies :else (cond ;; right is empty, but as this is a node, left is never empty (and (instance? PQEntry l) (> k (.k ^PQEntry l))) (cont (->PQNode l kv r lvl)), (and (instance? PQNode l) (> k (.k ^PQEntry (.ntry ^PQNode l)))) (recur l (fn [npq] (cont (->PQNode (.ntry ^PQNode l) npq r lvl)))), :else (cont (->PQNode kv l r lvl))))) ;; just replace contents, leave same (cont kv))))) ;; if was empty or just an entry, just use current entry
</lang>
Note that the above code is written partially using continuation passing style so as to leave the "recur" calls in tail call position as required for efficient looping in Clojure; for practical sieving ranges, the algorithm could likely use just raw function recursion as recursion depth is unlikely to be used beyond a depth of about ten or so, but raw recursion is said to be less code efficient.
The actual incremental sieve using the Priority Queue is as follows, which code uses the same optimizations of postponing the addition of prime composite streams to the queue until the square root of the currently sieved number is reached and using a secondary base primes stream to generate the primes composite stream markers in the queue as was used for the Hash Map version: <lang clojure> (defn primes-pq
"Infinite sequence of primes using an incremental Sieve or Eratosthenes with a Priority Queue" [] (letfn [(nxtoddprm [c q bsprms cmpsts] (if (>= c q) ;; only ever equal (let [p2 (* (first bsprms) 2), nbps (next bsprms), nbp (first nbps)] (recur (+ c 2) (* nbp nbp) nbps (insert-pq cmpsts (+ q p2) p2))) (let [mn (getMin-pq cmpsts)] (if (and mn (>= c (.k mn))) ;; never greater than (recur (+ c 2) q bsprms (loop [adv (.v mn), cmps cmpsts] ;; advance repeat composites for value (let [ncmps (replaceMinAs-pq cmps (+ c adv) adv), nmn (getMin-pq ncmps)] (if (and nmn (>= c (.k nmn))) (recur (.v nmn) ncmps) ncmps)))) (cons c (lazy-seq (nxtoddprm (+ c 2) q bsprms cmpsts)))))))] (do (def baseoddprms (cons 3 (lazy-seq (nxtoddprm 5 9 baseoddprms (empty-pq))))) (cons 2 (lazy-seq (nxtoddprm 3 9 baseoddprms (empty-pq)))))))
</lang>
The above code is faster than the Hash Map version up to about a sieving range of fifteen million or so, but gets progressively slower for larger ranges due to having (n log n log log n) computational complexity rather than the (n log log n) for the Hash Map version, which has a higher constant factor overhead that is overtaken by the extra "log n" factor.
It is slower that the fastest of the tree folding versions (which has the same computational complexity) due to the higher constant factor overhead of the Priority Queue operations (although perhaps a more efficient implementation of the MinHeap Priority Queue could be developed).
Again, these non-mutable array based sieves are about a hundred times slower than even the "one large memory buffer array" version as implemented in the bounded section; a page segmented version of the mutable bit-packed memory array would be several times faster.
All of these algorithms will respond to maximum wheel factorization, getting up to approximately four times faster if this is applied as compared to the the "odds-only" versions.
It is difficult if not impossible to apply efficient multi-processing to the above versions of the unbounded sieves as the next values of the primes sequence are dependent on previous changes of state for the Bird and Tree Folding versions; however, with the addition of a "update the whole Priority Queue (and reheapify)" or "update the Hash Map" to a given page start state functions, it is possible to do for these letter two algorithms; however, even though it is possible and there is some benefit for these latter two implementations, the benefit is less than using mutable arrays due to that the results must be enumerated into a data structure of some sort in order to be passed out of the page function whereas they can be directly enumerated from the array for the mutable array versions.
CMake
<lang cmake>function(eratosthenes var limit)
# Check for integer overflow. With CMake using 32-bit signed integer, # this check fails when limit > 46340. if(NOT limit EQUAL 0) # Avoid division by zero. math(EXPR i "(${limit} * ${limit}) / ${limit}") if(NOT limit EQUAL ${i}) message(FATAL_ERROR "limit is too large, would cause integer overflow") endif() endif()
# Use local variables prime_2, prime_3, ..., prime_${limit} as array. # Initialize array to y => yes it is prime. foreach(i RANGE 2 ${limit}) set(prime_${i} y) endforeach(i)
# Gather a list of prime numbers. set(list) foreach(i RANGE 2 ${limit}) if(prime_${i}) # Append this prime to list. list(APPEND list ${i})
# For each multiple of i, set n => no it is not prime. # Optimization: start at i squared. math(EXPR square "${i} * ${i}") if(NOT square GREATER ${limit}) # Avoid fatal error. foreach(m RANGE ${square} ${limit} ${i}) set(prime_${m} n) endforeach(m) endif() endif(prime_${i}) endforeach(i) set(${var} ${list} PARENT_SCOPE)
endfunction(eratosthenes)</lang>
# Print all prime numbers through 100. eratosthenes(primes 100) message(STATUS "${primes}")
COBOL
<lang cobol>*> Please ignore the asterisks in the first column of the next comments,
- > which are kludges to get syntax highlighting to work.
IDENTIFICATION DIVISION. PROGRAM-ID. Sieve-Of-Eratosthenes.
DATA DIVISION. WORKING-STORAGE SECTION.
01 Max-Number USAGE UNSIGNED-INT. 01 Max-Prime USAGE UNSIGNED-INT.
01 Num-Group. 03 Num-Table PIC X VALUE "P" OCCURS 1 TO 10000000 TIMES DEPENDING ON Max-Number INDEXED BY Num-Index. 88 Is-Prime VALUE "P" FALSE "N". 01 Current-Prime USAGE UNSIGNED-INT.
01 I USAGE UNSIGNED-INT.
PROCEDURE DIVISION. DISPLAY "Enter the limit: " WITH NO ADVANCING ACCEPT Max-Number DIVIDE Max-Number BY 2 GIVING Max-Prime
- *> Set Is-Prime of all non-prime numbers to false.
SET Is-Prime (1) TO FALSE PERFORM UNTIL Max-Prime < Current-Prime
- *> Set current-prime to next prime.
ADD 1 TO Current-Prime PERFORM VARYING Num-Index FROM Current-Prime BY 1 UNTIL Is-Prime (Num-Index) END-PERFORM MOVE Num-Index TO Current-Prime
- *> Set Is-Prime of all multiples of current-prime to
- *> false, starting from current-prime sqaured.
COMPUTE Num-Index = Current-Prime ** 2 PERFORM UNTIL Max-Number < Num-Index SET Is-Prime (Num-Index) TO FALSE SET Num-Index UP BY Current-Prime END-PERFORM END-PERFORM
- *> Display the prime numbers.
PERFORM VARYING Num-Index FROM 1 BY 1 UNTIL Max-Number < Num-Index IF Is-Prime (Num-Index) DISPLAY Num-Index END-IF END-PERFORM
GOBACK .</lang>
Common Lisp
<lang lisp>(defun sieve-of-eratosthenes (maximum)
(let ((sieve (make-array (1+ maximum) :element-type 'bit :initial-element 0))) (loop for candidate from 2 to maximum when (zerop (bit sieve candidate)) collect candidate and do (loop for composite from (expt candidate 2) to maximum by candidate do (setf (bit sieve composite) 1)))))</lang>
Working with odds only (above twice speedup), and only test divide up to the square root of the maximum:
<lang lisp>(defun sieve-odds (maximum) "sieve for odd numbers"
(cons 2 (let ((maxi (ash (1- maximum) -1)) (stop (ash (isqrt maximum) -1))) (let ((sieve (make-array (1+ maxi) :element-type 'bit :initial-element 0))) (loop for i from 1 to maxi when (zerop (sbit sieve i)) collect (1+ (ash i 1)) and when (<= i stop) do (loop for j from (ash (* i (1+ i)) 1) to maxi by (1+ (ash i 1)) do (setf (sbit sieve j) 1)))))))</lang>
While formally a wheel, odds are uniformly spaced and do not require any special processing except for value translation. Wheels proper aren't uniformly spaced and are thus trickier.
D
Simpler Version
Prints all numbers less than the limit.<lang d>import std.stdio, std.algorithm, std.range, std.functional;
uint[] sieve(in uint limit) nothrow @safe {
if (limit < 2) return []; auto composite = new bool[limit];
foreach (immutable uint n; 2 .. cast(uint)(limit ^^ 0.5) + 1) if (!composite[n]) for (uint k = n * n; k < limit; k += n) composite[k] = true;
//return iota(2, limit).filter!(not!composite).array; return iota(2, limit).filter!(i => !composite[i]).array;
}
void main() {
50.sieve.writeln;
}</lang>
- Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
Faster Version
This version uses an array of bits (instead of booleans, that are represented with one byte), and skips even numbers. The output is the same. <lang d>import std.stdio, std.math, std.array;
size_t[] sieve(in size_t m) pure nothrow @safe {
if (m < 3) return null; immutable size_t n = m - 1; enum size_t bpc = size_t.sizeof * 8; auto F = new size_t[((n + 2) / 2) / bpc + 1]; F[] = size_t.max;
size_t isSet(in size_t i) nothrow @safe @nogc { immutable size_t offset = i / bpc; immutable size_t mask = 1 << (i % bpc); return F[offset] & mask; }
void resetBit(in size_t i) nothrow @safe @nogc { immutable size_t offset = i / bpc; immutable size_t mask = 1 << (i % bpc); if ((F[offset] & mask) != 0) F[offset] = F[offset] ^ mask; }
for (size_t i = 3; i <= sqrt(real(n)); i += 2) if (isSet((i - 3) / 2)) for (size_t j = i * i; j <= n; j += 2 * i) resetBit((j - 3) / 2);
Appender!(size_t[]) result; result ~= 2; for (size_t i = 3; i <= n; i += 2) if (isSet((i - 3) / 2)) result ~= i; return result.data;
}
void main() {
50.sieve.writeln;
}</lang>
Extensible Version
(This version is used in the task Extensible prime generator.) <lang d>/// Extensible Sieve of Eratosthenes. struct Prime {
uint[] a = [2];
private void grow() pure nothrow @safe { immutable p0 = a[$ - 1] + 1; auto b = new bool[p0];
foreach (immutable di; a) { immutable uint i0 = p0 / di * di; uint i = (i0 < p0) ? i0 + di - p0 : i0 - p0; for (; i < b.length; i += di) b[i] = true; }
foreach (immutable uint i, immutable bi; b) if (!b[i]) a ~= p0 + i; }
uint opCall(in uint n) pure nothrow @safe { while (n >= a.length) grow; return a[n]; }
}
version (sieve_of_eratosthenes3_main) {
void main() { import std.stdio, std.range, std.algorithm;
Prime prime; uint.max.iota.map!prime.until!q{a > 50}.writeln; }
}</lang>
To see the output (that is the same), compile with -version=sieve_of_eratosthenes3_main
.
Dart
<lang dart>// helper function to pretty print an Iterable String iterableToString(Iterable seq) {
String str = "["; Iterator i = seq.iterator; if (i.moveNext()) str += i.current.toString(); while(i.moveNext()) { str += ", " + i.current.toString(); } return str + "]";
}
main() {
int limit = 1000; int strt = new DateTime.now().millisecondsSinceEpoch; Set<int> sieve = new Set<int>(); for(int i = 2; i <= limit; i++) { sieve.add(i); } for(int i = 2; i * i <= limit; i++) { if(sieve.contains(i)) { for(int j = i * i; j <= limit; j += i) { sieve.remove(j); } } } var sortedValues = new List<int>.from(sieve); int elpsd = new DateTime.now().millisecondsSinceEpoch - strt; print("Found " + sieve.length.toString() + " primes up to " + limit.toString() + " in " + elpsd.toString() + " milliseconds."); print(iterableToString(sortedValues)); // expect sieve.length to be 168 up to 1000...
// Expect.equals(168, sieve.length); }</lang>
- Output:
Found 168 primes up to 1000 in 9 milliseconds. [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, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]
Although it has the characteristics of a true Sieve of Eratosthenes, the above code isn't very efficient due to the remove/modify operations on the Set. Due to these, the computational complexity isn't close to linear with increasing range and it is quite slow for larger sieve ranges compared to compiled languages, taking about four seconds to sieve to ten million.
faster bit-packed array odds-only solution
<lang dart>import 'dart:math';
List<int> SoEOdds(int limit) {
List<int> prms = new List(); if (limit < 2) return prms; prms.add(2); if (limit < 3) return prms; int lmt = (limit - 3) >> 1; int bfsz = (lmt >> 5) + 1; int sqrtlmt = (sqrt(limit) - 3).floor() >> 1; var buf = new List<int>(); for (int i = 0; i < bfsz; i++) buf.add(0); for (int i = 0; i <= sqrtlmt; i++) if ((buf[i >> 5] & (1 << (i & 31))) == 0) { int p = i + i + 3; for (int j = (p * p - 3) >> 1; j <= lmt; j += p) buf[j >> 5] |= 1 << (j & 31); } for (int i = 0; i <= lmt; i++) if ((buf[i >> 5] & (1 << (i & 31))) == 0) prms.add(i + i + 3); return prms;
}
void main() {
int limit = 10000000; int strt = new DateTime.now().millisecondsSinceEpoch; List<int> primes = SoEOdds(limit); int count = primes.length; int elpsd = new DateTime.now().millisecondsSinceEpoch - strt; print("Found " + count.toString() + " primes up to " + limit.toString() + " in " + elpsd.toString() + " milliseconds.");
// print(iterableToString(primes)); // expect sieve.length to be 168 up to 1000... }</lang> The above code is somewhat faster at about ten seconds using the Dart VM to sieve to 100 million, although much faster at about 1.5 seconds run conventionally in Google Chrome using the JavaScript V8 engine, likely due to JavaScript using double floating point numbers for int's whereas the Dart VM uses arbitrary precision integers.
fast page segmented array infinite iterator (sieves odds-only)
<lang dart>import 'dart:collection';
class _SoEPagedIterator implements Iterator<int> {
static const int _BFSZ = 1 << 16; static const int _BFBTS = _BFSZ * 32; static const int _BFRNG = _BFBTS * 2; int _prime = null; int _bi = -1; int _lowi = 0; List<int> _bpa = new List<int>(); Iterator<int> _bps; List<int> _buf = new List<int>(); int get current => this._prime; bool moveNext() { // the following redundant local variable declaration is necessary to // prevent the dart2js compiler from "tree-shaking" and eliminating some // essential code from the below, which doesn't happen with the Dart VM compiler. int lowi = this._lowi; while (true) { if (this._bi < 1) { if (this._bi < 0) { this._bi++; this._prime = 2; break; } int nxt = 3 + (this._lowi << 1) + _BFRNG; this._buf.clear(); for (int i = 0; i < _BFSZ; i++) this._buf.add(0); // faster initialization: if (lowi <= 0) { // special culling for first page as no base primes yet: for (int i = 0, p = 3, sqr = 9; sqr < nxt; i++, p += 2, sqr = p * p) if ((this._buf[i >> 5] & (1 << (i & 31))) == 0) for (int j = (sqr - 3) >> 1; j < _BFBTS; j += p) this._buf[j >> 5] |= 1 << (j & 31); } else { // after the first page: if (this._bpa.length == 0) { // if this is the first page after the zero one: this._bps = new _SoEPagedIterator(); // initialize separate base primes stream: this._bps.moveNext(); // advance to the only even prime of two this._bps.moveNext(); // advance past 2 to the next prime of 3 } // get enough base primes for the page range... for (var lp = this._bps.current, sqr = lp * lp; sqr < nxt; this._bps.moveNext(), lp = this._bps.current, sqr = lp * lp) this._bpa.add(lp); for (var i = 0; i < this._bpa.length; i++) { int p = this._bpa[i]; int s = (p * p - 3) >> 1; if (s >= this._lowi) // adjust start index based on page lower limit... s -= this._lowi; else { int r = (this._lowi - s) % p; s = (r != 0) ? p - r : 0; } for (var j = s; j < _BFBTS; j += p) this._buf[j >> 5] |= 1 << (j & 31); } } } while (this._bi < _BFBTS && ((this._buf[this._bi >> 5] & (1 << (this._bi & 31))) != 0)) this._bi++; // find next marker still with prime status if (this._bi < _BFBTS) { // within buffer: output computed prime this._prime = 3 + ((this._lowi + this._bi++) << 1); break; } else { // beyond buffer range: advance buffer this._bi = 0; this._lowi += _BFBTS; lowi = this._lowi; } } return true; }
}
class SoEPagedOddsInfGen extends IterableBase<int> {
Iterator<int> get iterator { return new _SoEPagedIterator(); }
}
void main() {
int n = 1000000000; int strt = new DateTime.now().millisecondsSinceEpoch; int count = new SoEPagedOddsInfGen().takeWhile((p) => p <= n).length; int elpsd = new DateTime.now().millisecondsSinceEpoch - strt; print("For a range of " + n.toString() + ", " + count.toString() + " primes found in " + elpsd.toString() + " milliseconds.");
}</lang> This version calculates the 50,847,534 primes up to one billion in about 20 seconds under the Dart Virtual Machine (VM). Under the Google Chrome V8 JavaScript engine it should take the same time as the JavaScript from which it was translated of about five seconds, but takes about 14 seconds due to the dart2js compiler adding extra run time array buffer range checks to the innermost culling loops, even though the "check" compiler option was not selected.
Also note the comment at the beginning of the "moveNext()" method about the redundant local variable needed to be added in order for the code to run under JavaScript using Dart 1.5.1 (and possible other versions), which shouldn't happen when it runs fine under the Dart VM without that extra local variable (based only on the private class field _lowi).
Delphi
<lang delphi>program erathostenes;
{$APPTYPE CONSOLE}
type
TSieve = class private fPrimes: TArray<boolean>; procedure InitArray; procedure Sieve; function getNextPrime(aStart: integer): integer; function getPrimeArray: TArray<integer>; public function getPrimes(aMax: integer): TArray<integer>; end;
{ TSieve }
function TSieve.getNextPrime(aStart: integer): integer; begin
result := aStart; while not fPrimes[result] do inc(result);
end;
function TSieve.getPrimeArray: TArray<integer>; var
i, n: integer;
begin
n := 0; setlength(result, length(fPrimes)); // init array with maximum elements for i := 2 to high(fPrimes) do begin if fPrimes[i] then begin result[n] := i; inc(n); end; end; setlength(result, n); // reduce array to actual elements
end;
function TSieve.getPrimes(aMax: integer): TArray<integer>; begin
setlength(fPrimes, aMax); InitArray; Sieve; result := getPrimeArray;
end;
procedure TSieve.InitArray; begin
for i := 2 to high(fPrimes) do fPrimes[i] := true;
end;
procedure TSieve.Sieve; var
i, n, max: integer;
begin
max := length(fPrimes); i := 2; while i < sqrt(max) do begin n := sqr(i); while n < max do begin fPrimes[n] := false; inc(n, i); end; i := getNextPrime(i + 1); end;
end;
var
i: integer; Sieve: TSieve;
begin
Sieve := TSieve.Create; for i in Sieve.getPrimes(100) do write(i, ' '); Sieve.Free; readln;
end.</lang> Output:
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
DWScript
<lang delphi>function Primes(limit : Integer) : array of Integer; var
n, k : Integer; sieve := new Boolean[limit+1];
begin
for n := 2 to Round(Sqrt(limit)) do begin if not sieve[n] then begin for k := n*n to limit step n do sieve[k] := True; end; end; for k:=2 to limit do if not sieve[k] then Result.Add(k);
end;
var r := Primes(50); var i : Integer; for i:=0 to r.High do
PrintLn(r[i]);</lang>
Dylan
With outer to sqrt and inner to p^2 optimizations: <lang dylan>define method primes(n)
let limit = floor(n ^ 0.5) + 1; let sieve = make(limited(<simple-vector>, of: <boolean>), size: n + 1, fill: #t); let last-prime = 2;
while (last-prime < limit) for (x from last-prime ^ 2 to n by last-prime) sieve[x] := #f; end for; block (found-prime) for (n from last-prime + 1 below limit) if (sieve[n] = #f) last-prime := n; found-prime() end; end; last-prime := limit; end block; end while;
for (x from 2 to n) if (sieve[x]) format-out("Prime: %d\n", x); end; end;
end;</lang>
E
E's standard library doesn't have a step-by-N numeric range, so we'll define one, implementing the standard iteration protocol.
def rangeFromBelowBy(start, limit, step) { return def stepper { to iterate(f) { var i := start while (i < limit) { f(null, i) i += step } } } }
The sieve itself:
def eratosthenes(limit :(int > 2), output) { def composite := [].asSet().diverge() for i ? (!composite.contains(i)) in 2..!limit { output(i) composite.addAll(rangeFromBelowBy(i ** 2, limit, i)) } }
Example usage:
? eratosthenes(12, println) # stdout: 2 # 3 # 5 # 7 # 11
eC
Note: this is not a Sieve of Eratosthenes; it is just trial division. <lang cpp> public class FindPrime {
Array<int> primeList { [ 2 ], minAllocSize = 64 }; int index;
index = 3;
bool HasPrimeFactor(int x) { int max = (int)floor(sqrt((double)x)); for(i : primeList) { if(i > max) break; if(x % i == 0) return true; } return false; }
public int GetPrime(int x) { if(x > primeList.count - 1) { for (; primeList.count != x; index += 2) if(!HasPrimeFactor(index)) { if(primeList.count >= primeList.minAllocSize) primeList.minAllocSize *= 2; primeList.Add(index); } } return primeList[x-1]; }
}
class PrimeApp : Application {
FindPrime fp { }; void Main() { int num = argc > 1 ? atoi(argv[1]) : 1; PrintLn(fp.GetPrime(num)); }
} </lang>
Eiffel
<lang eiffel>class
APPLICATION
create
make
feature
make -- Run application. do across primes_through (100) as ic loop print (ic.item.out + " ") end end primes_through (a_limit: INTEGER): LINKED_LIST [INTEGER] -- Prime numbers through `a_limit' require valid_upper_limit: a_limit >= 2 local l_tab: ARRAY [BOOLEAN] do create Result.make create l_tab.make_filled (True, 2, a_limit) across l_tab as ic loop if ic.item then Result.extend (ic.target_index) across ((ic.target_index * ic.target_index) |..| l_tab.upper).new_cursor.with_step (ic.target_index) as id loop l_tab [id.item] := False end end end end
end</lang>
Output:
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
Elixir
<lang elixir>defmodule Prime do
def eratosthenes(limit \\ 1000) do seq = for i <- 2..limit, do: i sieve(seq, []) end defp sieve([], sieved), do: Enum.reverse(sieved) defp sieve([h | t], sieved) do list = for x <- t, rem(x,h)!=0, do: x sieve(list, [h | sieved]) end
end
IO.inspect Prime.eratosthenes(100)</lang>
- Output:
[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]
Emacs Lisp
<lang lisp> (defun sieve-set (limit)
(let ((xs (make-vector (1+ limit) 0))) (loop for i from 2 to limit when (zerop (aref xs i)) collect i and do (loop for m from (* i i) to limit by i do (aset xs m 1)))))
</lang>
Straightforward implementation of sieve of Eratosthenes, 2 times faster:
<lang lisp> (defun sieve (limit)
(let ((xs (vconcat [0 0] (number-sequence 2 limit)))) (loop for i from 2 to (sqrt limit) when (aref xs i) do (loop for m from (* i i) to limit by i do (aset xs m 0))) (remove 0 xs)))
</lang>
Erlang
<lang Erlang> -module( sieve_of_eratosthenes ).
-export( [primes_upto/1] ).
primes_upto( N ) -> Ns = lists:seq( 2, N ), Dict = dict:from_list( [{X, potential_prime} || X <- Ns] ), {Upto_sqrt_ns, _T} = lists:split( erlang:round(math:sqrt(N)), Ns ), {N, Prime_dict} = lists:foldl( fun find_prime/2, {N, Dict}, Upto_sqrt_ns ), lists:sort( dict:fetch_keys(Prime_dict) ).
find_prime( N, {Max, Dict} ) -> find_prime( dict:find(N, Dict), N, {Max, Dict} ).
find_prime( error, _N, Acc ) -> Acc; find_prime( {ok, _Value}, N, {Max, Dict} ) -> {Max, lists:foldl( fun dict:erase/2, Dict, lists:seq(N*N, Max, N) )}. </lang>
- Output:
35> sieve_of_eratosthenes:primes_upto( 20 ). [2,3,5,7,11,13,17,19]
ERRE
<lang ERRE> PROGRAM SIEVE_ORG
! -------------------------------------------------- ! Eratosthenes Sieve Prime Number Program in BASIC ! (da 3 a SIZE*2) from Byte September 1981 !--------------------------------------------------- CONST SIZE%=8190
DIM FLAGS%[SIZE%]
BEGIN
PRINT("Only 1 iteration") COUNT%=0 FOR I%=0 TO SIZE% DO IF FLAGS%[I%]=TRUE THEN !$NULL ELSE PRIME%=I%+I%+3 K%=I%+PRIME% WHILE NOT (K%>SIZE%) DO FLAGS%[K%]=TRUE K%=K%+PRIME% END WHILE PRINT(PRIME%;) COUNT%=COUNT%+1 END IF END FOR PRINT PRINT(COUNT%;" PRIMES")
END PROGRAM </lang>
- Output:
last lines of the output screen
15749 15761 15767 15773 15787 15791 15797 15803 15809 15817 15823 15859 15877 15881 15887 15889 15901 15907 15913 15919 15923 15937 15959 15971 15973 15991 16001 16007 16033 16057 16061 16063 16067 16069 16073 16087 16091 16097 16103 16111 16127 16139 16141 16183 16187 16189 16193 16217 16223 16229 16231 16249 16253 16267 16273 16301 16319 16333 16339 16349 16361 16363 16369 16381 1899 PRIMES
Euphoria
<lang euphoria>constant limit = 1000 sequence flags,primes flags = repeat(1, limit) for i = 2 to 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
primes = {} for i = 2 to limit do
if flags[i] = 1 then primes &= i end if
end for ? primes</lang>
Output:
{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,131,137,139,149,151,157,163,167,173,179, 181,191,193,197,199,211,223,227,229,233,239,241,251,257,263,269,271, 277,281,283,293,307,311,313,317,331,337,347,349,353,359,367,373,379, 383,389,397,401,409,419,421,431,433,439,443,449,457,461,463,467,479, 487,491,499,503,509,521,523,541,547,557,563,569,571,577,587,593,599, 601,607,613,617,619,631,641,643,647,653,659,661,673,677,683,691,701, 709,719,727,733,739,743,751,757,761,769,773,787,797,809,811,821,823, 827,829,839,853,857,859,863,877,881,883,887,907,911,919,929,937,941, 947,953,967,971,977,983,991,997}
F#
Functional
Richard Bird Sieve
This is the idea behind Richard Bird's unbounded code presented in the Epilogue of M. O'Neill's article in Haskell. It is about twice as much code as the Haskell code because F# does not have a built-in lazy list so that the effect must be constructed using a Co-Inductive Stream (CIS) type since no memoization is required, along with the use of recursive functions in combination with sequences. The type inference needs some help with the new CIS type (including selecting the generic type for speed). Note the use of recursive functions to implement multiple non-sharing delayed generating base primes streams, which along with these being non-memoizing means that the entire primes stream is not held in memory as for the original Bird code: <lang fsharp>type CIS<'T> = struct val v:'T val cont:unit->CIS<'T> //'Co Inductive Stream for laziness
new (v,cont) = { v = v; cont = cont } end
type Primes = uint32
let primesBird() =
let rec (^^) (xs: CIS<Prime>) (ys: CIS<Prime>) = // stream merge function let x = xs.v in let y = ys.v if x < y then CIS(x, fun() -> xs.cont() ^^ ys) elif y < x then CIS(y, fun() -> xs ^^ ys.cont()) else CIS(x, fun() -> xs.cont() ^^ ys.cont()) // no duplications let pmltpls p = let rec nxt c = CIS(c, fun() -> nxt (c + p)) in nxt (p * p) let rec allmltps (ps: CIS<Prime>) = CIS(pmltpls ps.v, fun() -> allmltps (ps.cont())) let rec cmpsts (css: CIS<CIS<Prime>>) = CIS(css.v.v, fun() -> (css.v.cont()) ^^ (cmpsts (css.cont()))) let rec minusat n (cs: CIS<Prime>) = if n < cs.v then CIS(n, fun() -> minusat (n + 1u) cs) else minusat (n + 1u) (cs.cont()) let rec baseprms() = CIS(2u, fun() -> minusat 3u (cmpsts (allmltps (baseprms())))) Seq.unfold (fun (ps: CIS<Prime>) -> Some(ps.v, ps.cont())) (minusat 2u (cmpsts (allmltps (baseprms()))))</lang>
The above code sieves all numbers of two and up including all even numbers as per the page specification; the following code makes the very minor changes for an odds-only sieve, with a speedup of over a factor of two: <lang fsharp>type CIS<'T> = struct val v:'T val cont:unit->CIS<'T> //'Co Inductive Stream for laziness
new (v,cont) = { v = v; cont = cont } end
type Prime = uint32
let primesBirdOdds() =
let rec (^^) (xs: CIS<Prime>) (ys: CIS<Prime>) = // stream merge function let x = xs.v in let y = ys.v if x < y then CIS(x, fun() -> xs.cont() ^^ ys) elif y < x then CIS(y, fun() -> xs ^^ ys.cont()) else CIS(x, fun() -> xs.cont() ^^ ys.cont()) // no duplications let pmltpls p = let adv = p + p let rec nxt c = CIS(c, fun() -> nxt (c + adv)) in nxt (p * p) let rec allmltps (ps: CIS<Prime>) = CIS(pmltpls ps.v, fun() -> allmltps (ps.cont())) let rec cmpsts (css: CIS<CIS<Prime>>) = CIS(css.v.v, fun() -> (css.v.cont()) ^^ (cmpsts (css.cont()))) let rec minusat n (cs: CIS<Prime>) = if n < cs.v then CIS(n, fun() -> minusat (n + 2u) cs) else minusat (n + 2u) (cs.cont()) let rec oddprms() = CIS(3u, fun() -> minusat 5u (cmpsts (allmltps (oddprms())))) Seq.unfold (fun (ps: CIS<Prime>) -> Some(ps.v, ps.cont())) (CIS(2u, fun() -> minusat 3u (cmpsts (allmltps (oddprms())))))</lang>
Tree Folding Sieve
The above code is still somewhat inefficient as it operates on a linear right extending structure that deepens linearly with increasing base primes (those up to the square root of the currently sieved number); the following code changes the structure into an infinite binary tree-like folding by combining each pair of prime composite streams before further processing as usual - this decreases the processing by approximately a factor of log n: <lang fsharp>type CIS<'T> = struct val v:'T val cont:unit->CIS<'T> //'Co Inductive Stream for laziness
new (v,cont) = { v = v; cont = cont } end
type Prime = uint32
let primesTreeFold() =
let rec (^^) (xs: CIS<Prime>) (ys: CIS<Prime>) = // merge streams; no duplicates let x = xs.v in let y = ys.v if x < y then CIS(x, fun() -> xs.cont() ^^ ys) elif y < x then CIS(y, fun() -> xs ^^ ys.cont()) else CIS(x, fun() -> xs.cont() ^^ ys.cont()) let pmltpls p = let adv = p + p let rec nxt c = CIS(c, fun() -> nxt (c + adv)) in nxt (p * p) let rec allmltps (ps: CIS<Prime>) = CIS(pmltpls ps.v, fun() -> allmltps (ps.cont())) let rec pairs (css: CIS<CIS<Prime>>) = let ncss = css.cont() CIS(css.v ^^ ncss.v, fun() -> pairs (ncss.cont())) let rec cmpsts (css: CIS<CIS<Prime>>) = CIS(css.v.v, fun() -> (css.v.cont()) ^^ (cmpsts << pairs << css.cont)()) let rec minusat n (cs: CIS<Prime>) = if n < cs.v then CIS(n, fun() -> minusat (n + 2u) cs) else minusat (n + 2u) (cs.cont()) let rec oddprms() = CIS(3u, fun() -> (minusat 5u << cmpsts << allmltps) (oddprms())) Seq.unfold (fun (ps: CIS<Prime>) -> Some(ps.v, ps.cont())) (CIS(2u, fun() -> (minusat 3u << cmpsts << allmltps) (oddprms())))</lang>
The above code is over four times faster than the "BirdOdds" version and is moderately useful for a range of the first million primes or so.
Priority Queue Sieve
In order to investigate Priority Queue Sieves as espoused by O'Neill in the referenced article, one must find an equivalent implementation of a Min Heap Priority Queue as used by her. There is such an purely functional implementation in RosettaCode translated from the Haskell code she used, from which the essential parts are duplicated here (Note that the key value is given an integer type in order to avoid the inefficiency of F# in generic comparison): <lang fsharp>[<RequireQualifiedAccess>] module MinHeap =
type HeapEntry<'V> = struct val k:uint32 val v:'V new(k,v) = {k=k;v=v} end [<CompilationRepresentation(CompilationRepresentationFlags.UseNullAsTrueValue)>] [<NoEquality; NoComparison>] type PQ<'V> = | Mt | Br of HeapEntry<'V> * PQ<'V> * PQ<'V>
let empty = Mt
let peekMin = function | Br(kv, _, _) -> Some(kv.k, kv.v) | _ -> None
let rec push wk wv = function | Mt -> Br(HeapEntry(wk, wv), Mt, Mt) | Br(vkv, ll, rr) -> if wk <= vkv.k then Br(HeapEntry(wk, wv), push vkv.k vkv.v rr, ll) else Br(vkv, push wk wv rr, ll)
let private siftdown wk wv pql pqr = let rec sift pl pr = match pl with | Mt -> Br(HeapEntry(wk, wv), Mt, Mt) | Br(vkvl, pll, plr) -> match pr with | Mt -> if wk <= vkvl.k then Br(HeapEntry(wk, wv), pl, Mt) else Br(vkvl, Br(HeapEntry(wk, wv), Mt, Mt), Mt) | Br(vkvr, prl, prr) -> if wk <= vkvl.k && wk <= vkvr.k then Br(HeapEntry(wk, wv), pl, pr) elif vkvl.k <= vkvr.k then Br(vkvl, sift pll plr, pr) else Br(vkvr, pl, sift prl prr) sift pql pqr
let replaceMin wk wv = function | Mt -> Mt | Br(_, ll, rr) -> siftdown wk wv ll rr</lang>
Except as noted for any individual code, all of the following codes need the following prefix code in order to implement the non-memoizing Co-Inductive Streams (CIS's) and to set the type of particular constants used in the codes to the same time as the "Prime" type: <lang fsharp>type CIS<'T> = struct val v: 'T val cont: unit -> CIS<'T> new(v,cont) = {v=v;cont=cont} end type Prime = uint32 let frstprm = 2u let frstoddprm = 3u let inc1 = 1u let inc = 2u</lang>
The F# equivalent to O'Neill's "odds-only" code is then implemented as follows, which needs the included changed prefix in order to change the primes type to a larger one to prevent overflow (as well the key type for the MinHeap needs to be changed from uint32 to uint64); it is functionally the same as the O'Neill code other than for minor changes to suit the use of CIS streams and the option output of the "peekMin" function: <lang fsharp>type CIS<'T> = struct val v: 'T val cont: unit -> CIS<'T> new(v,cont) = {v=v;cont=cont} end type Prime = uint64 let frstprm = 2UL let frstoddprm = 3UL let inc = 2UL
let primesPQ() =
let pmult p (xs: CIS<Prime>) = // does map (* p) xs let rec nxtm (cs: CIS<Prime>) = CIS(p * cs.v, fun() -> nxtm (cs.cont())) in nxtm xs let insertprime p xs table = MinHeap.push (p * p) (pmult p xs) table let rec sieve' (ns: CIS<Prime>) table = let nextcomposite = match MinHeap.peekMin table with | None -> ns.v // never happens | Some (k, _) -> k let rec adjust table = let (n, advs) = match MinHeap.peekMin table with | None -> (ns.v, ns.cont()) // never happens | Some kv -> kv if n <= ns.v then adjust (MinHeap.replaceMin advs.v (advs.cont()) table) else table if nextcomposite <= ns.v then sieve' (ns.cont()) (adjust table) else let n = ns.v in CIS(n, fun() -> let nxtns = ns.cont() in sieve' nxtns (insertprime n nxtns table)) let rec sieve (ns: CIS<Prime>) = let n = ns.v in CIS(n, fun() -> let nxtns = ns.cont() in sieve' nxtns (insertprime n nxtns MinHeap.empty)) let odds = // is the odds CIS from 3 up let rec nxto i = CIS(i, fun() -> nxto (i + inc)) in nxto frstoddprm Seq.unfold (fun (cis: CIS<Prime>) -> Some(cis.v, cis.cont())) (CIS(frstprm, fun() -> (sieve odds)))</lang>
However, that algorithm suffers in speed and memory use due to over-eager adding of prime composite streams to the queue such that the queue used is much larger than it needs to be and a much larger range of primes number must be used in order to avoid numeric overflow on the square of the prime added to the queue. The following code corrects that by using a secondary (actually a multiple of) base primes streams which are constrained to be based on a prime that is no larger than the square root of the currently sieved number - this permits the use of much smaller Prime types as per the default prefix: <lang fsharp>let primesPQx() =
let rec nxtprm n pq q (bps: CIS<Prime>) = if n >= q then let bp = bps.v in let adv = bp + bp let nbps = bps.cont() in let nbp = nbps.v nxtprm (n + inc) (MinHeap.push (n + adv) adv pq) (nbp * nbp) nbps else let ck, cv = match MinHeap.peekMin pq with | None -> (q, inc) // only happens until first insertion | Some kv -> kv if n >= ck then let rec adjpq ck cv pq = let npq = MinHeap.replaceMin (ck + cv) cv pq match MinHeap.peekMin npq with | None -> npq // never happens | Some(nk, nv) -> if n >= nk then adjpq nk nv npq else npq nxtprm (n + inc) (adjpq ck cv pq) q bps else CIS(n, fun() -> nxtprm (n + inc) pq q bps) let rec oddprms() = CIS(frstoddprm, fun() -> nxtprm (frstoddprm + inc) MinHeap.empty (frstoddprm * frstoddprm) (oddprms())) Seq.unfold (fun (cis: CIS<Prime>) -> Some(cis.v, cis.cont())) (CIS(frstprm, fun() -> (oddprms())))</lang>
The above code is well over five times faster than the previous translated O'Neill version for the given variety of reasons.
Although slightly faster than the Tree Folding code, this latter code is also limited in practical usefulness to about the first one to ten million primes or so.
All of the above codes can be tested in the F# REPL with the following to produce the millionth prime (the "nth" function is zero based):
> primesXXX() |> Seq.nth 999999;;
where primesXXX() is replaced by the given primes generating function to be tested, and which all produce the following output (after a considerable wait in some cases):
- Output:
val it : Prime = 15485863u
Imperative
The following code is written in functional style other than it uses a mutable bit array to sieve the composites:
<lang fsharp>let primes limit =
let buf = System.Collections.BitArray(int limit + 1, true) let cull p = { p * p .. p .. limit } |> Seq.iter (fun c -> buf.[int c] <- false) { 2u .. uint32 (sqrt (double limit)) } |> Seq.iter (fun c -> if buf.[int c] then cull c) { 2u .. limit } |> Seq.map (fun i -> if buf.[int i] then i else 0u) |> Seq.filter ((<>) 0u)
[<EntryPoint>] let main argv =
if argv = null || argv.Length = 0 then failwith "no command line argument for limit!!!" printfn "%A" (primes (System.UInt32.Parse argv.[0]) |> Seq.length) 0 // return an integer exit code</lang>
Substituting the following minor changes to the code for the "primes" function will only deal with the odd prime candidates for a speed up of over a factor of two as well as a reduction of the buffer size by a factor of two:
<lang fsharp>let primes limit =
let lmtb,lmtbsqrt = (limit - 3u) / 2u, (uint32 (sqrt (double limit)) - 3u) / 2u let buf = System.Collections.BitArray(int lmtb + 1, true) let cull i = let p = i + i + 3u in let s = p * (i + 1u) + i in { s .. p .. lmtb } |> Seq.iter (fun c -> buf.[int c] <- false) { 0u .. lmtbsqrt } |> Seq.iter (fun i -> if buf.[int i] then cull i ) let oddprms = { 0u .. lmtb } |> Seq.map (fun i -> if buf.[int i] then i + i + 3u else 0u) |> Seq.filter ((<>) 0u) seq { yield 2u; yield! oddprms }</lang>
The following code uses other functional forms for the inner culling loops of the "primes function" to reduce the use of inefficient sequences so as to reduce the execution time by another factor of almost three:
<lang fsharp>let primes limit =
let lmtb,lmtbsqrt = (limit - 3u) / 2u, (uint32 (sqrt (double limit)) - 3u) / 2u let buf = System.Collections.BitArray(int lmtb + 1, true) let rec culltest i = if i <= lmtbsqrt then let p = i + i + 3u in let s = p * (i + 1u) + i in let rec cullp c = if c <= lmtb then buf.[int c] <- false; cullp (c + p) (if buf.[int i] then cullp s); culltest (i + 1u) in culltest 0u seq {yield 2u; for i = 0u to lmtb do if buf.[int i] then yield i + i + 3u }</lang>
Now much of the remaining execution time is just the time to enumerate the primes as can be seen by turning "primes" into a primes counting function by substituting the following for the last line in the above code doing the enumeration; this makes the code run about a further five times faster:
<lang fsharp> let rec count i acc =
if i > int lmtb then acc else if buf.[i] then count (i + 1) (acc + 1) else count (i + 1) acc count 0 1</lang>
Since the final enumeration of primes is the main remaining bottleneck, it is worth using a "roll-your-own" enumeration implemented as an object expression so as to save many inefficiencies in the use of the built-in seq computational expression by substituting the following code for the last line of the previous codes, which will decrease the execution time by a factor of over three (instead of almost five for the counting-only version, making it almost as fast):
<lang fsharp> let nmrtr() =
let i = ref -2 let rec nxti() = i:=!i + 1;if !i <= int lmtb && not buf.[!i] then nxti() else !i <= int lmtb let inline curr() = if !i < 0 then (if !i= -1 then 2u else failwith "Enumeration not started!!!") else let v = uint32 !i in v + v + 3u { new System.Collections.Generic.IEnumerator<_> with member this.Current = curr() interface System.Collections.IEnumerator with member this.Current = box (curr()) member this.MoveNext() = if !i< -1 then i:=!i+1;true else nxti() member this.Reset() = failwith "IEnumerator.Reset() not implemented!!!"a interface System.IDisposable with member this.Dispose() = () } { new System.Collections.Generic.IEnumerable<_> with member this.GetEnumerator() = nmrtr() interface System.Collections.IEnumerable with member this.GetEnumerator() = nmrtr() :> System.Collections.IEnumerator }</lang>
The various optimization techniques shown here can be used "jointly and severally" on any of the basic versions for various trade-offs between code complexity and performance. Not shown here are other techniques of making the sieve faster, including extending wheel factorization to much larger wheels such as 2/3/5/7, pre-culling the arrays, page segmentation, and multi-processing.
Almost functional Unbounded
the following odds-only implmentations are written in an almost functional style avoiding the use of mutability except for the contents of the data structures uses to hold the state of the and any mutability necessary to implement a "roll-your-own" IEnumberable iterator interface for speed.
Unbounded Dictionary (Mutable Hash Table) Based Sieve
The following code uses the DotNet Dictionary class instead of the above functional Priority Queue to implement the sieve; as average (amortized) hash table access is O(1) rather than O(log n) as for the priority queue, this implementation is slightly faster than the priority queue version for the first million primes and will always be faster for any range above some low range value: <lang fsharp>type Prime = uint32 let frstprm = 2u let frstoddprm = 3u let inc = 2u let primesDict() =
let dct = System.Collections.Generic.Dictionary() let rec nxtprm n q (bps: CIS<Prime>) = if n >= q then let bp = bps.v in let adv = bp + bp let nbps = bps.cont() in let nbp = nbps.v dct.Add(n + adv, adv) nxtprm (n + inc) (nbp * nbp) nbps else if dct.ContainsKey(n) then let adv = dct.[n] dct.Remove(n) |> ignore
// let mutable nn = n + adv // ugly imperative code // while dct.ContainsKey(nn) do nn <- nn + adv // dct.Add(nn, adv)
let rec nxtmt k = // advance to next empty spot if dct.ContainsKey(k) then nxtmt (k + adv) else dct.Add(k, adv) in nxtmt (n + adv) nxtprm (n + inc) q bps else CIS(n, fun() -> nxtprm (n + inc) q bps) let rec oddprms() = CIS(frstoddprm, fun() -> nxtprm (frstoddprm + inc) (frstoddprm * frstoddprm) (oddprms())) Seq.unfold (fun (cis: CIS<Prime>) -> Some(cis.v, cis.cont())) (CIS(frstprm, fun() -> (oddprms())))</lang>
The above code uses functional forms of code (with the imperative style commented out to show how it could be done imperatively) and also uses a recursive non-sharing secondary source of base primes just as for the Priority Queue version. As for the functional codes, the Primes type can easily be changed to "uint64" for wider range of sieving.
In spite of having true O(n log log n) Sieve of Eratosthenes computational complexity where n is the range of numbers to be sieved, the above code is still not particularly fast due to the time required to compute the hash values and manipulations of the hash table.
Unbounded Page Segmented Mutable Array Sieve
All of the above unbounded implementations including the above Dictionary based version are quite slow due to their large constant factor computational overheads, making them more of an intellectual exercise than something practical, especially when larger sieving ranges are required. The following code implements an unbounded page segmented version of the sieve in not that many more lines of code, yet runs about 25 times faster than the Dictionary version for larger ranges of sieving such as to one billion; it uses functional forms without mutability other than for the contents of the arrays and a reference cell used to implement the "roll-your-own" IEnumerable/IEnumerator interfaces for speed: <lang fsharp>let private PGSZBTS = (1 <<< 14) * 8 // sieve buffer size in bits type private PS = class
val i:int val p:uint64 val cmpsts:uint32[] new(i,p,c) = { i=i; p=p; cmpsts=c } end
let rec primesPaged(): System.Collections.Generic.IEnumerable<_> =
let lbpse = lazy (primesPaged().GetEnumerator()) // lazy to prevent race let bpa = ResizeArray() // fills from above sequence as needed let makePg low = let nxt = low + (uint64 PGSZBTS <<< 1) let cmpsts = Array.zeroCreate (PGSZBTS >>> 5) let inline notprm c = cmpsts.[c >>> 5] &&& (1u <<< c) <> 0u let rec nxti c = if c < PGSZBTS && notprm c then nxti (c + 1) else c let inline mrkc c = let w = c >>> 5 cmpsts.[w] <- cmpsts.[w] ||| (1u <<< c) let rec cullf i = if notprm i then cullf (i + 1) else let p = 3 + i + i in let sqr = p * p if uint64 sqr < nxt then let rec cullp c = if c < PGSZBTS then mrkc c; cullp (c + p) else cullf (i + 1) in cullp ((sqr - 3) >>> 1) if low <= 3UL then cullf 0 // special culling for the first page else // cull rest based on a secondary base prime stream let bpse = lbpse.Force() if bpa.Count <= 0 then // move past 2 to 3 bpse.MoveNext() |> ignore; bpse.MoveNext() |> ignore let rec fill np = if np * np >= nxt then let bpasz = bpa.Count let rec cull i = if i < bpasz then let p = bpa.[i] in let sqr = p * p in let pi = int p let strt = if sqr >= low then int (sqr - low) >>> 1 else let r = int (((low - sqr) >>> 1) % p) if r = 0 then 0 else int p - r let rec cullp c = if c < PGSZBTS then mrkc c; cullp (c + pi) cullp strt; cull (i + 1) in cull 0 else bpa.Add(np); bpse.MoveNext() |> ignore fill bpse.Current fill bpse.Current // fill pba as necessary and do cull let ni = nxti 0 in let np = low + uint64 (ni <<< 1) PS(ni, np, cmpsts) let nmrtr() = let ps = ref (PS(0, 0UL, Array.zeroCreate 0)) { new System.Collections.Generic.IEnumerator<_> with member this.Current = (!ps).p interface System.Collections.IEnumerator with member this.Current = box ((!ps).p) member this.MoveNext() = let drps = !ps in let i = drps.i in let p = drps.p let cmpsts = drps.cmpsts in let lmt = cmpsts.Length <<< 5 if p < 3UL then (if p < 2UL then ps := PS(0, 2UL, cmpsts); true else ps := makePg 3UL; true) else let inline notprm c = cmpsts.[c >>> 5] &&& (1u <<< c) <> 0u let rec nxti c = if c < lmt && notprm c then nxti (c + 1) else c let ni = nxti (i + 1) in let np = p + uint64 ((ni - i) <<< 1) if ni < lmt then ps := PS(ni, np, cmpsts); true else ps := makePg np; true member this.Reset() = failwith "IEnumerator.Reset() not implemented!!!" interface System.IDisposable with member this.Dispose() = () } { new System.Collections.Generic.IEnumerable<_> with member this.GetEnumerator() = nmrtr() interface System.Collections.IEnumerable with member this.GetEnumerator() = nmrtr() :> System.Collections.IEnumerator }</lang>
As with all of the efficient unbounded sieves, the above code uses a secondary enumerator of the base primes less than the square root of the currently culled range ("lbpse"), which is this case is a lazy (deffered evaluation) binding so as to avoid a race condition.
The above code is written to output the "uint64" type for very large ranges of primes since there is little computational cost to doing this for this algorithm. As written, the practical range for this sieve is about 16 billion, however, it can be extended to about 10^14 (a week or two of execution time) by setting the "PGSZBTS" constant to the size of the CPU L2 cache rather than the L1 cache (L2 is up to about two Megabytes for modern high end desktop CPU's) at a slight loss of efficiency (a factor of up to two or so) per composite number culling operation due to the slower memory access time.
Even with the custom IEnumerable/IEnumerator interfaces using an object expression (the F# built-in sequence operators are terribly inefficient), the time to enumerate the resulting primes takes longer than the time to actually cull the composite numbers from the sieving arrays. The time to do the actual culling is thus over 50 times faster than done using the Dictionary version. The slowness of enumeration, no matter what further tweaks are done to improve it (each value enumerated will always take function calls and a scan loop that will always take something in the order of 100 CPU clock cycles per value), means that further gains in speed using extreme wheel factorization and multi-processing have little point unless the actual work on the resulting primes is done through use of auxiliary functions not using iteration.
Forth
: prime? ( n -- ? ) here + c@ 0= ; : composite! ( n -- ) here + 1 swap c! ; : sieve ( n -- ) here over erase 2 begin 2dup dup * > while dup prime? if 2dup dup * do i composite! dup +loop then 1+ repeat drop ." Primes: " 2 do i prime? if i . then loop ; 100 sieve
Fortran
<lang fortran>program sieve
implicit none integer, parameter :: i_max = 100 integer :: i logical, dimension (i_max) :: is_prime
is_prime = .true. is_prime (1) = .false. do i = 2, int (sqrt (real (i_max))) if (is_prime (i)) is_prime (i * i : i_max : i) = .false. end do do i = 1, i_max if (is_prime (i)) write (*, '(i0, 1x)', advance = 'no') i end do write (*, *)
end program sieve</lang> Output: <lang>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</lang> Optimised using a pre-computed wheel based on 2: <lang fortran>program sieve_wheel_2
implicit none integer, parameter :: i_max = 100 integer :: i logical, dimension (i_max) :: is_prime
is_prime = .true. is_prime (1) = .false. is_prime (4 : i_max : 2) = .false. do i = 3, int (sqrt (real (i_max))), 2 if (is_prime (i)) is_prime (i * i : i_max : 2 * i) = .false. end do do i = 1, i_max if (is_prime (i)) write (*, '(i0, 1x)', advance = 'no') i end do write (*, *)
end program sieve_wheel_2</lang> Output: <lang>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</lang>
GAP
<lang gap>Eratosthenes := function(n)
local a, i, j; a := ListWithIdenticalEntries(n, true); if n < 2 then return []; else for i in [2 .. n] do if a[i] then j := i*i; if j > n then return Filtered([2 .. n], i -> a[i]); else while j <= n do a[j] := false; j := j + i; od; fi; fi; od; fi;
end;
Eratosthenes(100);
[ 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 ]</lang>
GLBasic
<lang GLBasic>// Sieve of Eratosthenes (find primes) // GLBasic implementation
GLOBAL n%, k%, limit%, flags%[]
limit = 100 // search primes up to this number
DIM flags[limit+1] // GLBasic arrays start at 0
FOR n = 2 TO SQR(limit)
IF flags[n] = 0 FOR k = n*n TO limit STEP n flags[k] = 1 NEXT ENDIF
NEXT
// Display the primes FOR n = 2 TO limit
IF flags[n] = 0 THEN STDOUT n + ", "
NEXT
KEYWAIT </lang>
Go
<lang go>package main
import (
"fmt"
)
func main() {
const limit = 201 // means sieve numbers < 201
// sieve c := make([]bool, limit) // c for composite. false means prime candidate c[1] = true // 1 not considered prime p := 2 for { // first allowed optimization: outer loop only goes to sqrt(limit) p2 := p * p if p2 >= limit { break } // second allowed optimization: inner loop starts at sqr(p) for i := p2; i < limit; i += p { c[i] = true // it's a composite
} // scan to get next prime for outer loop for { p++ if !c[p] { break } } }
// sieve complete. now print a representation. for n := 1; n < limit; n++ { if c[n] { fmt.Print(" .") } else { fmt.Printf("%3d", n) } if n%20 == 0 { fmt.Println("") } }
}</lang> Output:
. 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 . . .131 . . . . .137 .139 . . . . . . . . .149 .151 . . . . .157 . . . . .163 . . .167 . . . . .173 . . . . .179 . 181 . . . . . . . . .191 .193 . . .197 .199 .
A fairly odd sieve tree method: <lang go>package main import "fmt"
type xint uint64 type xgen func()(xint)
func primes() func()(xint) { pp, psq := make([]xint, 0), xint(25)
var sieve func(xint, xint)xgen sieve = func(p, n xint) xgen { m, next := xint(0), xgen(nil) return func()(r xint) { if next == nil { r = n if r <= psq { n += p return }
next = sieve(pp[0] * 2, psq) // chain in pp = pp[1:] psq = pp[0] * pp[0]
m = next() } switch { case n < m: r, n = n, n + p case n > m: r, m = m, next() default: r, n, m = n, n + p, next() } return } }
f := sieve(6, 9) n, p := f(), xint(0)
return func()(xint) { switch { case p < 2: p = 2 case p < 3: p = 3 default: for p += 2; p == n; { p += 2 if p > n { n = f() } } pp = append(pp, p) } return p } }
func main() {
for i, p := 0, primes(); i < 100000; i++ {
fmt.Println(p())
}
}</lang>
See also the concurrent prime sieve example in the "Try Go" window at http://golang.org/
Groovy
This solution uses a BitSet for compactness and speed, but in Groovy, BitSet has full List semantics. It also uses both the "square root of the boundary" shortcut and the "square of the prime" shortcut. <lang groovy>def sievePrimes = { bound ->
def isPrime = new BitSet(bound) isPrime[0..1] = false isPrime[2..bound] = true (2..(Math.sqrt(bound))).each { pc -> if (isPrime[pc]) { ((pc**2)..bound).step(pc) { isPrime[it] = false } } } (0..bound).findAll { isPrime[it] }
}</lang>
Test: <lang groovy>println sievePrimes(100)</lang>
Output:
[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]
GW-BASIC
<lang qbasic>10 INPUT "ENTER NUMBER TO SEARCH TO: ";LIMIT 20 DIM FLAGS(LIMIT) 30 FOR N = 2 TO SQR (LIMIT) 40 IF FLAGS(N) < > 0 GOTO 80 50 FOR K = N * N TO LIMIT STEP N 60 FLAGS(K) = 1 70 NEXT K 80 NEXT N 90 REM DISPLAY THE PRIMES 100 FOR N = 2 TO LIMIT 110 IF FLAGS(N) = 0 THEN PRINT N;", "; 120 NEXT N</lang>
Haskell
Mutable unboxed arrays, odds only
Mutable array of unboxed Bool
s indexed by Int
s, representing odds only:
<lang haskell>import Control.Monad (forM_, when) import Control.Monad.ST import Data.Array.ST import Data.Array.Unboxed
sieveUO :: Int -> UArray Int Bool sieveUO top = runSTUArray $ do
let m = (top-1) `div` 2 r = floor . sqrt $ fromIntegral top + 1 sieve <- newArray (1,m) True -- :: ST s (STUArray s Int Bool) forM_ [1..r `div` 2] $ \i -> do -- prime(i) = 2i+1 isPrime <- readArray sieve i -- ((2i+1)^2-1)`div`2 = 2i(i+1) when isPrime $ do forM_ [2*i*(i+1), 2*i*(i+2)+1..m] $ \j -> do writeArray sieve j False return sieve
primesToUO :: Int -> [Int] primesToUO top | top > 1 = 2 : [2*i + 1 | (i,True) <- assocs $ sieveUO top]
| otherwise = []</lang>
This represents odds only in the array. Empirical orders of growth is ~ in n primes produced, and improving for bigger n s. Memory consumption is low (array seems to be packed) and growing about linearly with n. Can further be significantly sped up by re-writing the forM_
loops with direct recursion, and using unsafeRead
and unsafeWrite
operations.
Immutable arrays
This is the most straightforward: <lang haskell>import Data.Array.Unboxed
primesToA m = sieve 3 (array (3,m) [(i,odd i) | i<-[3..m]] :: UArray Int Bool)
where sieve p a | p*p > m = 2 : [i | (i,True) <- assocs a] | a!p = sieve (p+2) $ a//[(i,False) | i <- [p*p, p*p+2*p..m]] | otherwise = sieve (p+2) a</lang>
Its performance sharply depends on compiler optimizations. Compiled with -O2 flag in the presence of the explicit type signature, it is very fast in producing first few million primes. (//)
is an array update operator.
Immutable arrays, by segments
Works by segments between consecutive primes' squares. Should be the fastest of non-monadic code: <lang haskell>import Data.Array.Unboxed
primesSA = 2 : prs ()
where prs () = 3 : sieve 3 [] (prs ()) sieve x fs (p:ps) = [i*2 + x | (i,True) <- assocs a] ++ sieve (p*p) fs2 ps where q = (p*p-x)`div`2 fs2 = (p,0) : [(s, rem (y-q) s) | (s,y) <- fs] a :: UArray Int Bool a = accumArray (\ b c -> False) True (1,q-1) [(i,()) | (s,y) <- fs, i <- [y+s, y+s+s..q]]</lang>
Basic list-based sieve
Straightforward implementation of the sieve of Eratosthenes in its original bounded form. This finds primes in gaps between the composites, and composites as an enumeration of each prime's multiples. <lang haskell>primesTo m = 2 : eratos [3,5..m] where
eratos (p : xs) | p*p>m = p : xs | otherwise = p : eratos (xs `minus` [p*p, p*p+2*p..m])
minus a@(x:xs) b@(y:ys) = case compare x y of
LT -> x : minus xs b EQ -> minus xs ys GT -> minus a ys
minus a b = a </lang>
Its time complexity is similar to that of optimal trial division because of limitations of Haskell linked lists, where (minus a b)
takes time proportional to length(union a b)
and not (length b)
, as achieved in imperative setting with direct-access memory. Uses ordered list representation of sets.
This is reasonably useful up to ranges of fifteen million or about the first million primes.
Unbounded list based sieve
Unbounded, "naive", too eager to subtract (see above for the definition of minus
):
<lang haskell>primesE = sieve [2..]
where sieve (p:xs) = p : sieve (minus xs [p,p+p..])</lang>
This is slow, with complexity increasing as a square law or worse so that it is only moderately useful up to the first 10,000 primes or so.
The number of active streams can be limited to what's strictly necessary by postponement until after the square of a prime is seen: <lang haskell>primesEQ = 2 : sieve [3..] 4 primesEQ
where sieve (x:xs) q (p:t) | x < q = x : sieve xs q (p:t) | otherwise = sieve (minus xs [q, q+p..]) (head t^2) t</lang>
The above code is useful to a range of the first million primes or so.
The basic gradually-deepening left-leaning (((a-b)-c)- ... )
workflow above can be rearranged into the right-leaning (a-(b+(c+ ... )))
. This is the idea behind Richard Bird's unbounded code presented in M. O'Neill's article.
<lang haskell>primes = _Y $ ((2:) . minus [3..]
. foldr (\x-> (x*x :) . union [x*x+x, x*x+2*x..]) [])
_Y g = g (_Y g) -- = g . g . g . ... non-sharing multistage fixpoint combinator -- = let x = g x in g x -- = g (fix g) two-stage fixpoint combinator -- = let x = g x in x -- = fix g sharing fixpoint combinator
union a@(x:xs) b@(y:ys) = case compare x y of
LT -> x : union xs b EQ -> x : union xs ys GT -> y : union a ys</lang>
Using _Y
is meant to guarantee the separate supply of primes to be independently calculated, recursively, instead of the same one being reused, corecursively; thus the memory footprint is drastically reduced.
The above code is also useful to a range of the first million primes or so. The code can be further optimized by fusing minus [3..]
into one function, preventing a space leak with the newer GHC versions, getting the function gaps
defined below.
List-based tree-merging incremental sieve
Linear merging structure can further be replaced with an indefinitely deepening to the right tree-like structure, (a-(b+((c+d)+( ((e+f)+(g+h)) + ... ))))
.
This merges primes' multiples streams in a tree-like fashion, achieving theoretical time complexity which is only a log n factor above the optimal n log n log (log n), for n primes produced. Indeed, empirically it runs at about ~ n1.2 (for producing first few million primes), similarly to priority-queue–based version of M. O'Neill's, and with very low space complexity too (not counting the produced sequence of course): <lang haskell>primes :: [Int] primes = 2 : _Y ((3 :) . gaps 5 . _U . map(\p-> [p*p, p*p+2*p..]))
gaps k s@(x:xs) | k < x = k : gaps (k+2) s -- ~= ([k,k+2..] \\ s)
| otherwise = gaps (k+2) xs -- when null(s\\[k,k+2..])
_U ((x:xs):t) = x : (union xs . _U . pairs) t -- ~= nub . sort . concat
where pairs (xs:ys:t) = union xs ys : pairs t</lang>
Here's the test entry on Ideone.com, a comparison with more versions, and a similar code with wheel optimization.
Priority Queue based incremental sieve
The above work is derived from the Epilogue of the Melissa E. O'Neill paper which is much referenced with respect to incremental functional sieves; however, that paper is now dated and her comments comparing list based sieves to her original work leading up to a Priority Queue based implementation is no longer current given more recent work such as the above Tree Merging version. Accordingly, a modern "odd's-only" Priority Queue version is developed here for more current comparisons between the above list based incremental sieves and a continuation of O'Neill's work.
In order to implement a Priority Queue version with Haskell, an efficient Priority Queue, which is not part of the standard Haskell libraries is required. A Min Heap implementation is likely best suited for this task in providing the most efficient frequently used peeks of the next item in the queue and replacement of the first item in the queue (not using a "pop" followed by a "push) with "pop" operations then not used at all, and "push" operations used relatively infrequently. Judging by O'Neill's use of an efficient "deleteMinAndInsert" operation which she states "(We provide deleteMinAndInsert becausea heap-based implementation can support this operation with considerably less rearrangement than a deleteMin followed by an insert.)", which statement is true for a Min Heap Priority Queue and not others, and her reference to a priority queue by (Paulson, 1996), the queue she used is likely the one as provided as a simple true functional Min Heap implementation on RosettaCode, from which the essential functions are reproduced here: <lang haskell>data PriorityQ k v = Mt
| Br !k v !(PriorityQ k v) !(PriorityQ k v) deriving (Eq, Ord, Read, Show)
emptyPQ :: PriorityQ k v emptyPQ = Mt
peekMinPQ :: PriorityQ k v -> Maybe (k, v) peekMinPQ Mt = Nothing peekMinPQ (Br k v _ _) = Just (k, v)
pushPQ :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v pushPQ wk wv Mt = Br wk wv Mt Mt pushPQ wk wv (Br vk vv pl pr)
| wk <= vk = Br wk wv (pushPQ vk vv pr) pl | otherwise = Br vk vv (pushPQ wk wv pr) pl
siftdown :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v -> PriorityQ k v siftdown wk wv Mt _ = Br wk wv Mt Mt siftdown wk wv (pl @ (Br vk vv _ _)) Mt
| wk <= vk = Br wk wv pl Mt | otherwise = Br vk vv (Br wk wv Mt Mt) Mt
siftdown wk wv (pl @ (Br vkl vvl pll plr)) (pr @ (Br vkr vvr prl prr))
| wk <= vkl && wk <= vkr = Br wk wv pl pr | vkl <= vkr = Br vkl vvl (siftdown wk wv pll plr) pr | otherwise = Br vkr vvr pl (siftdown wk wv prl prr)
replaceMinPQ :: Ord k => k -> v -> PriorityQ k v -> PriorityQ k v replaceMinPQ wk wv Mt = Mt replaceMinPQ wk wv (Br _ _ pl pr) = siftdown wk wv pl pr </lang>
The "peekMin" function retrieves both of the key and value in a tuple so processing is required to access whichever is required for further processing. As well, the output of the peekMin function is a Maybe with the case of an empty queue providing a Nothing output.
The following code is O'Neill's original odds-only code (without wheel factorization) from her paper slightly adjusted as per the requirements of this Min Heap implementation as laid out above; note the `seq` adjustments to the "adjust" function to make the evaluation of the entry tuple more strict for better efficiency: <lang haskell> -- (c) 2006-2007 Melissa O'Neill. Code may be used freely so long as -- this copyright message is retained and changed versions of the file -- are clearly marked. -- the only changes are the names of the called PQ functions and the -- included processing for the result of the peek function being a maybe tuple.
primesPQ() = 2 : sieve [3,5..]
where sieve [] = [] sieve (x:xs) = x : sieve' xs (insertprime x xs emptyPQ) where insertprime p xs table = pushPQ (p*p) (map (* p) xs) table sieve' [] table = [] sieve' (x:xs) table | nextComposite <= x = sieve' xs (adjust table) | otherwise = x : sieve' xs (insertprime x xs table) where nextComposite = case peekMinPQ table of Just (c, _) -> c adjust table | n <= x = adjust (replaceMinPQ n' ns table) | otherwise = table where (n, n':ns) = case peekMinPQ table of Just tpl -> tpl</lang>
The above code is almost four times slower than the version of the Tree Merging sieve above for the first million primes although it is about the same speed as the original Richard Bird sieve with the "odds-only" adaptation as above. It is slow and uses a huge amount of memory for primarily one reason: over eagerness in adding prime composite streams to the queue, which are added as the primes are listed rather than when they are required as the output primes stream reaches the square of a given base prime; this over eagerness also means that the processed numbers must have a large range in order to not overflow when squared (as in the default Integer = infinite precision integers as used here and by O'Neill, but Int64's or Word64's would give a practical range) which processing of wide range numbers adds processing and memory requirement overhead. Although O'Neill's code is elegant, it also loses some efficiency due to the extensive use of lazy list processing, not all of which is required even for a wheel factorization implementation.
The following code is adjusted to reduce the amount of lazy list processing and to add a secondary base primes stream (or a succession of streams when the combinator is used) so as to overcome the above problems and reduce memory consumption to only that required for the primes below the square root of the currently sieved number; using this means that 32-bit Int's are sufficient for a reasonable range and memory requirements become relatively negligible: <lang haskell>primesPQx :: () -> [Int] primesPQx() = 2 : _Y ((3 :) . sieve 5 emptyPQ 9) -- initBasePrms
where _Y g = g (_Y g) -- non-sharing multi-stage fixpoint combinator OR
-- initBasePrms = 3 : sieve 5 emptyPQ 9 initBasePrms -- single stage
insertprime p table = let adv = 2 * p in let nv = p * p + adv in nv `seq` pushPQ nv adv table sieve n table q bps@(bp:bps') | n >= q = let nbp = head bps' in sieve (n + 2) (insertprime bp table) (nbp * nbp) bps' | n >= nextComposite = sieve (n + 2) (adjust table) q bps | otherwise = n : sieve (n + 2) table q bps where nextComposite = case peekMinPQ table of Nothing -> q -- at beginning when queue empty Just (c, _) -> c adjust table | c <= n = let nc = c + adv in nc `seq` adjust (replaceMinPQ nc adv table) | otherwise = table where (c, adv) = case peekMinPQ table of Just ct -> ct</lang>
The above code is over five times faster than the previous (O'Neill) Priority Queue code and about half again faster than the Tree Merging code for a range of a million primes, and will always be faster as the Min Heap is slightly more efficient than Tree Merging due to better tree balancing.
All of these codes including the list based ones would enjoy about the same constant factor improvement of up to about four times the speed with the application of maximum wheel factorization.
Page Segmented Sieve using a mutable unboxed array
All of the above unbounded sieves are quite limited in practical sieving range due to the large constant factor overheads in computation, making them mostly just interesting intellectual exercises other than for small ranges of about the first million to ten million primes; the following "odds-only page-segmented version using (bit-packed internally) mutable unboxed arrays is about 50 times faster than the fastest of the above algorithms for ranges of about that and higher, making it practical for the first several hundred million primes: <lang haskell>import Data.Bits import Data.Array.Base import Control.Monad.ST import Data.Array.ST (runSTUArray, STUArray(..))
type PrimeType = Int szPGBTS = (2^14) * 8 :: PrimeType -- CPU L1 cache in bits
primesPaged :: () -> [PrimeType] primesPaged() = 2 : _Y (listPagePrms . pagesFrom 0) where
_Y g = g (_Y g) -- non-sharing multi-stage fixpoint combinator listPagePrms (hdpg @ (UArray lowi _ rng _) : tlpgs) = let loop i = if i >= rng then listPagePrms tlpgs else if unsafeAt hdpg i then loop (i + 1) else let ii = lowi + fromIntegral i in case 3 + ii + ii of p -> p `seq` p : loop (i + 1) in loop 0 makePg lowi bps = runSTUArray $ do let limi = lowi + szPGBTS - 1 let nxt = 3 + limi + limi -- last candidate in range cmpsts <- newArray (lowi, limi) False let pbts = fromIntegral szPGBTS let cull (p:ps) = let sqr = p * p in if sqr > nxt then return cmpsts else let pi = fromIntegral p in let cullp c = if c > pbts then return () else do unsafeWrite cmpsts c True cullp (c + pi) in let a = (sqr - 3) `shiftR` 1 in let s = if a >= lowi then fromIntegral (a - lowi) else let r = fromIntegral ((lowi - a) `rem` p) in if r == 0 then 0 else pi - r in do { cullp s; cull ps} if lowi == 0 then do pg0 <- unsafeFreezeSTUArray cmpsts cull $ listPagePrms [pg0] else cull bps pagesFrom lowi bps = let cf lwi = case makePg lwi bps of pg -> pg `seq` pg : cf (lwi + szPGBTS) in cf lowi</lang>
The above code is currently implemented to use "Int" as the prime type but one can change the "PrimeType" to "Int64" (importing Data.Int) or "Word64" (importing Data.Word) to extend the range to its maximum practical range of above 10^14 or so. Note that for larger ranges that one will want to set the "szPGBTS" to something close to the CPU L2 or even L3 cache size (up to 8 Megabytes = 2^23 for an Intel i7) for a slight cost in speed (about a factor of 1.5) but so that it still computes fairly efficiently as to memory access up to those large ranges. It would be quite easy to modify the above code to make the page array size automatically increase in size with increasing range.
The above code takes only a few tens of milliseconds to compute the first million primes and a few seconds to calculate the first 50 million primes, with over half of those times expended in just enumerating the result lazy list, with even worse times when using 64-bit list processing (especially with 32-bit versions of GHC). A further improvement to reduce the computational cost of repeated list processing across the base pages for every page segment would be to store the required base primes (or base prime gaps) in an array that gets extended in size by factors of two (by copying the old array to the new extended array) as the number of base primes increases; in that way the scans across base primes per page segment would just be array accesses which are much faster than list enumeration.
Unlike many other other unbounded examples, this algorithm has the true Sieve of Eratosthenes computational time complexity of O(n log log n) where n is the sieving range with no extra "log n" factor while having a very low computational time cost per composite number cull of less than ten CPU clock cycles per cull (well under as in under 4 clock cycles for the Intel i7 using a page buffer size of the CPU L1 cache size).
There are other ways to make the algorithm faster including high degrees of wheel factorization, which can reduce the number of composite culling operations by a factor of about four for practical ranges, and multi-processing which can reduce the computation time proportionally to the number of available independent CPU cores, but there is little point to these optimizations as long as the lazy list enumeration is the bottleneck as it is starting to be in the above code. To take advantage of those optimizations, functions need to be provided that can compute the desired results without using list processing.
For ranges above about 10^14 where culling spans begin to exceed even an expanded size page array, other techniques need to be adapted such as the use of a "bucket sieve" which tracks the next page that larger base prime culling sequences will "hit" to avoid redundant (and time expensive) start address calculations for base primes that don't "hit" the current page.
However, even with the above code and its limitations for large sieving ranges, the speeds will never come close to as slow as the other "incremental" sieve algorithms, as the time will never exceed about 100 CPU clock cycles per composite number cull, where the fastest of those other algorithms takes many hundreds of CPU clock cycles per cull.
APL-style
Rolling set subtraction over the rolling element-wise addition on integers. Basic: <lang haskell>zipWith (flip (!!)) [0..]
. scanl1 minus . scanl1 (zipWith (+)) $ repeat [2..]</lang>
A bit optimized: <lang haskell>tail . concatMap fst
. (\(x:xs)-> scanl (\(_,a) b-> span (< head b) $ minus a b) ([],x) xs) . scanl1 (zipWith (+) . tail) $ tails [1..]</lang>
An illustration: <lang haskell>> mapM_ (print . take 15) $ take 10 $ scanl1 (zipWith(+)) $ repeat [2..] [ 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15] [ 4, 6, 8,10,12,14,16,18, 20, 22, 24, 26, 28, 30] [ 6, 9,12,15,18,21,24,27, 30, 33, 36, 39, 42, 45] [ 8,12,16,20,24,28,32,36, 40, 44, 48, 52, 56, 60] [10,15,20,25,30,35,40,45, 50, 55, 60, 65, 70, 75] [12,18,24,30,36,42,48,54, 60, 66, 72, 78, 84, 90] [14,21,28,35,42,49,56,63, 70, 77, 84, 91, 98,105] [16,24,32,40,48,56,64,72, 80, 88, 96,104,112,120] [18,27,36,45,54,63,72,81, 90, 99,108,117,126,135] [20,30,40,50,60,70,80,90,100,110,120,130,140,150]
> mapM_ (print . take 15) $ take 10 $ scanl1 (zipWith(+) . tail) $ tails [1..] [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15] [ 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32] [ 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39, 42, 45, 48, 51] [ 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64, 68, 72] [ 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95] [ 36, 42, 48, 54, 60, 66, 72, 78, 84, 90, 96,102,108,114,120] [ 49, 56, 63, 70, 77, 84, 91, 98,105,112,119,126,133,140,147] [ 64, 72, 80, 88, 96,104,112,120,128,136,144,152,160,168,176] [ 81, 90, 99,108,117,126,135,144,153,162,171,180,189,198,207] [100,110,120,130,140,150,160,170,180,190,200,210,220,230,240]</lang>
HicEst
<lang hicest>REAL :: N=100, sieve(N)
sieve = $ > 1 ! = 0 1 1 1 1 ... DO i = 1, N^0.5
IF( sieve(i) ) THEN DO j = i^2, N, i sieve(j) = 0 ENDDO ENDIF
ENDDO
DO i = 1, N
IF( sieve(i) ) WRITE() i
ENDDO </lang>
Icon and Unicon
<lang Icon> procedure main()
sieve(100) end
procedure sieve(n) local p,i,j p:=list(n, 1) every i:=2 to sqrt(n) & j:= i+i to n by i & p[i] == 1 do p[j] := 0 every write(i:=2 to n & p[i] == 1 & i) end</lang>
Alternatively using sets <lang Icon> procedure main()
sieve(100) end
procedure sieve(n) primes := set() every insert(primes,1 to n) every member(primes,i := 2 to n) do every delete(primes,i + i to n by i) delete(primes,1) every write(!sort(primes))
end</lang>
J
This problem is a classic example of how J can be used to represent mathematical concepts.
J uses x|y (residue) to represent the operation of finding the remainder during integer division of y divided by x
<lang J> 10|13 3</lang>
And x|/y (table) gives us a table with all possibilities from x and all possibilities from y.
<lang J> 2 3 4 |/ 2 3 4 0 1 0 2 0 1 2 3 0</lang>
Meanwhile, |/~y (reflex) copies the right argument and uses it as the left argment.
<lang J> |/~ 0 1 2 3 4 0 1 2 3 4 0 0 0 0 0 0 1 0 1 0 0 1 2 0 1 0 1 2 3 0</lang>
(Bigger examples might make the patterns more obvious but they also take up more space.)
By the way, we can ask J to count out the first N integers for us using i. (integers):
<lang J> i. 5 0 1 2 3 4</lang>
Anyways, the 0s in that last table represent the Sieve of Eratosthenes (in a symbolic or mathematical sense), and we can use = (equal) to find them.
<lang J> 0=|/~ i.5 1 0 0 0 0 1 1 1 1 1 1 0 1 0 1 1 0 0 1 0 1 0 0 0 1</lang>
Now all we need to do is add them up, using / (insert) in its single argument role to insert + between each row of that last table.
<lang J> +/0=|/~ i.5 5 1 2 2 3</lang>
The sieve wants the cases where we have two divisors:
<lang J> 2=+/0=|/~ i.5 0 0 1 1 0</lang>
And we just want to know the positions of the 1s in that list, which we can find using I. (indices):
<lang J> I.2=+/0=|/~ i.5 2 3
I.2=+/0=|/~ i.100
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</lang>
And we might want to express this sentence as a definition of a word that lets us use it with an arbitrary argument. There are a variety of ways of doing this. For example:
<lang J>sieve0=: verb def 'I.2=+/0=|/~ i.y'</lang>
Of course, we can also express this in an even more elaborate fashion. The elaboration makes more efficient use of resources for large arguments, at the expense of less efficient use of resources for small arguments:
<lang J>sieve1=: 3 : 0
m=. <.%:y z=. $0 b=. y{.1 while. m>:j=. 1+b i. 0 do. b=. b+.y$(-j){.1 z=. z,j end. z,1+I.-.b )</lang>
"Wheels" may be implemented as follows:
<lang J>sieve2=: 3 : 0
m=. <.%:y z=. y (>:#]) 2 3 5 7 b=. 1,}.y$+./(*/z)$&>(-z){.&.>1 while. m>:j=. 1+b i. 0 do. b=. b+.y$(-j){.1 z=. z,j end. z,1+I.-.b
)</lang>
The use of 2 3 5 7 as wheels provides a 20% time improvement for n=1000 and 2% for n=1e6 but note that sieve2 is still 25 times slower than i.&.(p:inv) for n=1e6. Then again, the value of the sieve of eratosthenes was not efficiency but simplicity. So perhaps we should ignore resource consumption issues and instead focus on intermediate results for reasonably sized example problems?
<lang J> 0=|/~ i.8 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 0 1 0 0 1 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0 0 0 1</lang>
Columns with two "1" values correspond to prime numbers.
Java
<lang java5>import java.util.LinkedList;
public class Sieve{
public static LinkedList<Integer> sieve(int n){ if(n < 2) return new LinkedList<Integer>(); LinkedList<Integer> primes = new LinkedList<Integer>(); LinkedList<Integer> nums = new LinkedList<Integer>();
for(int i = 2;i <= n;i++){ //unoptimized nums.add(i); }
while(nums.size() > 0){ int nextPrime = nums.remove(); for(int i = nextPrime * nextPrime;i <= n;i += nextPrime){ nums.removeFirstOccurrence(i); } primes.add(nextPrime); } return primes; }
}</lang>
To optimize by testing only odd numbers, replace the loop marked "unoptimized" with these lines: <lang java5>nums.add(2); for(int i = 3;i <= n;i += 2){
nums.add(i);
}</lang>
Version using a BitSet: <lang java5>import java.util.LinkedList; import java.util.BitSet;
public class Sieve{
public static LinkedList<Integer> sieve(int n){ LinkedList<Integer> primes = new LinkedList<Integer>(); BitSet nonPrimes = new BitSet(n+1); for (int p = 2; p <= n ; p = nonPrimes.nextClearBit(p+1)) { for (int i = p * p; i <= n; i += p) nonPrimes.set(i); primes.add(p); } return primes; }
}</lang>
Infinite iterator
An iterator that will generate primes indefinitely (perhaps until it runs out of memory), but very slowly.
<lang java5>import java.util.Iterator; import java.util.PriorityQueue; import java.math.BigInteger;
// generates all prime numbers public class InfiniteSieve implements Iterator<BigInteger> {
private static class NonPrimeSequence implements Comparable<NonPrimeSequence> {
BigInteger currentMultiple; BigInteger prime;
public NonPrimeSequence(BigInteger p) { prime = p; currentMultiple = p.multiply(p); // start at square of prime } @Override public int compareTo(NonPrimeSequence other) { // sorted by value of current multiple return currentMultiple.compareTo(other.currentMultiple); }
}
private BigInteger i = BigInteger.valueOf(2); // priority queue of the sequences of non-primes // the priority queue allows us to get the "next" non-prime quickly final PriorityQueue<NonPrimeSequence> nonprimes = new PriorityQueue<NonPrimeSequence>();
@Override public boolean hasNext() { return true; } @Override public BigInteger next() {
// skip non-prime numbers for ( ; !nonprimes.isEmpty() && i.equals(nonprimes.peek().currentMultiple); i = i.add(BigInteger.ONE)) {
// for each sequence that generates this number, // have it go to the next number (simply add the prime) // and re-position it in the priority queue
while (nonprimes.peek().currentMultiple.equals(i)) { NonPrimeSequence x = nonprimes.poll(); x.currentMultiple = x.currentMultiple.add(x.prime); nonprimes.offer(x); } } // prime
// insert a NonPrimeSequence object into the priority queue
nonprimes.offer(new NonPrimeSequence(i)); BigInteger result = i; i = i.add(BigInteger.ONE); return result;
}
public static void main(String[] args) {
Iterator<BigInteger> sieve = new InfiniteSieve(); for (int i = 0; i < 25; i++) { System.out.println(sieve.next()); }
}
}</lang>
- Output:
2 3 5 7 11 13 17 19
Infinite iterator with a faster algorithm (sieves odds-only)
The adding of each discovered prime's incremental step information to the mapping should be postponed until the candidate number reaches the primes square, as it is useless before that point. This drastically reduces the space complexity from O(n/log(n)) to O(sqrt(n/log(n))), in n primes produced, and also lowers the run time complexity due to the use of the hash table based HashMap, which is much more efficient for large ranges.
<lang java5>import java.util.Iterator; import java.util.HashMap;
// generates all prime numbers up to about 10 ^ 19 if one can wait 1000's of years or so... public class SoEInfHashMap implements Iterator<Long> {
long candidate = 2; Iterator<Long> baseprimes = null; long basep = 3; long basepsqr = 9; // HashMap of the sequences of non-primes // the hash map allows us to get the "next" non-prime reasonably quickly // but further allows re-insertions to take amortized constant time final HashMap<Long,Long> nonprimes = new HashMap<>();
@Override public boolean hasNext() { return true; } @Override public Long next() { // do the initial primes separately to initialize the base primes sequence if (this.candidate <= 5L) if (this.candidate++ == 2L) return 2L; else { this.candidate++; if (this.candidate == 5L) return 3L; else { this.baseprimes = new SoEInfHashMap(); this.baseprimes.next(); this.baseprimes.next(); // throw away 2 and 3 return 5L; } } // skip non-prime numbers including squares of next base prime for ( ; this.candidate >= this.basepsqr || //equals nextbase squared => not prime nonprimes.containsKey(this.candidate); candidate += 2) { // insert a square root prime sequence into hash map if to limit if (candidate >= basepsqr) { // if square of base prime, always equal long adv = this.basep << 1; nonprimes.put(this.basep * this.basep + adv, adv); this.basep = this.baseprimes.next(); this.basepsqr = this.basep * this.basep; } // else for each sequence that generates this number, // have it go to the next number (simply add the advance) // and re-position it in the hash map at an emply slot else { long adv = nonprimes.remove(this.candidate); long nxt = this.candidate + adv; while (this.nonprimes.containsKey(nxt)) nxt += adv; //unique keys this.nonprimes.put(nxt, adv); } } // prime long tmp = candidate; this.candidate += 2; return tmp; }
public static void main(String[] args) { int n = 100000000; long strt = System.currentTimeMillis(); SoEInfHashMap sieve = new SoEInfHashMap(); int count = 0; while (sieve.next() <= n) count++; long elpsd = System.currentTimeMillis() - strt; System.out.println("Found " + count + " primes up to " + n + " in " + elpsd + " milliseconds."); }
}</lang>
- Output:
Found 5761455 primes up to 100000000 in 4297 milliseconds.
Infinite iterator with a very fast page segmentation algorithm (sieves odds-only)
Although somewhat faster than the previous infinite iterator version, the above code is still over 10 times slower than an infinite iterator based on array paged segmentation as in the following code, where the time to enumerate/iterate over the found primes (common to all the iterators) is now about half of the total execution time:
<lang java5>import java.util.Iterator; import java.util.ArrayList;
// generates all prime numbers up to about 10 ^ 19 if one can wait 100's of years or so... // practical range is about 10^14 in a week or so... public class SoEPagedOdds implements Iterator<Long> {
private final int BFSZ = 1 << 16; private final int BFBTS = BFSZ * 32; private final int BFRNG = BFBTS * 2; private long bi = -1; private long lowi = 0; private final ArrayList<Integer> bpa = new ArrayList<>(); private Iterator<Long> bps; private final int[] buf = new int[BFSZ]; @Override public boolean hasNext() { return true; } @Override public Long next() { if (this.bi < 1) { if (this.bi < 0) { this.bi = 0; return 2L; } //this.bi muxt be 0 long nxt = 3 + (this.lowi << 1) + BFRNG; if (this.lowi <= 0) { // special culling for first page as no base primes yet: for (int i = 0, p = 3, sqr = 9; sqr < nxt; i++, p += 2, sqr = p * p) if ((this.buf[i >>> 5] & (1 << (i & 31))) == 0) for (int j = (sqr - 3) >> 1; j < BFBTS; j += p) this.buf[j >>> 5] |= 1 << (j & 31); } else { // after the first page: for (int i = 0; i < this.buf.length; i++) this.buf[i] = 0; // clear the sieve buffer if (this.bpa.isEmpty()) { // if this is the first page after the zero one: this.bps = new SoEPagedOdds(); // initialize separate base primes stream: this.bps.next(); // advance past the only even prime of two this.bpa.add(this.bps.next().intValue()); // get the next prime (3 in this case) } // get enough base primes for the page range... for (long p = this.bpa.get(this.bpa.size() - 1), sqr = p * p; sqr < nxt; p = this.bps.next(), this.bpa.add((int)p), sqr = p * p) ; for (int i = 0; i < this.bpa.size() - 1; i++) { long p = this.bpa.get(i); long s = (p * p - 3) >>> 1; if (s >= this.lowi) // adjust start index based on page lower limit... s -= this.lowi; else { long r = (this.lowi - s) % p; s = (r != 0) ? p - r : 0; } for (int j = (int)s; j < BFBTS; j += p) this.buf[j >>> 5] |= 1 << (j & 31); } } } while ((this.bi < BFBTS) && ((this.buf[(int)this.bi >>> 5] & (1 << ((int)this.bi & 31))) != 0)) this.bi++; // find next marker still with prime status if (this.bi < BFBTS) // within buffer: output computed prime return 3 + ((this.lowi + this.bi++) << 1); else { // beyond buffer range: advance buffer this.bi = 0; this.lowi += BFBTS; return this.next(); // and recursively loop } }
public static void main(String[] args) { long n = 1000000000; long strt = System.currentTimeMillis(); Iterator<Long> gen = new SoEPagedOdds(); int count = 0; while (gen.next() <= n) count++; long elpsd = System.currentTimeMillis() - strt; System.out.println("Found " + count + " primes up to " + n + " in " + elpsd + " milliseconds."); }
}</lang>
- Output:
Found 50847534 primes up to 1000000000 in 3201 milliseconds.
JavaScript
<lang javascript>function eratosthenes(limit) {
var primes = []; if (limit >= 2) { var sqrtlmt = Math.sqrt(limit) - 2; var nums = new Array(); // start with an empty Array... for (var i = 2; i <= limit; i++) // and nums.push(i); // only initialize the Array once... for (var i = 0; i <= sqrtlmt; i++) { var p = nums[i] if (p) for (var j = p * p - 2; j < nums.length; j += p) nums[j] = 0; } for (var i = 0; i < nums.length; i++) { var p = nums[i]; if (p) primes.push(p); } } return primes;
}
var primes = eratosthenes(100);
if (typeof print == "undefined")
print = (typeof WScript != "undefined") ? WScript.Echo : alert;
print(primes);</lang> outputs:
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
Substituting the following code for the function for an odds-only algorithm using bit packing for the array produces code that is many times faster than the above:
<lang javascript>function eratosthenes(limit) {
var prms = []; if (limit >= 2) prms = [2]; if (limit >= 3) { var sqrtlmt = (Math.sqrt(limit) - 3) >> 1; var lmt = (limit - 3) >> 1; var bfsz = (lmt >> 5) + 1 var buf = []; for (var i = 0; i < bfsz; i++) buf.push(0); for (var i = 0; i <= sqrtlmt; i++) if ((buf[i >> 5] & (1 << (i & 31))) == 0) { var p = i + i + 3; for (var j = (p * p - 3) >> 1; j <= lmt; j += p) buf[j >> 5] |= 1 << (j & 31); } for (var i = 0; i <= lmt; i++) if ((buf[i >> 5] & (1 << (i & 31))) == 0) prms.push(i + i + 3); } return prms;
}</lang>
While the above code is quite fast especially using an efficient JavaScript engine such as Google Chrome's V8, it isn't as elegant as it could be using the features of the new EcmaScript6 specification when it comes out about the end of 2014 and when JavaScript engines including those of browsers implement that standard in that we might choose to implement an incremental algorithm iterators or generators similar to as implemented in Python or F# (yield). Meanwhile, we can emulate some of those features by using a simulation of an iterator class (which is easier than using a call-back function) for an "infinite" generator based on an Object dictionary as in the following odds-only code written as a JavaScript class:
<lang javascript>var SoEIncClass = (function () {
function SoEIncClass() { this.n = 0; } SoEIncClass.prototype.next = function () { this.n += 2; if (this.n < 7) { // initialization of sequence to avoid runaway: if (this.n < 3) { // only even of two: this.n = 1; // odds from here... return 2; } if (this.n < 5) return 3; this.dict = {}; // n must be 5... this.bps = new SoEIncClass(); // new source of base primes this.bps.next(); // advance past the even prime of two... this.p = this.bps.next(); // first odd prime (3 in this case) this.q = this.p * this.p; // set guard return 5; } else { // past initialization: var s = this.dict[this.n]; // may or may not be defined... if (!s) { // not defined: if (this.n < this.q) // haven't reached the guard: return this.n; // found a prime else { // n === q => not prime but at guard, so: var p2 = this.p << 1; // the span odds-only is twice prime this.dict[this.n + p2] = p2; // add next composite of prime to dict this.p = this.bps.next(); this.q = this.p * this.p; // get next base prime guard return this.next(); // not prime so advance... } } else { // is a found composite of previous base prime => not prime delete this.dict[this.n]; // advance to next composite of this prime: var nxt = this.n + s; while (this.dict[nxt]) nxt += s; // find unique empty slot in dict this.dict[nxt] = s; // to put the next composite for this base prime return this.next(); // not prime so advance... } } }; return SoEIncClass;
})();</lang>
The above code can be used to find the nth prime (which would require estimating the required range limit using the previous fixed range code) by using the following code:
<lang javascript>var gen = new SoEIncClass(); for (var i = 1; i < 1000000; i++, gen.next()); var prime = gen.next();
if (typeof print == "undefined")
print = (typeof WScript != "undefined") ? WScript.Echo : alert;
print(prime);</lang>
to produce the following output (about five seconds using Google Chrome's V8 JavaScript engine):
15485863
The above code is considerably slower than the fixed range code due to the multiple method calls and the use of an object as a dictionary, which (using a hash table as its basis for most implementations) will have about a constant O(1) amortized time per operation but has quite a high constant overhead to convert the numeric indices to strings which are then hashed to be used as table keys for the look-up operations as compared to doing this more directly in implementations such as the Python dict with Python's built-in hashing functions for every supported type.
This can be implemented as an "infinite" odds-only generator using page segmentation for a considerable speed-up with the alternate JavaScript class code as follows:
<lang javascript>var SoEPgClass = (function () {
function SoEPgClass() { this.bi = -1; // constructor resets the enumeration to start... } SoEPgClass.prototype.next = function () { if (this.bi < 1) { if (this.bi < 0) { this.bi++; this.lowi = 0; // other initialization done here... this.bpa = []; return 2; } else { // bi must be zero: var nxt = 3 + (this.lowi << 1) + 262144; this.buf = new Array(); for (var i = 0; i < 4096; i++) // faster initialization: this.buf.push(0); if (this.lowi <= 0) { // special culling for first page as no base primes yet: for (var i = 0, p = 3, sqr = 9; sqr < nxt; i++, p += 2, sqr = p * p) if ((this.buf[i >> 5] & (1 << (i & 31))) === 0) for (var j = (sqr - 3) >> 1; j < 131072; j += p) this.buf[j >> 5] |= 1 << (j & 31); } else { // after the first page: if (!this.bpa.length) { // if this is the first page after the zero one: this.bps = new SoEPgClass(); // initialize separate base primes stream: this.bps.next(); // advance past the only even prime of two this.bpa.push(this.bps.next()); // get the next prime (3 in this case) } // get enough base primes for the page range... for (var p = this.bpa[this.bpa.length - 1], sqr = p * p; sqr < nxt; p = this.bps.next(), this.bpa.push(p), sqr = p * p) ; for (var i = 0; i < this.bpa.length; i++) { var p = this.bpa[i]; var s = (p * p - 3) >> 1; if (s >= this.lowi) // adjust start index based on page lower limit... s -= this.lowi; else { var r = (this.lowi - s) % p; s = (r != 0) ? p - r : 0; } for (var j = s; j < 131072; j += p) this.buf[j >> 5] |= 1 << (j & 31); } } } } while (this.bi < 131072 && this.buf[this.bi >> 5] & (1 << (this.bi & 31))) this.bi++; // find next marker still with prime status if (this.bi < 131072) // within buffer: output computed prime return 3 + ((this.lowi + this.bi++) << 1); else { // beyond buffer range: advance buffer this.bi = 0; this.lowi += 131072; return this.next(); // and recursively loop } }; return SoEPgClass;
})();</lang>
The above code is about fifty times faster (about five seconds to calculate 50 million primes to about a billion on the Google Chrome V8 JavaScript engine) than the above dictionary based code.
The speed for both of these "infinite" solutions will also respond to further wheel factorization techniques, especially for the dictionary based version where any added overhead to deal with the factorization wheel will be negligible compared to the dictionary overhead. The dictionary version would likely speed up about a factor of three or a little more with maximum wheel factorization applied; the page segmented version probably won't gain more than a factor of two and perhaps less due to the overheads of array look-up operations.
jq
Short and sweet ...
<lang jq># Denoting the input by $n, which is assumed to be a positive integer,
- eratosthenes/0 produces an array of primes less than or equal to $n:
def eratosthenes:
# erase(i) sets .[i*j] to false for integral j > 1 def erase(i): if .[i] then reduce range(2; (1 + length) / i) as $j (.; .[i * $j] = false) else . end;
(. + 1) as $n | (($n|sqrt) / 2) as $s | [null, null, range(2; $n)] | reduce (2, 1 + (2 * range(1; $s))) as $i (.; erase($i)) | map(select(.));</lang>
Examples: <lang jq>100 | eratosthenes</lang>
- Output:
[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] <lang jq>1e7 | eratosthenes | length</lang>
- Output:
664579
Liberty BASIC
<lang lb> 'Notice that arrays are globally visible to functions.
'The sieve() function uses the flags() array. 'This is a Sieve benchmark adapted from BYTE 1985 ' May, page 286
size = 7000 dim flags(7001) start = time$("ms") print sieve(size); " primes found." print "End of iteration. Elapsed time in milliseconds: "; time$("ms")-start end
function sieve(size) for i = 0 to size if flags(i) = 0 then prime = i + i + 3 k = i + prime while k <= size flags(k) = 1 k = k + prime wend sieve = sieve + 1 end if next i end function</lang>
Logo
to sieve :limit make "a (array :limit 2) ; initialized to empty lists make "p [] for [i 2 :limit] [ if empty? item :i :a [ queue "p :i for [j [:i * :i] :limit :i] [setitem :j :a :i] ] ] output :p end print sieve 100 ; 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
Lua
<lang lua>function erato(n)
if n < 2 then return {} end local t = {0} -- clears '1' local sqrtlmt = math.sqrt(n) for i = 2, n do t[i] = 1 end for i = 2, sqrtlmt do if t[i] ~= 0 then for j = i*i, n, i do t[j] = 0 end end end local primes = {} for i = 2, n do if t[i] ~= 0 then table.insert(primes, i) end end return primes
end</lang>
The following changes the code to odds-only using the same large array-based algorithm: <lang lua>function erato2(n)
if n < 2 then return {} end if n < 3 then return {2} end local t = {} local lmt = (n - 3) / 2 local sqrtlmt = (math.sqrt(n) - 3) / 2 for i = 0, lmt do t[i] = 1 end for i = 0, sqrtlmt do if t[i] ~= 0 then local p = i + i + 3 for j = (p*p - 3) / 2, lmt, p do t[j] = 0 end end end local primes = {2} for i = 0, lmt do if t[i] ~= 0 then table.insert(primes, i + i + 3) end end return primes
end</lang>
The following code implements an odds-only "infinite" generator style using a table as a hash table, including postponing adding base primes to the table:
<lang lua>function newEratoInf()
local _cand = 0; local _lstbp = 3; local _lstsqr = 9 local _composites = {}; local _bps = nil local _self = {} function _self.next() if _cand < 9 then if _cand < 1 then _cand = 1; return 2 elseif _cand >= 7 then --advance aux source base primes to 3... _bps = newEratoInf() _bps.next(); _bps.next() end end _cand = _cand + 2 if _composites[_cand] == nil then -- may be prime if _cand >= _lstsqr then -- if not the next base prime local adv = _lstbp + _lstbp -- if next base prime _composites[_lstbp * _lstbp + adv] = adv -- add cull seq _lstbp = _bps.next(); _lstsqr = _lstbp * _lstbp -- adv next base prime return _self.next() else return _cand end -- is prime else local v = _composites[_cand] _composites[_cand] = nil local nv = _cand + v while _composites[nv] ~= nil do nv = nv + v end _composites[nv] = v return _self.next() end end return _self
end
gen = newEratoInf() count = 0 while gen.next() <= 10000000 do count = count + 1 end -- sieves to 10 million print(count) </lang>
which outputs "664579" in about three seconds. As this code uses much less memory for a given range than the previous ones and retains efficiency better with range, it is likely more appropriate for larger sieve ranges.
Lucid
prime where prime = 2 fby (n whenever isprime(n)); n = 3 fby n+2; isprime(n) = not(divs) asa divs or prime*prime > N where N is current n; divs = N mod prime eq 0; end; end
recursive
sieve( N ) where N = 2 fby N + 1; sieve( i ) = i fby sieve ( i whenever i mod first i ne 0 ) ; end
M4
<lang M4>define(`lim',100)dnl define(`for',
`ifelse($#,0, ``$0, `ifelse(eval($2<=$3),1, `pushdef(`$1',$2)$5`'popdef(`$1')$0(`$1',eval($2+$4),$3,$4,`$5')')')')dnl
for(`j',2,lim,1,
`ifdef(a[j], `', `j for(`k',eval(j*j),lim,j, `define(a[k],1)')')')
</lang>
Output:
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
Mathematica
<lang Mathematica>Eratosthenes[n_] := Module[{numbers = Range[n]},
Do[If[numbersi != 0, Do[numbersi j = 0, {j, 2, n/i}]], {i, 2, Sqrt[n]}]; Select[numbers, # > 1 &]]
Eratosthenes[100]</lang>
Slightly Optimized Version
The below has been improved to not require so many operations per composite number cull for about two thirds the execution time: <lang Mathematica>Eratosthenes[n_] := Module[{numbers = Range[n]},
Do[If[numbersi != 0, Do[numbersj = 0, {j,i i,n,i}]],{i,2,Sqrt[n]}]; Select[numbers, # > 1 &]]
Eratosthenes[100]</lang>
Sieving Odds-Only Version
The below has been further improved to only sieve odd numbers for a further reduction in execution time by a factor of over two: <lang Mathematica>Eratosthenes2[n_] := Module[{numbers = Range[3, n, 2], limit = (n - 1)/2},
Do[c = numbersi; If[c != 0, Do[numbersj = 0, {j,(c c - 1)/2,limit,c}]], {i,1,(Sqrt[n] - 1)/2}]; Prepend[Select[numbers, # > 1 &], 2]]
Eratosthenes2[100]</lang>
MATLAB
Somewhat optimized true Sieve of Eratosthenes
<lang MATLAB> function P = erato(x) % Sieve of Eratosthenes: returns all primes between 2 and x
P = [0 2:x] ; % Create vector with all ints between 2 and x where % position 1 is hard-coded as 0 since 1 is not a prime.
for (n=2:sqrt(x)) % All primes factors lie between 2 and sqrt(x). if P(n) % If the current value is not 0 (i.e. a prime), P((2*n):n:x) = 0 ; % then replace all further multiples of it with 0. end end % At this point P is a vector with only primes and zeroes.
P = P(P ~= 0) ; % Remove all zeroes from P, leaving only the primes.
return</lang>The optimization lies in fewer steps in the for loop, use of MATLAB's built-in array operations and no modulo calculation.
Limitation: your machine has to be able to allocate enough memory for an array of length x.
A more efficient Sieve
A more efficient Sieve avoids creating a large double precision vector P, instead using a logical array (which consumes 1/8 the memory of a double array of the same size) and only converting to double those values corresponding to primes.
<lang MATLAB> function P = sieveOfEratosthenes(x) ISP = [false true(1, x-1)]; % 1 is not prime, but we start off assuming all numbers between 2 and x are for n = 2:sqrt(x)
if ISP(n) ISP((2*n):n:x) = false; % Multiples of n that are greater than n are not primes end
end % The ISP vector that we have calculated is essentially the output of the ISPRIME function called on 1:x P = find(ISP); % Convert the ISPRIME output to the values of the primes by finding the locations
% of the TRUE values in S.
</lang>
You can compare the output of this function against the PRIMES function included in MATLAB, which performs a somewhat more memory-efficient Sieve (by not storing even numbers, at the expense of a more complicated indexing expression inside the IF statement.)
Maxima
<lang maxima>sieve(n):=block(
[a:makelist(true,n),i:1,j], a[1]:false, do ( i:i+1, unless a[i] do i:i+1, if i*i>n then return(sublist_indices(a,identity)), for j from i*i step i while j<=n do a[j]:false )
)$</lang>
MAXScript
fn eratosthenes n = ( multiples = #() print 2 for i in 3 to n do ( if (findItem multiples i) == 0 then ( print i for j in (i * i) to n by i do ( append multiples j ) ) ) ) eratosthenes 100
Modula-3
Regular version
<lang modula3>MODULE Prime EXPORTS Main;
IMPORT IO;
CONST LastNum = 1000;
VAR a: ARRAY [2..LastNum] OF BOOLEAN;
BEGIN
FOR i := FIRST(a) TO LAST(a) DO a[i] := TRUE; END;
FOR i := FIRST(a) TO LAST(a) DO IF a[i] THEN IO.PutInt(i); IO.Put(" "); FOR j := FIRST(a) TO LAST(a) DO IF j MOD i = 0 THEN a[j] := FALSE; END; END; END; END; IO.Put("\n");
END Prime.</lang>
Advanced version
This version uses more "advanced" types. <lang modula3>(* From the CM3 examples folder (comments removed). *)
MODULE Sieve EXPORTS Main;
IMPORT IO;
TYPE
Number = [2..1000]; Set = SET OF Number;
VAR
prime: Set := Set {FIRST(Number) .. LAST(Number)};
BEGIN
FOR i := FIRST(Number) TO LAST(Number) DO IF i IN prime THEN IO.PutInt(i); IO.Put(" ");
FOR j := i TO LAST(Number) BY i DO prime := prime - Set{j}; END; END; END; IO.Put("\n");
END Sieve.</lang>
MUMPS
<lang MUMPS>ERATO1(HI)
;performs the Sieve of Erotosethenes up to the number passed in. ;This version sets an array containing the primes SET HI=HI\1 KILL ERATO1 ;Don't make it new - we want it to remain after we quit the function NEW I,J,P FOR I=2:1:(HI**.5)\1 FOR J=I*I:I:HI SET P(J)=1 FOR I=2:1:HI S:'$DATA(P(I)) ERATO1(I)=I KILL I,J,P QUIT</lang>
Example:
USER>SET MAX=100,C=0 DO ERATO1^ROSETTA(MAX) USER>WRITE !,"PRIMES BETWEEN 1 AND ",MAX,! FOR SET I=$ORDER(ERATO1(I)) Q:+I<1 WRITE I,", " PRIMES BETWEEN 1 AND 100 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,
NetRexx
Version 1 (slow)
<lang Rexx>/* NetRexx */
options replace format comments java crossref savelog symbols binary
parse arg loWatermark hiWatermark . if loWatermark = | loWatermark = '.' then loWatermark = 1 if hiWatermark = | hiWatermark = '.' then hiWatermark = 200
do
if \loWatermark.datatype('w') | \hiWatermark.datatype('w') then - signal NumberFormatException('arguments must be whole numbers') if loWatermark > hiWatermark then - signal IllegalArgumentException('the start value must be less than the end value')
seive = sieveOfEratosthenes(hiWatermark) primes = getPrimes(seive, loWatermark, hiWatermark).strip
say 'List of prime numbers from' loWatermark 'to' hiWatermark 'via a "Sieve of Eratosthenes" algorithm:' say ' 'primes.changestr(' ', ',') say ' Count of primes:' primes.words
catch ex = Exception
ex.printStackTrace
end
return
method sieveOfEratosthenes(hn = long) public static binary returns Rexx
sv = Rexx(isTrue) sv[1] = isFalse ix = long jx = long
loop ix = 2 while ix * ix <= hn if sv[ix] then loop jx = ix * ix by ix while jx <= hn sv[jx] = isFalse end jx end ix
return sv
method getPrimes(seive = Rexx, lo = long, hi = long) private constant binary returns Rexx
primes = Rexx() loop p_ = lo to hi if \seive[p_] then iterate p_ primes = primes p_ end p_
return primes
method isTrue public constant binary returns boolean
return 1 == 1
method isFalse public constant binary returns boolean
return \isTrue
</lang>
- Output
List of prime numbers from 1 to 200 via a "Sieve of Eratosthenes" algorithm: 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,131,137,139,149,151,157,163,167,173,179,181,191,193,197,199 Count of primes: 46
Version 2 (significantly, i.e. 10 times faster)
<lang NetRexx>/* NetRexx ************************************************************
- Essential improvements:Use boolean instead of Rexx for sv
- and remove methods isTrue and isFalse
- 24.07.2012 Walter Pachl courtesy Kermit Kiser
- /
options replace format comments java crossref savelog symbols binary
parse arg loWatermark hiWatermark . if loWatermark = | loWatermark = '.' then loWatermark = 1 if hiWatermark = | hiWatermark = '.' then hiWatermark = 200000
startdate=Date Date() do
if \loWatermark.datatype('w') | \hiWatermark.datatype('w') then - signal NumberFormatException('arguments must be whole numbers') if loWatermark > hiWatermark then - signal IllegalArgumentException(- 'the start value must be less than the end value') sieve = sieveOfEratosthenes(hiWatermark) primes = getPrimes(sieve, loWatermark, hiWatermark).strip if hiWatermark = 200 Then do say 'List of prime numbers from' loWatermark 'to' hiWatermark say ' 'primes.changestr(' ', ',') end
catch ex = Exception
ex.printStackTrace
end enddate=Date Date() Numeric Digits 20 say (enddate.getTime-startdate.getTime)/1000 'seconds elapsed' say ' Count of primes:' primes.words
return
method sieveOfEratosthenes(hn = int) -
public static binary returns boolean[] true = boolean 1 false = boolean 0 sv = boolean[hn+1] sv[1] = false
ix = int jx = int
loop ix=2 to hn sv[ix]=true end ix
loop ix = 2 while ix * ix <= hn if sv[ix] then loop jx = ix * ix by ix while jx <= hn sv[jx] = false end jx end ix
return sv
method getPrimes(sieve = boolean[], lo = int, hi = int) -
private constant binary Returns Rexx p_ = int primes = Rexx() loop p_ = lo to hi if \sieve[p_] then iterate p_ primes = primes p_ end p_
return primes</lang>
Nial
primes is sublist [ each (2 = sum eachright (0 = mod) [pass,count]), pass ] rest count
Using it
|primes 10 =2 3 5 7
Nim
Based on one of Python solutions: <lang nim>import math, strutils
var is_prime: seq[Bool] = @[] is_prime.add(False) is_prime.add(False)
iterator iprimes_upto(limit: int): int =
for n in high(is_prime) .. limit+2: is_prime.add(True) for n in 2 .. limit + 1: if is_prime[n]: yield n for i in countup((n *% n), limit+1, n): # start at ``n`` squared try: is_prime[i] = False except EInvalidIndex: break
echo("Primes are:")
for x in iprimes_upto(200):
write(stdout, x, " ")
writeln(stdout,"")</lang>
- Output:
Primes are: 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 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199
Niue
<lang Niue>[ dup 2 < ] '<2 ; [ 1 + 'count ; [ <2 [ , ] when ] count times ] 'fill-stack ;
0 'n ; 0 'v ;
[ .clr 0 'n ; 0 'v ; ] 'reset ; [ len 1 - n - at 'v ; ] 'set-base ; [ n 1 + 'n ; ] 'incr-n ; [ mod 0 = ] 'is-factor ; [ dup * ] 'sqr ;
[ set-base
v sqr 2 at > not [ [ dup v = not swap v is-factor and ] remove-if incr-n run ] when ] 'run ;
[ fill-stack run ] 'sieve ;
( tests )
10 sieve .s ( => 2 3 5 7 9 ) reset newline 30 sieve .s ( => 2 3 5 7 11 13 17 19 23 29 ) </lang>
Oberon-2
<lang oberon2>MODULE Primes;
IMPORT Out, Math;
CONST N = 1000;
VAR a: ARRAY N OF BOOLEAN; i, j, m: INTEGER;
BEGIN
(* Set all elements of a to TRUE. *) FOR i := 1 TO N - 1 DO a[i] := TRUE; END;
(* Compute square root of N and convert back to INTEGER. *) m := ENTIER(Math.Sqrt(N));
FOR i := 2 TO m DO IF a[i] THEN FOR j := 2 TO (N - 1) DIV i DO a[i*j] := FALSE; END; END; END;
(* Print all the elements of a that are TRUE. *) FOR i := 2 TO N - 1 DO IF a[i] THEN Out.Int(i, 4); END; END; Out.Ln;
END Primes.</lang>
OCaml
Imperative
<lang ocaml>let sieve n =
let is_prime = Array.create n true in let limit = truncate(sqrt (float (n - 1))) in for i = 2 to limit do if is_prime.(i) then let j = ref (i*i) in while !j < n do is_prime.(!j) <- false; j := !j + i; done done; is_prime.(0) <- false; is_prime.(1) <- false; is_prime</lang>
<lang ocaml>let primes n =
let primes, _ = let sieve = sieve n in Array.fold_right (fun is_prime (xs, i) -> if is_prime then (i::xs, i-1) else (xs, i-1)) sieve ([], Array.length sieve - 1) in primes</lang>
in the top-level:
# primes 100 ;; - : int list = [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]
Functional
<lang ocaml>(* first define some iterators *)
- let fold_iter f init a b =
let rec aux acc i = if i > b then (acc) else aux (f acc i) (succ i) in aux init a ;;
val fold_iter : ('a -> int -> 'a) -> 'a -> int -> int -> 'a = <fun>
- let fold_step f init a b step =
let rec aux acc i = if i > b then (acc) else aux (f acc i) (i + step) in aux init a ;;
val fold_step : ('a -> int -> 'a) -> 'a -> int -> int -> int -> 'a = <fun>
(* remove a given value from a list *)
- let remove li v =
let rec aux acc = function | hd::tl when hd = v -> (List.rev_append acc tl) | hd::tl -> aux (hd::acc) tl | [] -> li in aux [] li ;;
val remove : 'a list -> 'a -> 'a list = <fun>
(* the main function *)
- let primes n =
let li = (* create a list [from 2; ... until n] *) List.rev(fold_iter (fun acc i -> (i::acc)) [] 2 n) in let limit = truncate(sqrt(float n)) in fold_iter (fun li i -> if List.mem i li (* test if (i) is prime *) then (fold_step remove li (i*i) n i) else li) li 2 (pred limit) ;;
val primes : int -> int list = <fun>
- primes 200 ;;
- : int list = [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; 131; 137; 139; 149; 151; 157; 163; 167; 173; 179; 181; 191; 193; 197; 199]</lang>
Another functional version
This uses zero to denote struck-out numbers. It is slightly inefficient as it strikes-out multiples above p rather than p2
<lang ocaml># let rec strike_nth k n l = match l with
| [] -> [] | h :: t -> if k = 0 then 0 :: strike_nth (n-1) n t else h :: strike_nth (k-1) n t;;
val strike_nth : int -> int -> int list -> int list = <fun>
- let primes n =
let limit = truncate(sqrt(float n)) in let rec range a b = if a > b then [] else a :: range (a+1) b in let rec sieve_primes l = match l with | [] -> [] | 0 :: t -> sieve_primes t | h :: t -> if h > limit then List.filter ((<) 0) l else h :: sieve_primes (strike_nth (h-1) h t) in sieve_primes (range 2 n) ;;
val primes : int -> int list = <fun>
- primes 200;;
- : int list = [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; 131; 137; 139; 149; 151; 157; 163; 167; 173; 179; 181; 191; 193; 197; 199]</lang>
Oforth
<lang Oforth>func: eratosthenes(n) { | i j |
ListBuffer newSize(n) dup add(null) seqFrom(2, n) over addAll 2 n sqrt asInteger for: i [ dup at(i) ifNotNull: [ i sq n i step: j [ dup put(j, null) ] ] ] filter(#notNull)
}</lang>
- Output:
>eratosthenes(100) println [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]
Oz
<lang oz>declare
fun {Sieve N} S = {Array.new 2 N true} M = {Float.toInt {Sqrt {Int.toFloat N}}} in for I in 2..M do
if S.I then for J in I*I..N;I do S.J := false end end
end S end
fun {Primes N} S = {Sieve N} in for I in 2..N collect:C do
if S.I then {C I} end
end end
in
{Show {Primes 30}}</lang>
PARI/GP
<lang parigp>Eratosthenes(lim)={
my(v=Vectorsmall(lim\1,unused,1)); forprime(p=2,sqrt(lim), forstep(i=p^2,lim,p, v[i]=0 ) ); for(i=1,lim,if(v[i],print1(i", ")))
};</lang>
An alternate version:
<lang parigp>Sieve(n)= { v=vector(n,unused,1); for(i=2,sqrt(n),
if(v[i], forstep(j=i^2,n,i,v[j]=0)));
for(i=2,n,if(v[i],print1(i))) };</lang>
Pascal
Note: Some Pascal implementations put quite low limits on the size of a set (e.g. Turbo Pascal doesn't allow more than 256 members). To compile on such an implementation, reduce the constant PrimeLimit accordingly. <lang pascal> program primes(output)
const
PrimeLimit = 1000;
var
primes: set of 1 .. PrimeLimit; n, k: integer; needcomma: boolean;
begin
{ calculate the primes } primes := [2 .. PrimeLimit]; for n := 1 to trunc(sqrt(PrimeLimit)) do begin if n in primes then begin k := n*n; while k < PrimeLimit do begin primes := primes - [k]; k := k + n end end end;
{ output the primes } needcomma := false; for n := 1 to PrimeLimit do if n in primes then begin if needcomma then write(', '); write(n); needcomma := true end
end. </lang>
alternative using wheel
Using growing wheel to fill array for sieving for minimal unmark operations. Sieving only with possible-prime factors. <lang pascal> program prim(output); //Sieve of Erathosthenes with fast elimination of multiples of small primes {$IFNDEF FPC}
{$APPTYPE CONSOLE}
{$ENDIF} const
PrimeLimit = 100*1000*1000;//1;
type
tLimit = 1..PrimeLimit;
var
//always initialized with 0 => false at startup primes: array [tLimit] of boolean;
function BuildWheel: longInt; //fill primfield with no multiples of small primes //returns next sieveprime //speedup ~1/3 var
//wheelprimes = 2,3,5,7,11... ; //wheelsize = product [i= 0..wpno-1]wheelprimes[i] > Uint64 i> 13 wheelprimes :array[0..13] of byte; wheelSize,wpno, pr,pw,i, k: LongWord;
begin
//the mother of all numbers 1 ;-) //the first wheel = generator of numbers //not divisible by the small primes first found primes pr := 1; primes[1]:= true; WheelSize := 1;
wpno := 0; repeat inc(pr); //pw = pr projected in wheel of wheelsize pw := pr; if pw > wheelsize then dec(pw,wheelsize); If Primes[pw] then begin
// writeln(pr:10,pw:10,wheelsize:16);
k := WheelSize+1; //turn the wheel (pr-1)-times for i := 1 to pr-1 do begin inc(k,WheelSize); if k<primeLimit then move(primes[1],primes[k-WheelSize],WheelSize) else begin move(primes[1],primes[k-WheelSize],PrimeLimit-WheelSize*i); break; end; end; dec(k); IF k > primeLimit then k := primeLimit; wheelPrimes[wpno] := pr; primes[pr] := false;
inc(wpno); //the new wheelsize WheelSize := k;
//sieve multiples of the new found prime i:= pr; i := i*i; while i <= k do begin primes[i] := false; inc(i,pr); end; end; until WheelSize >= PrimeLimit;
//re-insert wheel-primes // 1 still stays prime while wpno > 0 do begin dec(wpno); primes[wheelPrimes[wpno]] := true; end; BuildWheel := pr+1;
end;
procedure Sieve; var
sieveprime, fakt : LongWord;
begin //primes[1] = true is needed to stop for sieveprime = 2 // at //Search next smaller possible prime
sieveprime := BuildWheel;
//alternative here
//fillchar(primes,SizeOf(Primes),chr(ord(true)));sieveprime := 2; repeat if primes[sieveprime] then begin //eliminate 'possible prime' multiples of sieveprime //must go downwards //2*2 would unmark 4 -> 4*2 = 8 wouldnt be unmarked fakt := PrimeLimit DIV sieveprime; IF fakt < sieveprime then BREAK; repeat //Unmark primes[sieveprime*fakt] := false; //Search next smaller possible prime repeat dec(fakt); until primes[fakt]; until fakt < sieveprime; end; inc(sieveprime); until false; //remove 1 primes[1] := false;
end;
var
prCnt, i : LongWord;
Begin
Sieve; {count the primes } prCnt := 0; for i:= 1 to PrimeLimit do inc(prCnt,Ord(primes[i])); writeln(prCnt,' primes up to ',PrimeLimit);
end.</lang>
output: ( i3 4330 Haswell 3.5 Ghz fpc 2.6.4 -O3 )
5761455 primes up to 100000000 real 0m0.204s user 0m0.193s sys 0m0.013s
Perl
For highest performance and ease, typically a module would be used, such as Math::Prime::Util, Math::Prime::FastSieve, or Math::Prime::XS.
Classic Sieve
<lang perl>sub sieve {
my $n = shift; my @composite; for my $i (2 .. int(sqrt($n))) { if (!$composite[$i]) { for (my $j = $i*$i; $j <= $n; $j += $i) { $composite[$j] = 1; } } } my @primes; for my $i (2 .. $n) { $composite[$i] || push @primes, $i; } @primes;
}</lang>
Odds only (faster)
<lang perl>sub sieve2 {
my($n) = @_; return @{([],[],[2],[2,3],[2,3])[$n]} if $n <= 4;
my @composite; for (my $t = 3; $t*$t <= $n; $t += 2) { if (!$composite[$t]) { for (my $s = $t*$t; $s <= $n; $s += $t*2) { $composite[$s]++ } } } my @primes = (2); for (my $t = 3; $t <= $n; $t += 2) { $composite[$t] || push @primes, $t; } @primes;
}</lang>
Odds only, using vectors for lower memory use
<lang perl>sub dj_vector {
my($end) = @_; return @{([],[],[2],[2,3],[2,3])[$end]} if $end <= 4; $end-- if ($end & 1) == 0; # Ensure end is odd
my ($sieve, $n, $limit, $s_end) = ( , 3, int(sqrt($end)), $end >> 1 ); while ( $n <= $limit ) { for (my $s = ($n*$n) >> 1; $s <= $s_end; $s += $n) { vec($sieve, $s, 1) = 1; } do { $n += 2 } while vec($sieve, $n >> 1, 1) != 0; } my @primes = (2); do { push @primes, 2*$_+1 if !vec($sieve,$_,1) } for (1..int(($end-1)/2)); @primes;
}</lang>
Odds only, using strings for best performance
Compared to array versions, about 2x faster (with 5.16.0 or later) and lower memory. Much faster than the experimental versions below. It's possible a mod-6 or mod-30 wheel could give more improvement, though possibly with obfuscation. The best next step for performance and functionality would be segmenting. <lang perl>sub string_sieve {
my ($n, $i, $s, $d, @primes) = (shift, 7);
local $_ = '110010101110101110101110111110' . '101111101110101110101110111110' x ($n/30);
until (($s = $i*$i) > $n) { $d = $i<<1; do { substr($_, $s, 1, '1') } until ($s += $d) > $n; 1 while substr($_, $i += 2, 1); } $_ = substr($_, 1, $n); # For just the count: return ($_ =~ tr/0//); push @primes, pos while m/0/g; @primes;
}</lang>
This older version uses half the memory, but at the expense of a bit of speed and code complexity: <lang perl>sub dj_string {
my($end) = @_; return @{([],[],[2],[2,3],[2,3])[$end]} if $end <= 4; $end-- if ($end & 1) == 0; my $s_end = $end >> 1;
my $whole = int( ($end>>1) / 15); # prefill with 3 and 5 marked my $sieve = '100010010010110' . '011010010010110' x $whole; substr($sieve, ($end>>1)+1) = ; my ($n, $limit, $s) = ( 7, int(sqrt($end)), 0 ); while ( $n <= $limit ) { for ($s = ($n*$n) >> 1; $s <= $s_end; $s += $n) { substr($sieve, $s, 1) = '1'; } do { $n += 2 } while substr($sieve, $n>>1, 1); } # If you just want the count, it's very fast: # my $count = 1 + $sieve =~ tr/0//; my @primes = (2); push @primes, 2*pos($sieve)-1 while $sieve =~ m/0/g; @primes;
}</lang>
Experimental
These are examples of golfing or unusual styles.
Golfing a bit, at the expense of speed: <lang perl>sub sieve{ my (@s, $i); grep { not $s[ $i = $_ ] and do { $s[ $i += $_ ]++ while $i <= $_[0]; 1 } } 2..$_[0] }
print join ", " => sieve 100;</lang>
Or with bit strings (much slower than the vector version above): <lang perl>sub sieve{ my ($s, $i); grep { not vec $s, $i = $_, 1 and do { (vec $s, $i += $_, 1) = 1 while $i <= $_[0]; 1 } } 2..$_[0] }
print join ", " => sieve 100;</lang>
A short recursive version: <lang perl>sub erat {
my $p = shift; return $p, $p**2 > $_[$#_] ? @_ : erat(grep $_%$p, @_)
}
print join ', ' => erat 2..100000;</lang>
Regexp (purely an example -- the regex engine limits it to only 32769):<lang perl>sub sieve { my ($s, $p) = "." . ("x" x shift);
1 while ++$p and $s =~ /^(.{$p,}?)x/g and $p = length($1) and $s =~ s/(.{$p})./${1}./g and substr($s, $p, 1) = "x"; $s }
print sieve(1000);</lang>
Extensible sieves
Here are two incremental versions, which allows one to create a tied array of primes: <lang perl>use strict; use warnings; package Tie::SieveOfEratosthenes;
sub TIEARRAY { my $class = shift; bless \$class, $class; }
- If set to true, produces copious output. Observing this output
- is an excellent way to gain insight into how the algorithm works.
use constant DEBUG => 0;
- If set to true, causes the code to skip over even numbers,
- improving runtime. It does not alter the output content, only the speed.
use constant WHEEL2 => 0;
BEGIN {
# This is loosely based on the Python implementation of this task, # specifically the "Infinite generator with a faster algorithm"
my @primes = (2, 3); my $ps = WHEEL2 ? 1 : 0; my $p = $primes[$ps]; my $q = $p*$p; my $incr = WHEEL2 ? 2 : 1; my $candidate = $primes[-1] + $incr; my %sieve;
print "Initial: p = $p, q = $q, candidate = $candidate\n" if DEBUG;
sub FETCH { my $n = pop; return if $n < 0; return $primes[$n] if $n <= $#primes; OUTER: while( 1 ) {
# each key in %sieve is a composite number between # p and p-squared. Each value in %sieve is $incr x the prime # which acted as a 'seed' for that key. We use the value # to step through multiples of the seed-prime, until we find # an empty slot in %sieve. while( my $s = delete $sieve{$candidate} ) { print "$candidate a multiple of ".($s/$incr).";\t\t" if DEBUG; my $composite = $candidate + $s; $composite += $s while exists $sieve{$composite}; print "The next stored multiple of ".($s/$incr)." is $composite\n" if DEBUG; $sieve{$composite} = $s; $candidate += $incr; }
print "Candidate $candidate is not in sieve\n" if DEBUG;
while( $candidate < $q ) { print "$candidate is prime\n" if DEBUG; push @primes, $candidate; $candidate += $incr; next OUTER if exists $sieve{$candidate}; }
die "Candidate = $candidate, p = $p, q = $q" if $candidate > $q; print "Candidate $candidate is equal to $p squared;\t" if DEBUG;
# Thus, it is now time to add p to the sieve, my $step = $incr * $p; my $composite = $q + $step; $composite += $step while exists $sieve{$composite}; print "The next multiple of $p is $composite\n" if DEBUG; $sieve{$composite} = $step;
# and fetch out a new value for p from our primes array. $p = $primes[++$ps]; $q = $p * $p;
# And since $candidate was equal to some prime squared, # it's obviously composite, and we need to increment it. $candidate += $incr; print "p is $p, q is $q, candidate is $candidate\n" if DEBUG; } continue { return $primes[$n] if $n <= $#primes; } }
}
if( !caller ) { tie my (@prime_list), 'Tie::SieveOfEratosthenes'; my $limit = $ARGV[0] || 100; my $line = ""; for( my $count = 0; $prime_list[$count] < $limit; ++$count ) { $line .= $prime_list[$count]. ", "; next if length($line) <= 70; if( $line =~ tr/,// > 1 ) { $line =~ s/^(.*,) (.*, )/$2/; print $1, "\n"; } else { print $line, "\n"; $line = ""; } } $line =~ s/, \z//; print $line, "\n" if $line; }
1;</lang> This one is based on the vector sieve shown earlier, but adds to a list as needed, just sieving in the segment. Slightly faster and half the memory vs. the previous incremental sieve. It uses the same API -- arguably we should be offset by one so $primes[$n] returns the $n'th prime. <lang perl>use strict; use warnings; package Tie::SieveOfEratosthenes;
sub TIEARRAY {
my $class = shift; my @primes = (2,3,5,7); return bless \@primes, $class;
}
sub prextend { # Extend the given list of primes using a segment sieve
my($primes, $to) = @_; $to-- unless $to & 1; # Ensure end is odd return if $to < $primes->[-1]; my $sqrtn = int(sqrt($to)+0.001); prextend($primes, $sqrtn) if $primes->[-1] < $sqrtn; my($segment, $startp) = (, $primes->[-1]+1); my($s_beg, $s_len) = ($startp >> 1, ($to>>1) - ($startp>>1)); for my $p (@$primes) { last if $p > $sqrtn; if ($p >= 3) { my $p2 = $p*$p; if ($p2 < $startp) { # Bump up to next odd multiple of p >= startp my $f = 1+int(($startp-1)/$p); $p2 = $p * ($f | 1); } for (my $s = ($p2>>1)-$s_beg; $s <= $s_len; $s += $p) { vec($segment, $s, 1) = 1; # Mark composites in segment } } } # Now add all the primes found in the segment to the list do { push @$primes, 1+2*($_+$s_beg) if !vec($segment,$_,1) } for 0 .. $s_len;
}
sub FETCHSIZE { 0x7FFF_FFFF } # Allows foreach to work sub FETCH {
my($primes, $n) = @_; return if $n < 0; # Keep expanding the list as necessary, 5% larger each time. prextend($primes, 1000+int(1.05*$primes->[-1])) while $n > $#$primes; return $primes->[$n];
}
if( !caller ) {
tie my @prime_list, 'Tie::SieveOfEratosthenes'; my $limit = $ARGV[0] || 100; print $prime_list[0]; my $i = 1; while ($prime_list[$i] < $limit) { print " ", $prime_list[$i++]; } print "\n";
}
1;</lang>
Perl 6
<lang perl6>sub sieve( Int $limit ) {
my @is-prime = False, False, True xx $limit - 1;
gather for @is-prime.kv -> $number, $is-prime { if $is-prime { take $number; @is-prime[$_] = False if $_ %% $number for $number**2 .. $limit; } }
}
(sieve 100).join(",").say;</lang>
A recursive version: <lang perl6>multi erat(Int $N) { erat 2 .. $N } multi erat(@a where @a[0] > sqrt @a[*-1]) { @a } multi erat(@a) { @a[0], erat(@a.grep: * % @a[0]) } say erat 100;</lang>
Of course, upper limits are for wusses. Here's a version using infinite streams, that just keeps going until you ^C it (works under Niecza):
<lang perl6>role Filter[Int $factor] {
method next { repeat until $.value % $factor { callsame } }
}
class Stream {
has Int $.value is rw = 1; method next { ++$.value } method filter { self but Filter[$.value] }
}
.next, say .value for Stream.new, *.filter ... *;</lang>
PHP
<lang php> function iprimes_upto($limit) {
for ($i = 2; $i < $limit; $i++) {
$primes[$i] = true;
} for ($n = 2; $n < $limit; $n++) {
if ($primes[$n]) { for ($i = $n*$n; $i < $limit; $i += $n) { $primes[$i] = false; } }
} return $primes;
} </lang>
PicoLisp
<lang PicoLisp>(de sieve (N)
(let Sieve (range 1 N) (set Sieve) (for I (cdr Sieve) (when I (for (S (nth Sieve (* I I)) S (nth (cdr S) I)) (set S) ) ) ) (filter bool Sieve) ) )</lang>
Output:
: (sieve 100) -> (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)
PL/I
<lang pli>eratos: proc options (main) reorder;
dcl i fixed bin (31); dcl j fixed bin (31); dcl n fixed bin (31); dcl sn fixed bin (31);
dcl hbound builtin; dcl sqrt builtin;
dcl sysin file; dcl sysprint file;
get list (n); sn = sqrt(n);
begin;
dcl primes(n) bit (1) aligned init ((*)((1)'1'b));
i = 2;
do while(i <= sn); do j = i ** 2 by i to hbound(primes, 1); /* Adding a test would just slow down processing! */ primes(j) = '0'b; end;
do i = i + 1 to sn until(primes(i)); end; end;
do i = 2 to hbound(primes, 1); if primes(i) then put data(i); end;
end; end eratos;</lang>
Pop11
define eratostenes(n); lvars bits = inits(n), i, j; for i from 2 to n do if bits(i) = 0 then printf('' >< i, '%s\n'); for j from 2*i by i to n do 1 -> bits(j); endfor; endif; endfor; enddefine;
PowerShell
Basic procedure
It outputs immediately so that the number can be used by the pipeline. <lang PowerShell>function Sieve ( [int] $num ) {
$isprime = @{} 2..$num | Where-Object { $isprime[$_] -eq $null } | ForEach-Object { $_ $isprime[$_] = $true $i=$_*$_ for ( ; $i -le $num; $i += $_ ) { $isprime[$i] = $false } }
}</lang>
Processing
Calculate the primes up to 1000000 with Processing, including a visualisation of the process. As an additional visual effect, the layout of the pixel could be changed from the line-by-line layout to a spiral-like layout starting in the middle of the screen. <lang java>int maxx,maxy; int max; boolean[] sieve;
void plot(int pos, boolean active) {
set(pos%maxx,pos/maxx, active?#000000:#ffffff);
}
void setup() {
size(1000, 1000, P2D); frameRate(2); maxx=width; maxy=height; max=width*height; sieve=new boolean[max+1]; sieve[1]=false; plot(0,false); plot(1,false); for(int i=2;i<=max;i++) { sieve[i]=true; plot(i,true); }
}
int i=2;
void draw() {
if(!sieve[i]) { while(i*i<max && !sieve[i]) { i++; } } if(sieve[i]) { print(i+" "); for(int j=i*i;j<=max;j+=i) { if(sieve[j]) { sieve[j]=false; plot(j,false); } } } if(i*i<max) { i++; } else { noLoop(); println("finished"); }
}</lang>
Prolog
Using facts to record composite numbers
The first two solutions use Prolog "facts" to record the composite (i.e. already-visited) numbers.
Elementary approach: multiplication-free, division-free, mod-free, and cut-free
The basic Eratosthenes sieve depends on nothing more complex than counting. In celebration of this simplicity, the first approach to the problem taken here is free of multiplication and division, as well as Prolog's non-logical "cut".
It defines the predicate between/4 to avoid division, and composite/1 to record integers that are found to be composite.
<lang Prolog>% %sieve( +N, -Primes ) is true if Primes is the list of consecutive primes % that are less than or equal to N sieve( N, [2|Rest]) :-
retractall( composite(_) ), sieve( N, 2, Rest ) -> true. % only one solution
% sieve P, find the next non-prime, and then recurse: sieve( N, P, [I|Rest] ) :-
sieve_once(P, N), (P = 2 -> P2 is P+1; P2 is P+2), between(P2, N, I), (composite(I) -> fail; sieve( N, I, Rest )).
% It is OK if there are no more primes less than or equal to N: sieve( N, P, [] ).
sieve_once(P, N) :-
forall( between(P, N, P, IP), (composite(IP) -> true ; assertz( composite(IP) )) ).
% To avoid division, we use the iterator
% between(+Min, +Max, +By, -I)
% where we assume that By > 0
% This is like "for(I=Min; I <= Max; I+=By)" in C.
between(Min, Max, By, I) :-
Min =< Max, A is Min + By, (I = Min; between(A, Max, By, I) ).
% Some Prolog implementations require the dynamic predicates be
% declared:
- - dynamic( composite/1 ).
</lang> The above has been tested with SWI-Prolog and gprolog.
<lang Prolog>% SWI-Prolog:
?- time( (sieve(100000,P), length(P,N), writeln(N), last(P, LP), writeln(LP) )). % 1,323,159 inferences, 0.862 CPU in 0.921 seconds (94% CPU, 1534724 Lips) P = [2, 3, 5, 7, 11, 13, 17, 19, 23|...], N = 9592, LP = 99991. </lang>
Optimized approach
<lang Prolog>sieve(N, [2|PS]) :- % PS is list of odd primes up to N
retractall(mult(_)), sieve_O(3,N,PS).
sieve_O(I,N,PS) :- % sieve odds from I up to N to get PS
I =< N, !, I1 is I+2, ( mult(I) -> sieve_O(I1,N,PS) ; ( I =< N / I -> ISq is I*I, DI is 2*I, add_mults(DI,ISq,N) ; true ), PS = [I|T], sieve_O(I1,N,T) ).
sieve_O(I,N,[]) :- I > N.
add_mults(DI,I,N) :-
I =< N, !, ( mult(I) -> true ; assert(mult(I)) ), I1 is I+DI, add_mults(DI,I1,N).
add_mults(_,I,N) :- I > N.
main(N) :- current_prolog_flag(verbose,F),
set_prolog_flag(verbose,normal), time( sieve( N,P)), length(P,Len), last(P, LP), writeln([Len,LP]), set_prolog_flag(verbose,F).
- - dynamic( mult/1 ).
- - main(100000), main(1000000).</lang>
Running it produces
<lang Prolog>%% stdout copy [9592, 99991] [78498, 999983]
%% stderr copy % 293,176 inferences, 0.14 CPU in 0.14 seconds (101% CPU, 2094114 Lips) % 3,122,303 inferences, 1.63 CPU in 1.67 seconds (97% CPU, 1915523 Lips)</lang>
which indicates ~ N1.1 empirical orders of growth, which is consistent with the O(N log log N) theoretical runtime complexity.
Using a priority queue
Uses a ariority queue, from the paper "The Genuine Sieve of Eratosthenes" by Melissa O'Neill. Works with YAP (Yet Another Prolog)
<lang Prolog>?- use_module(library(heaps)).
prime(2). prime(N) :- prime_heap(N, _).
prime_heap(3, H) :- list_to_heap([9-6], H). prime_heap(N, H) :-
prime_heap(M, H0), N0 is M + 2, next_prime(N0, H0, N, H).
next_prime(N0, H0, N, H) :-
\+ min_of_heap(H0, N0, _), N = N0, Composite is N*N, Skip is N+N, add_to_heap(H0, Composite, Skip, H).
next_prime(N0, H0, N, H) :-
min_of_heap(H0, N0, _), adjust_heap(H0, N0, H1), N1 is N0 + 2, next_prime(N1, H1, N, H).
adjust_heap(H0, N, H) :-
min_of_heap(H0, N, _), get_from_heap(H0, N, Skip, H1), Composite is N + Skip, add_to_heap(H1, Composite, Skip, H2), adjust_heap(H2, N, H).
adjust_heap(H, N, H) :-
\+ min_of_heap(H, N, _).</lang>
PureBasic
Basic procedure
<lang PureBasic>For n=2 To Sqr(lim)
If Nums(n)=0 m=n*n While m<=lim Nums(m)=1 m+n Wend EndIf
Next n</lang>
Working example
<lang PureBasic>Dim Nums.i(0) Define l, n, m, lim
If OpenConsole()
; Ask for the limit to search, get that input and allocate a Array Print("Enter limit for this search: ") lim=Val(Input()) ReDim Nums(lim) ; Use a basic Sieve of Eratosthenes For n=2 To Sqr(lim) If Nums(n)=#False m=n*n While m<=lim Nums(m)=#True m+n Wend EndIf Next n ;Present the result to our user PrintN(#CRLF$+"The Prims up to "+Str(lim)+" are;") m=0: l=Log10(lim)+1 For n=2 To lim If Nums(n)=#False Print(RSet(Str(n),l)+" ") m+1 If m>72/(l+1) m=0: PrintN("") EndIf EndIf Next Print(#CRLF$+#CRLF$+"Press ENTER to exit"): Input() CloseConsole()
EndIf</lang>
Output may look like;
Enter limit for this search: 750 The Prims up to 750 are; 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 131 137 139 149 151 157 163 167 173 179 181 191 193 197 199 211 223 227 229 233 239 241 251 257 263 269 271 277 281 283 293 307 311 313 317 331 337 347 349 353 359 367 373 379 383 389 397 401 409 419 421 431 433 439 443 449 457 461 463 467 479 487 491 499 503 509 521 523 541 547 557 563 569 571 577 587 593 599 601 607 613 617 619 631 641 643 647 653 659 661 673 677 683 691 701 709 719 727 733 739 743 Press ENTER to exit
Python
Note that the examples use range instead of xrange for Python 3 and Python 2 compatability, but when using Python 2 xrange is the nearest equivalent to Python 3's implementation of range and should be substituted for ranges with a very large number of items.
Using set lookup
The version below uses a set to store the multiples. set objects are much faster (usually O(log n)) than lists (O(n)) for checking if a given object is a member. Using the set.update method avoids explicit iteration in the interpreter, giving a further speed improvement.
<lang python>def eratosthenes2(n):
multiples = set() for i in range(2, n+1): if i not in multiples: yield i multiples.update(range(i*i, n+1, i))
print(list(eratosthenes2(100)))</lang>
Using array lookup
The version below uses array lookup to test for primality. The function primes_upto() is a straightforward implementation of Sieve of Eratosthenesalgorithm. It returns prime numbers less than or equal to limit. <lang python>def primes_upto(limit):
is_prime = [False] * 2 + [True] * (limit - 1) for n in range(int(limit**0.5 + 1.5)): # stop at ``sqrt(limit)`` if is_prime[n]: for i in range(n*n, limit+1, n): is_prime[i] = False return [i for i, prime in enumerate(is_prime) if prime]</lang>
Using generator
The following code may be slightly slower than using the array/list as above, but uses no memory for output: <lang python>def iprimes_upto(limit):
is_prime = [False] * 2 + [True] * (limit - 1) for n in xrange(int(limit**0.5 + 1.5)): # stop at ``sqrt(limit)`` if is_prime[n]: for i in range(n * n, limit + 1, n): # start at ``n`` squared is_prime[i] = False for i in xrange(limit + 1):
if is_prime[i]: yield i</lang>
- Example:
<lang python>>>> list(iprimes_upto(15))
[2, 3, 5, 7, 11, 13]</lang>
Odds-only version of the array sieve above
The following code is faster than the above array version using only odd composite operations (for a factor of over two) and because it has been optimized to use slice operations for composite number culling to avoid extra work by the interpreter: <lang python>def primes2(limit):
if limit < 2: return [] if limit < 3: return [2] lmtbf = (limit - 3) // 2 buf = [True] * (lmtbf + 1) for i in range((int(limit ** 0.5) - 3) // 2 + 1): if buf[i]: p = i + i + 3 s = p * (i + 1) + i buf[s::p] = [False] * ((lmtbf - s) // p + 1) return [2] + [i + i + 3 for i, v in enumerate(buf) if v]</lang>
Note that "range" needs to be changed to "xrange" for maximum speed with Python 2.
Odds-only version of the generator version above
The following code is faster than the above generator version using only odd composite operations (for a factor of over two) and because it has been optimized to use slice operations for composite number culling to avoid extra work by the interpreter:
<lang python>def iprimes2(limit):
yield 2 if limit < 3: return lmtbf = (limit - 3) // 2 buf = [True] * (lmtbf + 1) for i in range((int(limit ** 0.5) - 3) // 2 + 1): if buf[i]: p = i + i + 3 s = p * (i + 1) + i buf[s::p] = [False] * ((lmtbf - s) // p + 1) for i in range(lmtbf + 1): if buf[i]: yield (i + i + 3)</lang>
Note that this version may actually run slightly faster than the equivalent array version with the advantage that the output doesn't require any memory.
Also note that "range" needs to be changed to "xrange" for maximum speed with Python 2.
Factorization wheel235 version of the generator version
This uses a 235 factorial wheel for further reductions in operations; the same techniques can be applied to the array version as well; it runs slightly faster and uses slightly less memory as compared to the odds-only algorithms:
<lang python>def primes235(limit):
yield 2; yield 3; yield 5 if limit < 7: return modPrms = [7,11,13,17,19,23,29,31] gaps = [4,2,4,2,4,6,2,6,4,2,4,2,4,6,2,6] # 2 loops for overflow ndxs = [0,0,0,0,1,1,2,2,2,2,3,3,4,4,4,4,5,5,5,5,5,5,6,6,7,7,7,7,7,7] lmtbf = (limit + 23) // 30 * 8 - 1 # integral number of wheels rounded up lmtsqrt = (int(limit ** 0.5) - 7) lmtsqrt = lmtsqrt // 30 * 8 + ndxs[lmtsqrt % 30] # round down on the wheel buf = [True] * (lmtbf + 1) for i in range(lmtsqrt + 1): if buf[i]: ci = i & 7; p = 30 * (i >> 3) + modPrms[ci] s = p * p - 7; p8 = p << 3 for j in range(8): c = s // 30 * 8 + ndxs[s % 30] buf[c::p8] = [False] * ((lmtbf - c) // p8 + 1) s += p * gaps[ci]; ci += 1 for i in range(lmtbf - 6 + (ndxs[(limit - 7) % 30])): # adjust for extras if buf[i]: yield (30 * (i >> 3) + modPrms[i & 7])</lang>
Note: Much of the time (almost two thirds for this last case for Python 2.7.6) for any of these array/list or generator algorithms is used in the computation and enumeration of the final output in the last line(s), so any slight changes to those lines can greatly affect execution time. For Python 3 this enumeration is about twice as slow as Python 2 (Python 3.3 slow and 3.4 slower) for an even bigger percentage of time spent just outputting the results. This slow enumeration means that there is little advantage to versions that use even further wheel factorization, as the composite number culling is a small part of the time to enumerate the results.
If just the count of the number of primes over a range is desired, then converting the functions to prime counting functions by changing the final enumeration lines to "return buf.count(True)" will save a lot of time.
Note that "range" needs to be changed to "xrange" for maximum speed with Python 2 where Python 2's "xrange" is a better choice for very large sieve ranges.
Timings were done primarily in Python 2 although source is Python 2/3 compatible (shows range and not xrange).
Using numpy
Below code adapted from literateprograms.org using numpy <lang python>import numpy def primes_upto2(limit):
is_prime = numpy.ones(limit + 1, dtype=numpy.bool) for n in xrange(2, int(limit**0.5 + 1.5)): if is_prime[n]: is_prime[n*n::n] = 0 return numpy.nonzero(is_prime)[0][2:]</lang>
Performance note: there is no point to add wheels here, due to execution of p[n*n::n] = 0 and nonzero() takes us almost all time.
Also see Prime numbers and Numpy – Python.
Using wheels with numpy
Version with wheel based optimization: <lang python>from numpy import array, bool_, multiply, nonzero, ones, put, resize
def makepattern(smallprimes):
pattern = ones(multiply.reduce(smallprimes), dtype=bool_) pattern[0] = 0 for p in smallprimes: pattern[p::p] = 0 return pattern
def primes_upto3(limit, smallprimes=(2,3,5,7,11)):
sp = array(smallprimes) if limit <= sp.max(): return sp[sp <= limit] # isprime = resize(makepattern(sp), limit + 1) isprime[:2] = 0; put(isprime, sp, 1) # for n in range(sp.max() + 2, int(limit**0.5 + 1.5), 2): if isprime[n]: isprime[n*n::n] = 0 return nonzero(isprime)[0]</lang>
Examples: <lang python>>>> primes_upto3(10**6, smallprimes=(2,3)) # Wall time: 0.17 array([ 2, 3, 5, ..., 999961, 999979, 999983]) >>> primes_upto3(10**7, smallprimes=(2,3)) # Wall time: 2.13 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto3(15) array([ 2, 3, 5, 7, 11, 13]) >>> primes_upto3(10**7, smallprimes=primes_upto3(15)) # Wall time: 1.31 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto2(10**7) # Wall time: 1.39 <-- version without wheels array([ 2, 3, 5, ..., 9999971, 9999973, 9999991]) >>> primes_upto3(10**7) # Wall time: 1.30 array([ 2, 3, 5, ..., 9999971, 9999973, 9999991])</lang> The above-mentioned examples demonstrate that the given wheel based optimization does not show significant performance gain.
Infinite generator
A generator that will generate primes indefinitely (perhaps until it runs out of memory). Used as a library here.
<lang python>import heapq
- generates all prime numbers
def sieve():
# priority queue of the sequences of non-primes # the priority queue allows us to get the "next" non-prime quickly nonprimes = [] i = 2 while True: if nonprimes and i == nonprimes[0][0]: # non-prime while nonprimes[0][0] == i: # for each sequence that generates this number, # have it go to the next number (simply add the prime) # and re-position it in the priority queue x = nonprimes[0] x[0] += x[1] heapq.heapreplace(nonprimes, x) else: # prime # insert a 2-element list into the priority queue: # [current multiple, prime] # the first element allows sorting by value of current multiple # we start with i^2 heapq.heappush(nonprimes, [i*i, i]) yield i i += 1</lang>
Example:
>>> foo = sieve() >>> for i in range(8): ... print(next(foo)) ... 2 3 5 7 11 13 17 19
Infinite generator with a faster algorithm
The adding of each discovered prime's incremental step info to the mapping should be postponed until the prime's square is seen amongst the candidate numbers, as it is useless before that point. This drastically reduces the space complexity from O(n) to O(sqrt(n/log(n))), in n
primes produced, and also lowers the run time complexity quite low (this test entry in Python 2.7 and this test entry in Python 3.x shows about ~ n1.08 empirical order of growth which is very close to the theoretical value of O(n log(n) log(log(n))), in n
primes produced):
<lang python>def primes():
yield 2; yield 3; yield 5; yield 7; bps = (p for p in primes()) # separate supply of "base" primes (b.p.) p = next(bps) and next(bps) # discard 2, then get 3 q = p * p # 9 - square of next base prime to keep track of, sieve = {} # in the sieve dict n = 9 # n is the next candidate number while True: if n not in sieve: # n is not a multiple of any of base primes, if n < q: # below next base prime's square, so yield n # n is prime else: p2 = p + p # n == p * p: for prime p, add p * p + 2 * p sieve[q + p2] = p2 # to the dict, with 2 * p as the increment step p = next(bps); q = p * p # pull next base prime, and get its square else: s = sieve.pop(n); nxt = n + s # n's a multiple of some b.p., find next multiple while nxt in sieve: nxt += s # ensure each entry is unique sieve[nxt] = s # nxt is next non-marked multiple of this prime n += 2 # work on odds only
import itertools def primes_up_to(limit):
return list(itertools.takewhile(lambda p: p <= limit, primes()))</lang>
Fast infinite generator using a wheel
Although theoretically over three times faster than odds-only, the following code using a 2/3/5/7 wheel is only about 1.5 times faster than the above odds-only code due to the extra overheads in code complexity. The test link for Python 2.7 and test link for Python 3.x show about the same empirical order of growth as the odds-only implementation above once the range grows enough so the dict operations become amortized to a constant factor.
<lang python>def primes():
for p in [2,3,5,7]: yield p # base wheel primes gaps1 = [ 2,4,2,4,6,2,6,4,2,4,6,6,2,6,4,2,6,4,6,8,4,2,4,2,4,8 ] gaps = gaps1 + [ 6,4,6,2,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10 ] # wheel2357 def wheel_prime_pairs(): yield (11,0); bps = wheel_prime_pairs() # additional primes supply p, pi = next(bps); q = p * p # adv to get 11 sqr'd is 121 as next square to put sieve = {}; n = 13; ni = 1 # into sieve dict; init cndidate, wheel ndx while True: if n not in sieve: # is not a multiple of previously recorded primes if n < q: yield (n, ni) # n is prime with wheel modulo index else: npi = pi + 1 # advance wheel index if npi > 47: npi = 0 sieve[q + p * gaps[pi]] = (p, npi) # n == p * p: put next cull position on wheel p, pi = next(bps); q = p * p # advance next prime and prime square to put else: s, si = sieve.pop(n) nxt = n + s * gaps[si] # move current cull position up the wheel si = si + 1 # advance wheel index if si > 47: si = 0 while nxt in sieve: # ensure each entry is unique by wheel nxt += s * gaps[si] si = si + 1 # advance wheel index if si > 47: si = 0 sieve[nxt] = (s, si) # next non-marked multiple of a prime nni = ni + 1 # advance wheel index if nni > 47: nni = 0 n += gaps[ni]; ni = nni # advance on the wheel for p, pi in wheel_prime_pairs(): yield p # strip out indexes</lang>
Further gains of about 1.5 times in speed can be made using the same code by only changing the tables and a few constants for a further constant factor gain of about 1.5 times in speed by using a 2/3/5/7/11/13/17 wheel (with the gaps list 92160 elements long) computed for a slight constant overhead time as per the test link for Python 2.7 and test link for Python 3.x. Further wheel factorization will not really be worth it as the gains will be small (if any and not losses) and the gaps table huge - it is already too big for efficient use by 32-bit Python 3 and the wheel should likely be stopped at 13: <lang python>def primes():
whlPrms = [2,3,5,7,11,13,17] # base wheel primes for p in whlPrms: yield p def makeGaps(): buf = [True] * (3 * 5 * 7 * 11 * 13 * 17 + 1) # all odds plus extra for o/f for p in whlPrms: if p < 3: continue # no need to handle evens strt = (p * p - 19) >> 1 # start position (divided by 2 using shift) while strt < 0: strt += p buf[strt::p] = [False] * ((len(buf) - strt - 1) // p + 1) # cull for p whlPsns = [i + i for i,v in enumerate(buf) if v] return [whlPsns[i + 1] - whlPsns[i] for i in range(len(whlPsns) - 1)] gaps = makeGaps() # big wheel gaps def wheel_prime_pairs(): yield (19,0); bps = wheel_prime_pairs() # additional primes supply p, pi = next(bps); q = p * p # adv to get 11 sqr'd is 121 as next square to put sieve = {}; n = 23; ni = 1 # into sieve dict; init cndidate, wheel ndx while True: if n not in sieve: # is not a multiple of previously recorded primes if n < q: yield (n, ni) # n is prime with wheel modulo index else: npi = pi + 1 # advance wheel index if npi > 92159: npi = 0 sieve[q + p * gaps[pi]] = (p, npi) # n == p * p: put next cull position on wheel p, pi = next(bps); q = p * p # advance next prime and prime square to put else: s, si = sieve.pop(n) nxt = n + s * gaps[si] # move current cull position up the wheel si = si + 1 # advance wheel index if si > 92159: si = 0 while nxt in sieve: # ensure each entry is unique by wheel nxt += s * gaps[si] si = si + 1 # advance wheel index if si > 92159: si = 0 sieve[nxt] = (s, si) # next non-marked multiple of a prime nni = ni + 1 # advance wheel index if nni > 92159: nni = 0 n += gaps[ni]; ni = nni # advance on the wheel for p, pi in wheel_prime_pairs(): yield p # strip out indexes
</lang>
R
This code is vectorised, so it runs quickly for n < one million. The allocation of the primes vector means memory usage is too high for n much larger than that.<lang R>sieve <- function(n) {
n <- as.integer(n) if(n > 1e6) stop("n too large") primes <- rep(TRUE, n) primes[1] <- FALSE last.prime <- 2L for(i in last.prime:floor(sqrt(n))) { primes[seq.int(2L*last.prime, n, last.prime)] <- FALSE last.prime <- last.prime + min(which(primes[(last.prime+1):n])) } which(primes)
}
sieve(1000)</lang>
Racket
Imperative versions
Ugly imperative version: <lang Racket>#lang racket
(define (sieve n)
(define non-primes '()) (define primes '()) (for ([i (in-range 2 (add1 n))]) (unless (member i non-primes) (set! primes (cons i primes)) (for ([j (in-range (* i i) (add1 n) i)]) (set! non-primes (cons j non-primes))))) (reverse primes))
(sieve 100)</lang>
A little nicer, but still imperative: <lang Racket>#lang racket (define (sieve n)
(define primes (make-vector (add1 n) #t)) (for* ([i (in-range 2 (add1 n))] #:when (vector-ref primes i) [j (in-range (* i i) (add1 n) i)]) (vector-set! primes j #f)) (for/list ([n (in-range 2 (add1 n))] #:when (vector-ref primes n)) n))
(sieve 100)</lang>
Imperative version using a bit vector: <lang Racket>#lang racket (require data/bit-vector)
- Returns a list of prime numbers up to natural number limit
(define (eratosthenes limit)
(define bv (make-bit-vector (+ limit 1) #f)) (bit-vector-set! bv 0 #t) (bit-vector-set! bv 1 #t) (for* ([i (in-range (sqrt limit))] #:unless (bit-vector-ref bv i) [j (in-range (* 2 i) (bit-vector-length bv) i)]) (bit-vector-set! bv j #t)) ;; translate to a list of primes (for/list ([i (bit-vector-length bv)] #:unless (bit-vector-ref bv i)) i))
(eratosthenes 100) </lang>
- Output:
'(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)
Infinite list of primes
Using laziness
These examples use infinite lists (streams) to implement the sieve of Eratosthenes in a functional way, and producing all prime numbers. The following functions are used as a prefix for pieces of code that follow:
<lang Racket>#lang lazy (define (ints-from i d) (cons i (ints-from (+ i d) d))) (define (after n l f)
(if (< (car l) n) (cons (car l) (after n (cdr l) f)) (f l)))
(define (diff l1 l2)
(let ([x1 (car l1)] [x2 (car l2)]) (cond [(< x1 x2) (cons x1 (diff (cdr l1) l2 ))] [(> x1 x2) (diff l1 (cdr l2)) ] [else (diff (cdr l1) (cdr l2)) ])))
(define (union l1 l2) ; union of two lists
(let ([x1 (car l1)] [x2 (car l2)]) (cond [(< x1 x2) (cons x1 (union (cdr l1) l2 ))] [(> x1 x2) (cons x2 (union l1 (cdr l2)))] [else (cons x1 (union (cdr l1) (cdr l2)))])))</lang>
Basic sieve
<lang Racket>(define (sieve l)
(define x (car l)) (cons x (sieve (diff (cdr l) (ints-from (+ x x) x)))))
(define primes (sieve (ints-from 2 1))) (!! (take 25 primes))</lang>
Runs at ~ n^2.1 empirically, for n <= 1500 primes produced.
With merged composites
Note that the first numbers, 2 and 3, are done manually outside of the sieve loop -- to ensure that the non-primes list is never empty, which simplifies the code since all functions assume non-empty infinite lists. This also means that the even numbers are ignored.
<lang Racket>(define (sieve l non-primes)
(let ([x (car l)] [np (car non-primes)]) (cond [(= x np) (sieve (cdr l) (cdr non-primes))] ; else x < np [else (cons x (sieve (cdr l) (union (ints-from (* x x) (* 2 x)) non-primes)))])))
(define primes (cons 2 (cons 3 (sieve (ints-from 5 2) (ints-from 9 6)))))</lang>
Runs at ~ n^2.4 and worse, empirically, producing up to n=700
primes.
Basic sieve Optimized with postponed processing
Since a prime's multiples that count start from its square, we should only start removing them when we reach that square. <lang Racket>(define (sieve l prs)
(define p (car prs)) (define q (* p p)) (after q l (λ(t) (sieve (diff t (ints-from q p)) (cdr prs)))))
(define primes (cons 2 (sieve (ints-from 3 1) primes)))</lang>
Runs at ~ n^1.4 up to n=10,000. The initial 2 in the self-referential primes definition is needed to prevent a "black hole".
Merged composites Optimized with postponed processing
Since a prime's multiples that matter start from its square, we should only add them when we reach that square. Improves the time complexity to ~ n^1.5, tested producing up to n=6,000
primes.
<lang Racket>(define (composites l q primes)
(after q l (λ(t) (let ([p (car primes)] [r (cadr primes)]) (composites (union t (ints-from q (* 2 p))) ; q = p*p (* r r) (cdr primes))))))
(define primes (cons 2 (cons 3 (diff (ints-from 5 2)
(composites (ints-from 9 6) 25 (cddr primes))))))</lang>
The processing starts from odds only, ignoring all evens above 2, so effectively it is employing the simplest kind of a wheel.
Implementation of Richard Bird's algorithm
Appears in M.O'Neill's paper. Achieves on its own the proper postponement that is specifically arranged for in the version above (with after
), and is yet more efficient, because it folds to the right and so builds the right-leaning structure of merges at run time, where the more frequently-producing streams of multiples appear higher in that structure, so the composite numbers produced by them have less merge
nodes to percolate through:
<lang Racket>(define primes
(cons 2 (diff (ints-from 3 1) (foldr (λ(p r) (define q (* p p)) (cons q (union (ints-from (+ q p) p) r))) '() primes))))</lang>
Using threads and channels
Same algorithm as "merged composites" above (without the optimization), but now using threads and channels to produce a channel of all prime numbers (similar to newsqueak). The macro at the top is a convenient wrapper around definitions of channels using a thread that feeds them.
<lang Racket>#lang racket (define-syntax (define-thread-loop stx)
(syntax-case stx () [(_ (name . args) expr ...) (with-syntax ([out! (datum->syntax stx 'out!)]) #'(define (name . args) (define out (make-channel)) (define (out! x) (channel-put out x)) (thread (λ() (let loop () expr ... (loop)))) out))]))
(define-thread-loop (ints-from i d) (out! i) (set! i (+ i d))) (define-thread-loop (merge c1 c2)
(let loop ([x1 (channel-get c1)] [x2 (channel-get c2)]) (cond [(< x1 x2) (out! x1) (loop (channel-get c1) x2)] [(> x1 x2) (out! x2) (loop x1 (channel-get c2))] [else (out! x1) (loop (channel-get c1) (channel-get c2))])))
(define-thread-loop (sieve l non-primes)
(let loop ([x (channel-get l)] [np (channel-get non-primes)]) (cond [(< np x) (loop x (channel-get non-primes))] [(= np x) (loop (channel-get l) (channel-get non-primes))] [else (out! x) (set! non-primes (merge (ints-from (* x x) x) non-primes))])))
(define-thread-loop (cons x l)
(out! x) (let loop () (out! (channel-get l)) (loop)))
(define primes (cons 2 (sieve (ints-from 3 2) (ints-from 2 2)))) (for/list ([i 25] [x (in-producer channel-get eof primes)]) x)</lang>
Using generators
Yet another variation of the same algorithm as above, this time using generators.
<lang Racket>#lang racket (require racket/generator) (define (ints-from i d)
(generator () (let loop ([i i]) (yield i) (loop (+ i d)))))
(define (merge g1 g2)
(generator () (let loop ([x1 (g1)] [x2 (g2)]) (cond [(< x1 x2) (yield x1) (loop (g1) x2)] [(> x1 x2) (yield x2) (loop x1 (g2))] [else (yield x1) (loop (g1) (g2))]))))
(define (sieve l non-primes)
(generator () (let loop ([x (l)] [np (non-primes)]) (cond [(< np x) (loop x (non-primes))] [(= np x) (loop (l) (non-primes))] [else (yield x) (set! non-primes (merge (ints-from (* x x) x) non-primes)) (loop (l) (non-primes))]))))
(define (cons x l) (generator () (yield x) (let loop () (yield (l)) (loop)))) (define primes (cons 2 (sieve (ints-from 3 2) (ints-from 2 2)))) (for/list ([i 25] [x (in-producer primes)]) x)</lang>
REXX
no wheel version
The first three REXX versions make use of a sparse stemmed array: [@.].
As the stemmed array gets heavily populated, the number of entries may slow down the REXX interpreter substantially,
depending upon the efficacy of the hashing technique being used for REXX variables (setting/retrieving).
<lang REXX>/*REXX program generates primes via the sieve of Eratosthenes algorithm.*/
parse arg H .; if H== then H=200 /*was the high limit specified? */
w=length(H); @prime=right('prime',20) /*W is used for formatting output*/
@.=. /*assume all numbers are prime. */
- =0 /*number of primes found so far. */
do j=2 for H-1 /*all integers up to H inclusive.*/ if @.j== then iterate /*Composite? Then skip this num.*/ #=#+1 /*bump the prime number counter. */ say @prime right(#,w) " ───► " right(j,w) /*show the prime.*/ do m=j*j to H by j; @.m=; end /*strike all multiples as ¬ prime*/ end /*j*/ /* ─── */ /*stick a fork in it, we're done.*/
say; say right(#,w+length(@prime)+1) 'primes found.'</lang> output when using the input default of: 200
prime 1 ───► 2 prime 2 ───► 3 prime 3 ───► 5 prime 4 ───► 7 prime 5 ───► 11 prime 6 ───► 13 prime 7 ───► 17 prime 8 ───► 19 prime 9 ───► 23 prime 10 ───► 29 prime 11 ───► 31 prime 12 ───► 37 prime 13 ───► 41 prime 14 ───► 43 prime 15 ───► 47 prime 16 ───► 53 prime 17 ───► 59 prime 18 ───► 61 prime 19 ───► 67 prime 20 ───► 71 prime 21 ───► 73 prime 22 ───► 79 prime 23 ───► 83 prime 24 ───► 89 prime 25 ───► 97 prime 26 ───► 101 prime 27 ───► 103 prime 28 ───► 107 prime 29 ───► 109 prime 30 ───► 113 prime 31 ───► 127 prime 32 ───► 131 prime 33 ───► 137 prime 34 ───► 139 prime 35 ───► 149 prime 36 ───► 151 prime 37 ───► 157 prime 38 ───► 163 prime 39 ───► 167 prime 40 ───► 173 prime 41 ───► 179 prime 42 ───► 181 prime 43 ───► 191 prime 44 ───► 193 prime 45 ───► 197 prime 46 ───► 199 46 primes found.
wheel version, optional prime list suppression
This version skips striking the even numbers (as being not prime), 2 is handled as a special case.
Also supported is the suppression of listing the primes if the H (high limit) is negative.
Also added is a final message indicating the number of primes found.
<lang rexx>/*REXX pgm gens primes via a wheeled sieve of Eratosthenes algorithm.*/
parse arg H .; if H== then H=200 /*let the highest # be specified.*/
tell=h>0; H=abs(H); w=length(H) /*neg H suppresses prime listing.*/
if 2<=H & tell then say right(1,w+20)'st prime ───► ' right(2,w)
- = w<=H /*number of primes found so far. */
@.=. /*assume all numbers are prime. */ !=0 /*skips top part of sieve marking*/
do j=3 by 2 for (H-2)%2 /*odd integers up to H inclusive.*/ if @.j== then iterate /*composite? Then skip this num.*/ #=#+1 /*bump the prime number counter. */ if tell then say right(#,w+20)th(#) 'prime ───► ' right(j,w) if ! then iterate /*should the top part be skipped?*/ jj=j*j /*compute the square of J. __ */ if jj>H then !=1 /*indicate skipping if j > √ H.*/ do m=jj to H by j+j; @.m=; end /*strike odd multiples as ¬ prime*/ end /*j*/ /* ─── */
say say right(#,w+20) 'prime's(#) "found." exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────one─liner subroutines────────────────────────────────────────*/ s: if arg(1)==1 then return arg(3); return word(arg(2) 's',1) /*pluralizer.*/ th: procedure; parse arg x; x=abs(x); return word('th st nd rd',1+x//10*(x//100%10\==1)*(x//10<4))</lang> output when using the input default of: 200
1st prime ───► 2 2nd prime ───► 3 3rd prime ───► 5 4th prime ───► 7 5th prime ───► 11 6th prime ───► 13 7th prime ───► 17 8th prime ───► 19 9th prime ───► 23 10th prime ───► 29 11th prime ───► 31 12th prime ───► 37 13th prime ───► 41 14th prime ───► 43 15th prime ───► 47 16th prime ───► 53 17th prime ───► 59 18th prime ───► 61 19th prime ───► 67 20th prime ───► 71 21st prime ───► 73 22nd prime ───► 79 23rd prime ───► 83 24th prime ───► 89 25th prime ───► 97 26th prime ───► 101 27th prime ───► 103 28th prime ───► 107 29th prime ───► 109 30th prime ───► 113 31st prime ───► 127 32nd prime ───► 131 33rd prime ───► 137 34th prime ───► 139 35th prime ───► 149 36th prime ───► 151 37th prime ───► 157 38th prime ───► 163 39th prime ───► 167 40th prime ───► 173 41st prime ───► 179 42nd prime ───► 181 43rd prime ───► 191 44th prime ───► 193 45th prime ───► 197 46th prime ───► 199 46 primes found.
output when using the input of: -1000
168 primes found.
output when using the input of: -10000
1229 primes found.
output when using the input of: -100000
9592 primes found.
output when using the input of: -1000000
78498 primes found.
output when using the input of: -10000000
664579 primes found.
output when using the input of: -100000000
16 +++ @.m= Error 5 running "C:\sieve_of_Eratosthenes.rex", line 16: System resources exhausted
The above (using Regina 3.82 under Windows/XP) shows one of the weaknesses of this implementation of the Sieve of Eratosthenes: it must keep an array of all (if not most) values which is used to strike out composite numbers.
The System resources exhausted error can be postponed by implementing further optimizations (expanding the wheel with low primes).
wheel version
This version skips striking the even numbers (as being not prime).
It also uses a short-circuit test for striking out composites ≤ √target.
<lang rexx>/*REXX pgm gens primes via a wheeled sieve of Eratosthenes algorithm.*/
parse arg H .; if H== then H=200 /*high# can be specified on C.L. */
w=length(H); @prime=right('prime', 20) /*W is used for formatting output*/
if 2<=H then say @prime right(1,w) " ───► " right(2,w)
- = 2<=H /*number of primes found so far. */
@.=. /*assume all numbers are prime. */ !=0 /*skips top part of sieve marking*/
do j=3 by 2 for (H-2)%2 /*odd integers up to H inclusive.*/ if @.j== then iterate /*composite? Then skip this num.*/ #=#+1 /*bump the prime number counter. */ say @prime right(#,w) " ───► " right(j,w) /*show the prime.*/ if ! then iterate /*should the top part be skipped?*/ jj=j*j /*compute the square of J. __ */ if jj>H then !=1 /*indicate skipping if j > √ H.*/ do m=jj to H by j+j; @.m=; end /*strike odd multiples as ¬ prime*/ end /*j*/ /* ─── */
say /*stick a fork in it, we're done.*/
say right(#, w+length(@prime)+1) 'primes found.'</lang>
output is identical to the first (non-wheel) version; program execution is over twice as fast.
The addition of the short-circuit test (variable !) makes it about another 20% faster.
Wheel Version restructured
<lang rexx>/*REXX program generates primes via sieve of Eratosthenes algorithm.
- 21.07.2012 Walter Pachl derived from above Rexx version
- avoid symbols @ and # (not supported by ooRexx)
- avoid negations (think positive)
- /
highest=200 /*define highest number to use. */ is_prime.=1 /*assume all numbers are prime. */ w=length(highest) /*width of the biggest number, */ /* it's used for aligned output.*/ Do j=3 To highest By 2, /*strike multiples of odd ints. */ While j*j<=highest /* up to sqrt(highest) */ If is_prime.j Then Do Do jm=j*3 To highest By j+j /*start with next odd mult. of J */ is_prime.jm=0 /*mark odd mult. of J not prime. */ End End End np=0 /*number of primes shown */ Call tell 2 Do n=3 To highest By 2 /*list all the primes found. */ If is_prime.n Then Do Call tell n End End Exit
tell: Parse Arg prime
np=np+1 Say ' prime number' right(np,w) " --> " right(prime,w) Return</lang>
output is identical to the above versions.
Ruby
eratosthenes starts with nums = [nil, nil, 2, 3, 4, 5, ..., n]
, then marks ( the nil setting ) multiples of 2, 3, 5, 7, ...
there, then returns all non-nil numbers which are the primes.
<lang ruby>def eratosthenes(n)
nums = [nil, nil, *2..n] (2..Math.sqrt(n)).each do |i| (i**2..n).step(i){|m| nums[m] = nil} if nums[i] end nums.compact
end
p eratosthenes(100)</lang>
[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]
With a wheel
eratosthenes2 adds more optimizations, but the code is longer.
- The array
nums
only tracks odd numbers (skips multiples of 2). - The array
nums
holds booleans instead of integers, and every multiple of 3 beginsfalse
. - The outer loop skips multiples of 2 and 3.
- Both inner loops skip multiples of 2 and 3.
<lang ruby>def eratosthenes2(n)
# For odd i, if i is prime, nums[i >> 1] is true. # Set false for all multiples of 3. nums = [true, false, true] * ((n + 5) / 6) nums[0] = false # 1 is not prime. nums[1] = true # 3 is prime.
# Outer loop and both inner loops are skipping multiples of 2 and 3. # Outer loop checks i * i > n, same as i > Math.sqrt(n). i = 5 until (m = i * i) > n if nums[i >> 1] i_times_2 = i << 1 i_times_4 = i << 2 while m <= n nums[m >> 1] = false m += i_times_2 nums[m >> 1] = false m += i_times_4 # When i = 5, skip 45, 75, 105, ... end end i += 2 if nums[i >> 1] m = i * i i_times_2 = i << 1 i_times_4 = i << 2 while m <= n nums[m >> 1] = false m += i_times_4 # When i = 7, skip 63, 105, 147, ... nums[m >> 1] = false m += i_times_2 end end i += 4 end
primes = [2] nums.each_index {|i| primes << (i * 2 + 1) if nums[i]} primes.pop while primes.last > n primes
end
p eratosthenes2(100)</lang>
This simple benchmark compares eratosthenes with eratosthenes2.
<lang ruby>require 'benchmark' Benchmark.bmbm {|x|
x.report("eratosthenes") { eratosthenes(1_000_000) } x.report("eratosthenes2") { eratosthenes2(1_000_000) }
}</lang>
eratosthenes2 runs about 4 times faster than eratosthenes.
With the standard library
MRI 1.9.x implements the sieve of Eratosthenes at file prime.rb, class EratosthensesSeive
(around line 421). This implementation optimizes for space, by packing the booleans into 16-bit integers. It also hardcodes all primes less than 256.
<lang ruby>require 'prime' p Prime::EratosthenesGenerator.new.take_while {|i| i <= 100}</lang>
Run BASIC
<lang runbasic>input "Gimme the limit:"; limit dim flags(limit) for i = 2 to limit
for k = i*i to limit step i flags(k) = 1 next k
if flags(i) = 0 then print i;", "; next i</lang>
Gimme the limit:?100 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,
Rust
Sieve of Eratosthenes - No optimization
<lang rust>fn simple_sieve(limit: usize) -> Vec<usize> {
let mut is_prime = vec![true; limit+1]; is_prime[0] = false; if limit >= 1 { is_prime[1] = false }
for num in 2..limit+1 { if is_prime[num] { let mut multiple = num*num; while multiple <= limit { is_prime[multiple] = false; multiple += num; } } }
is_prime.iter().enumerate() .filter_map(|(pr, &is_pr)| if is_pr {Some(pr)} else {None} ) .collect()
}
fn main() {
println!("{:?}", simple_sieve(100));
}</lang>
- Output:
[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]
Scala
Genuine Eratosthenes sieve
<lang Scala>import scala.annotation.tailrec import scala.collection.parallel.mutable import scala.compat.Platform
object GenuineEratosthenesSieve extends App {
def sieveOfEratosthenes(limit: Int) = { val (primes: mutable.ParSet[Int], sqrtLimit) = (mutable.ParSet.empty ++ (2 to limit), math.sqrt(limit).toInt) @tailrec def prim(candidate: Int): Unit = { if (candidate <= sqrtLimit) { if (primes contains candidate) primes --= candidate * candidate to limit by candidate prim(candidate + 1) } } prim(2) primes } // BitSet toList is shuffled when using the ParSet version. So it has to be sorted before using it as a sequence.
assert(sieveOfEratosthenes(15099480).size == 976729) println(s"Successfully completed without errors. [total ${Platform.currentTime - executionStart} ms]")
}</lang>
- Output:
Successfully completed without errors. [total 39807 ms] Process finished with exit code 0
While concise, the above code is quite slow but a little faster not using the ParSet (take out the '.par' for speed), in which case the sorting ('sorted') is not necessary for an additional small gain in speed; the above code is slow because of all the overhead in processing the bit packed "BitSet" bib-by-bit using complex "higher-order" method calls.
The following '''odds-only''' code is written in a very concise functional style (no mutable state other than the contents of the composites buffer and "higher order functions" for clarity), in this case using a Scala mutable BitSet:
<lang Scala>object SoEwithBitSet {
def makeSoE_PrimesTo(top: Int): Iterator[Int] = { val topNdx = (top - 3) / 2 //odds composite BitSet buffer offset down to 3 val cmpsts = new scala.collection.mutable.BitSet(topNdx + 1) //size includes topNdx @inline def cullPrmCmpsts(prmNdx: Int) = { val prm = prmNdx + prmNdx + 3; cmpsts ++= ((prm * prm - 3) >>> 1 to topNdx by prm) } (0 to (Math.sqrt(top).toInt - 3) / 2).filterNot { cmpsts }.foreach { cullPrmCmpsts } Iterator.single(2) ++ (0 to topNdx).filterNot { cmpsts }.map { pi => pi + pi + 3 } }
}</lang>
In spite of being very concise, it is very much faster than the above code converted to odds-only due to the use of the BitSet instead of the hash table based Set (or ParSet), taking only a few seconds to enumerate the primes to 100 million as compared to the 10's of seconds to count the primes to above 15 million above.
Using tail recursion
The below '''odds-only''' code using a primitive array (bit packed) and tail recursion to avoid some of the enumeration delays due to nested complex "higher order" function calls is almost eight times faster than the above more functional code:
<lang Scala>object SoEwithArray {
def makeSoE_PrimesTo(top: Int) = { import scala.annotation.tailrec val topNdx = (top - 3) / 2 + 1 //odds composite BitSet buffer offset down to 3 plus 1 for overflow val (cmpsts, sqrtLmtNdx) = (new Array[Int]((topNdx >>> 5) + 1), (Math.sqrt(top).toInt - 3) / 2)
@inline def isCmpst(ci: Int): Boolean = (cmpsts(ci >>> 5) & (1 << (ci & 31))) != 0
@inline def setCmpst(ci: Int): Unit = cmpsts(ci >>> 5) |= 1 << (ci & 31)
@tailrec def forCndNdxsFrom(cndNdx: Int): Unit = if (cndNdx <= sqrtLmtNdx) { if (!isCmpst(cndNdx)) { //is prime val p = cndNdx + cndNdx + 3 @tailrec def cullPrmCmpstsFrom(cmpstNdx: Int): Unit = if (cmpstNdx <= topNdx) { setCmpst(cmpstNdx); cullPrmCmpstsFrom(cmpstNdx + p) } cullPrmCmpstsFrom((p * p - 3) >>> 1) } forCndNdxsFrom(cndNdx + 1) }; forCndNdxsFrom(0)
@tailrec def getNxtPrmFrom(cndNdx: Int): Int = if ((cndNdx > topNdx) || !isCmpst(cndNdx)) cndNdx + cndNdx + 3 else getNxtPrmFrom(cndNdx + 1)
Iterator.single(2) ++ Iterator.iterate(3)(p => getNxtPrmFrom(((p + 2) - 3) >>> 1)).takeWhile(_ <= top) }
}</lang>
It can be tested with the following code:
<lang Scala>object Main extends App {
import SoEwithArray._ val top_num = 100000000 val strt = System.nanoTime() val count = makeSoE_PrimesTo(top_num).size val end = System.nanoTime() println(s"Successfully completed without errors. [total ${(end - strt) / 1000000} ms]") println(f"Found $count primes up to $top_num" + ".") println("Using one large mutable Array and tail recursive loops.")
}</lang>
To produce the following output:
- Output:
Successfully completed without errors. [total 661 ms] Found 5761455 primes up to 100000000. Using one large mutable Array and tail recursive loops.
Odds-only page-segmented "infinite" generator version using tail recursion
The above code still uses an amount of memory proportional to the range of the sieve (although bit-packed as 8 values per byte). As well as only sieving odd candidates, the following code uses a fixed range buffer that is about the size of the CPU L2 cache plus only storage for the base primes up to the square root of the range for a large potential saving in RAM memory used as well as greatly reducing memory access times. The use of innermost tail recursive loops for critical loops where the majority of the execution time is spent rather than "higher order" functions from iterators also greatly reduces execution time, with much of the remaining time used just to enumerate the primes output:
<lang Scala>object APFSoEPagedOdds {
import scala.annotation.tailrec private val CACHESZ = 1 << 18 //used cache buffer size private val PGSZ = CACHESZ / 4 //number of int's in cache private val PGBTS = PGSZ * 32 //number of bits in buffer //processing output type is a tuple of low bit (odds) address, // bit range size, and the actual culled page segment buffer. private type Chunk = (Long, Int, Array[Int]) //produces an iteration of all the primes from an iteration of Chunks private def enumChnkPrms(chnks: Stream[Chunk]): Iterator[Long] = { def iterchnk(chnk: Chunk) = { //iterating primes per Chunk val (lw, rng, bf) = chnk @tailrec def nxtpi(i: Int): Int = { //find next prime index not composite if (i < rng && (bf(i >>> 5) & (1 << (i & 31))) != 0) nxtpi(i + 1) else i } Iterator.iterate(nxtpi(0))(i => nxtpi(i + 1)).takeWhile { _ < rng } .map { i => ((lw + i) << 1) + 3 } } //map from index to prime value chnks.toIterator.flatMap { iterchnk } } //culls the composite number bit representations from the bit-packed page buffer //using a given source of a base primes iterator private def cullPg(bsprms: Iterator[Long], lowi: Long, buf: Array[Int]): Unit = { //cull for all base primes until square >= nxt val rng = buf.length * 32; val nxt = lowi + rng @tailrec def cull(bps: Iterator[Long]): Unit = { //given prime then calculate the base start address for prime squared val bp = bps.next(); val s = (bp * bp - 3) / 2 //almost all of the execution time is spent in the following tight loop @tailrec def cullp(j: Int): Unit = { //cull the buffer for given prime if (j < rng) { buf(j >>> 5) |= 1 << (j & 31); cullp(j + bp.toInt) } } if (s < nxt) { //only cull for primes squared less than max //calculate the start address within this given page segment val strt = if (s >= lowi) (s - lowi).toInt else { val b = (lowi - s) % bp if (b == 0) 0 else (bp - b).toInt } cullp(strt); if (bps.hasNext) cull(bps) } } //loop for all primes in range //for the first page, use own bit pattern as a source of base primes //if another source is not given if (lowi <= 0 && bsprms.isEmpty) cull(enumChnkPrms(Stream((0, buf.length << 5, buf)))) //otherwise use the given source of base primes else if (bsprms.hasNext) cull(bsprms) } //makes a chunk given a low address in (odds) bits private def mkChnk(lwi: Long): Chunk = { val rng = PGBTS; val buf = new Array[Int](rng / 32); val bps = if (lwi <= 0) Iterator.empty else enumChnkPrms(basePrms) cullPg(bps, lwi, buf); (lwi, rng, buf) } //new independent source of base primes in a stream of packed-bit arrays //memoized by converting it to a Stream and retaining a reference here private val basePrms: Stream[Chunk] = Stream.iterate(mkChnk(0)) { case (lw, rng, bf) => { mkChnk(lw + rng) } } //produces an infinite iterator over all the chunk results private def itrRslts[R](rsltf: Chunk => R): Iterator[R] = { def mkrslt(lwi: Long) = { //makes tuple of result and next low index val c = mkChnk(lwi); val (_, rng, _) = c; (rsltf(c), lwi + rng) } Iterator.iterate(mkrslt(0)) { case (_, nlwi) => mkrslt(nlwi) } .map { case (rslt, _) => rslt} } //infinite iteration of results //iterates across the "infinite" produced output primes def enumSoEPrimes(): Iterator[Long] = //use itrRsltsMP to produce Chunks iteration Iterator.single(2L) ++ enumChnkPrms(itrRslts { identity }.toStream) //counts the number of remaining primes in a page segment buffer //using a very fast bitcount per integer element //with special treatment for the last page private def countpgto(top: Long, b: Array[Int], nlwp: Long) = { val numbts = b.length * 32; val prng = numbts * 2 @tailrec def cnt(i: Int, c: Int): Int = { //tight int bitcount loop if (i >= b.length) c else cnt (i + 1, c - Integer.bitCount(b(i))) } if (nlwp > top) { //for top in page, calculate int address containing top val bi = ((top - nlwp + prng) >>> 1).toInt val w = bi >>> 5; b(w) |= -2 << (bi & 31) //mark all following as composite for (i <- w + 1 until b.length) b(i) = -1 } //for all int's to end of buffer cnt(0, numbts) } //counting the entire buffer in every case //counts all the primes up to a top value def countSoEPrimesTo(top: Long): Long = { if (top < 2) return 0L else if (top < 3) return 1L //no work necessary //count all Chunks using multi-processing val gen = itrRslts { case (lwi, rng, bf) => val nlwp = (lwi + rng) * 2 + 3; (countpgto(top, bf, nlwp), nlwp) } //a loop to take Chunk's up to including top limit but not past it @tailrec def takeUpto(acc: Long): Long = { val (cnt, nlwp) = gen.next(); val nacc = acc + cnt if (nlwp <= top) takeUpto(nacc) else nacc }; takeUpto(1) }
}</lang>
As the above and all following sieves are "infinite", they all require an extra range limiting condition to produce a finite output, such as the addition of ".takeWhile(_ <= topLimit)" where "topLimit" is the specified range as is done in the following code:
<lang Scala>object MainSoEPagedOdds extends App {
import APFSoEPagedOdds._ countSoEPrimesTo(100) val top = 1000000000 val strt = System.currentTimeMillis() val cnt = enumSoEPrimes().takeWhile { _ <= top }.length
// val cnt = countSoEPrimesTo(top)
val elpsd = System.currentTimeMillis() - strt println(f"Found $cnt primes up to $top in $elpsd milliseconds.")
}</lang>
which outputs the following:
- Output:
Found 50847534 primes up to 1000000000 in 5867 milliseconds.
While the above code is reasonably fast, much of the execution time is consumed by the use of the built-in functions and iterators for concise code, especially in the use of iterators for primes output. To show this, the code includes a "countSoEPrimesTo" function/method that can be uncommented in the above code (commenting out the "takeWhile" line) to produce the following output:
- Output:
Found 50847534 primes up to 1000000000 in 2623 milliseconds.
This shows that it takes somewhat longer to enumerate the primes than it does to actually produce them; this could be improved with a "roll-your-own" enumeration Iterator implementation at considerable increased complexity, but enumeration time will still be a significant portion of the execution time. Further improvements to the code using extreme wheel factorization and multi-processing will make enumeration time an even higher percentage of the total; this is why for large ranges one writes functions/methods similar to "countSoEPrimesTo" to (say) sum the primes, to find the nth prime, etc.
Odds-Only "infinite" generator sieve using Streams and Co-Inductive Streams
The following code uses delayed recursion via Streams to implement the Richard Bird algorithm mentioned in the last part (the Epilogue) of M.O'Neill's paper, which is a true incremental Sieve of Eratosthenes. It is nowhere near as fast as the array based solutions due to the overhead of functionally chasing the merging of the prime multiple streams; this also means that the empirical performance is not according to the usual Sieve of Eratosthenes approximations due to this overhead increasing as the log of the sieved range, but it is much better than the "unfaithful" sieve.
<lang Scala> def birdPrimes() = {
def oddPrimes: Stream[Int] = { def merge(xs: Stream[Int], ys: Stream[Int]): Stream[Int] = { val (x, y) = (xs.head, ys.head) if (y > x) x #:: merge(xs.tail, ys) else if (x > y) y #:: merge(xs, ys.tail) else x #:: merge(xs.tail, ys.tail) } def primeMltpls(p: Int): Stream[Int] = Stream.iterate(p * p)(_ + p + p) def allMltpls(ps: Stream[Int]): Stream[Stream[Int]] = primeMltpls(ps.head) #:: allMltpls(ps.tail) def join(ams: Stream[Stream[Int]]): Stream[Int] = ams.head.head #:: merge(ams.head.tail, join(ams.tail)) def oddPrms(n: Int, composites: Stream[Int]): Stream[Int] = if (n >= composites.head) oddPrms(n + 2, composites.tail) else n #:: oddPrms(n + 2, composites) //following uses a new recursive source of odd base primes 3 #:: oddPrms(5, join(allMltpls(oddPrimes))) } 2 #:: oddPrimes }</lang>
Now this algorithm doesn't really need the memoization and full laziness as offered by Streams, so an implementation and use of a Co-Inductive Stream (CIS) class is sufficient and reduces execution time by almost a factor of two:<lang scala> class CIS[A](val start: A, val continue: () => CIS[A])
def primesBirdCIS: Iterator[Int] = { def merge(xs: CIS[Int], ys: CIS[Int]): CIS[Int] = { val (x, y) = (xs.start, ys.start)
if (y > x) new CIS(x, () => merge(xs.continue(), ys)) else if (x > y) new CIS(y, () => merge(xs, ys.continue())) else new CIS(x, () => merge(xs.continue(), ys.continue())) }
def primeMltpls(p: Int): CIS[Int] = { def nextCull(cull: Int): CIS[Int] = new CIS[Int](cull, () => nextCull(cull + 2 * p)) nextCull(p * p) }
def allMltpls(ps: CIS[Int]): CIS[CIS[Int]] = new CIS[CIS[Int]](primeMltpls(ps.start), () => allMltpls(ps.continue())) def join(ams: CIS[CIS[Int]]): CIS[Int] = { new CIS[Int](ams.start.start, () => merge(ams.start.continue(), join(ams.continue()))) }
def oddPrimes(): CIS[Int] = { def oddPrms(n: Int, composites: CIS[Int]): CIS[Int] = { //"minua" if (n >= composites.start) oddPrms(n + 2, composites.continue()) else new CIS[Int](n, () => oddPrms(n + 2, composites)) } //following uses a new recursive source of odd base primes new CIS(3, () => oddPrms(5, join(allMltpls(oddPrimes())))) }
Iterator.single(2) ++ Iterator.iterate(oddPrimes())(_.continue()).map(_.start) }</lang>
Further gains in performance for these last two implementations can be had by using further wheel factorization and "tree folding/merging" as per this Haskell implementation.
Odds-Only "infinite" generator sieve using a hash table (HashMap)
As per the "unfaithful sieve" article linked above, the incremental "infinite" Sieve of Eratosthenes can be implemented using a hash table instead of a Priority Queue or Map (Binary Heap) as were used in that article. The following implementation postpones the adding of base prime representations to the hash table until necessary to keep the hash table small:
<lang scala> def SoEInc: Iterator[Int] = {
val nextComposites = scala.collection.mutable.HashMap[Int, Int]() def oddPrimes: Iterator[Int] = { val basePrimes = SoEInc basePrimes.next() basePrimes.next() // skip the two and three prime factors @tailrec def makePrime(state: (Int, Int, Int)): (Int, Int, Int) = { val (candidate, nextBasePrime, nextSquare) = state if (candidate >= nextSquare) { val adv = nextBasePrime << 1 nextComposites += ((nextSquare + adv) -> adv) val np = basePrimes.next() makePrime((candidate + 2, np, np * np)) } else if (nextComposites.contains(candidate)) { val adv = nextComposites(candidate) nextComposites -= (candidate) += (Iterator.iterate(candidate + adv)(_ + adv) .dropWhile(nextComposites.contains(_)).next() -> adv) makePrime((candidate + 2, nextBasePrime, nextSquare)) } else (candidate, nextBasePrime, nextSquare) } Iterator.iterate((5, 3, 9)) { case (c, p, q) => makePrime((c + 2, p, q)) } .map { case (p, _, _) => p } } List(2, 3).toIterator ++ oddPrimes }</lang>
The above could be implemented using Streams or Co-Inductive Streams to pass the continuation parameters as passed here in a tuple but there would be no real difference in speed and there is no need to use the implied laziness. As compared to the versions of the Bird (or tree folding) Sieve of Eratosthenes, this has the expected same computational complexity as the array based versions, but is about 20 times slower due to the constant overhead of processing the key value hashing. Memory use is quite low, only being the hash table entries for each of the base prime values less than the square root of the last prime enumerated multiplied by the size of each hash entry (about 12 bytes in this case) plus a "load factor" percentage overhead in hash table size to minimize hash collisions (about twice as large as entries actually used by default on average).
The Scala implementable of a mutable HashMap is slower than the java.util.HashMap one by a factor of almost two, but the Scala version is used here to keep the code more portable (as to CLR). One can also quite easily convert this code to use the immutable Scala HashMap, but the code runs about four times slower due to the required "copy on update" operations for immutable objects.
This algorithm is very responsive to further application of wheel factorization, which can make it run up to about four times faster for the composite number culling operations; however, that is not enough to allow it to catch up to the array based sieves.
Scheme
Tail-recursive solution
<lang scheme>; Tail-recursive solution : (define (sieve n)
(define (aux u v) (let ((p (car v))) (if (> (* p p) n) (let rev-append ((u u) (v v)) (if (null? u) v (rev-append (cdr u) (cons (car u) v)))) (aux (cons p u) (let wheel ((u '()) (v (cdr v)) (a (* p p))) (cond ((null? v) (reverse u)) ((= (car v) a) (wheel u (cdr v) (+ a p))) ((> (car v) a) (wheel u v (+ a p))) (else (wheel (cons (car v) u) (cdr v) a)))))))) (aux '(2) (let range ((v '()) (k (if (odd? n) n (- n 1)))) (if (< k 3) v (range (cons k v) (- k 2))))))
- > (sieve 100)
- (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)
- > (length (sieve 10000000))
- 664579</lang>
Simpler, non-tail-recursive solution
<lang scheme>; Simpler solution, with the penalty that none of 'iota, 'strike or 'sieve is tail-recursive : (define (iota start stop stride)
(if (> start stop) (list) (cons start (iota (+ start stride) stop stride))))
(define (strike lst start stride)
(cond ((null? lst) lst) ((= (car lst) start) (strike (cdr lst) (+ start stride) stride)) ((> (car lst) start) (strike lst (+ start stride) stride)) (else (cons (car lst) (strike (cdr lst) start stride)))))
(define (primes limit)
(let ((stop (sqrt limit))) (define (sieve lst) (let ((p (car lst))) (if (> p stop) lst (cons p (sieve (strike (cdr lst) (* p p) p)))))) (sieve (iota 2 limit 1))))
(display (primes 100)) (newline)</lang> Output: <lang>(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)</lang>
Optimised using an odds-wheel
Optimised using a pre-computed wheel based on 2 (i.e. odds only): <lang scheme>(define (primes-wheel-2 limit)
(let ((stop (sqrt limit))) (define (sieve lst) (let ((p (car lst))) (if (> p stop) lst (cons p (sieve (strike (cdr lst) (* p p) (* 2 p))))))) (cons 2 (sieve (iota 3 limit 2)))))
(display (primes-wheel-2 100)) (newline)</lang> Output: <lang>(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)</lang>
Vector-based
Vector-based (faster), works with RRS: <lang scheme>; initialize v to vector of sequential integers (define (initialize! v)
(define (iter v n) (if (>= n (vector-length v)) (values) (begin (vector-set! v n n) (iter v (+ n 1))))) (iter v 0))
- set every nth element of vector v to 0,
- starting with element m
(define (strike! v m n)
(cond ((>= m (vector-length v)) (values)) (else (begin (vector-set! v m 0) (strike! v (+ m n) n)))))
- lowest non-zero index of vector v >= n
(define (nextprime v n)
(if (zero? (vector-ref v n)) (nextprime v (+ n 1)) (vector-ref v n)))
- remove elements satisfying pred? from list lst
(define (remove pred? lst)
(cond ((null? lst) '()) ((pred? (car lst))(remove pred? (cdr lst))) (else (cons (car lst) (remove pred? (cdr lst))))))
- the sieve itself
(define (sieve n)
(define stop (sqrt n)) (define (iter v p) (cond ((> p stop) v) (else (begin (strike! v (* p p) p) (iter v (nextprime v (+ p 1))))))) (let ((v (make-vector (+ n 1)))) (initialize! v) (vector-set! v 1 0) ; 1 is not a prime (remove zero? (vector->list (iter v 2)))))</lang>
Using SICP-style streams
Using SICP-style head-forced streams. Works with MIT-Scheme – or any other Scheme, if writing out by hand the expansion of the only macro here, s-cons
, with explicit lambda. Common functions:
<lang scheme> ;;;; Stream Implementation
(define (head s) (car s)) (define (tail s) ((cdr s))) (define-syntax s-cons (syntax-rules () ((s-cons h t) (cons h (lambda () t)))))
;;;; Stream Utility Functions (define (from-By x s) (s-cons x (from-By (+ x s) s))) (define (take n s) (cond ((> n 1) (cons (head s) (take (- n 1) (tail s)))) ((= n 1) (list (head s))) ;; don't force it too soon!! (else ()))) ;; so (take 4 (s-map / (from-By 4 -1))) works (define (drop n s) (cond ((> n 0) (drop (- n 1) (tail s))) (else s))) (define (s-map f s) (s-cons (f (head s)) (s-map f (tail s)))) (define (s-diff s1 s2) (let ((h1 (head s1)) (h2 (head s2))) (cond ((< h1 h2) (s-cons h1 (s-diff (tail s1) s2 ))) ((< h2 h1) (s-diff s1 (tail s2))) (else (s-diff (tail s1) (tail s2)))))) (define (s-union s1 s2) (let ((h1 (head s1)) (h2 (head s2))) (cond ((< h1 h2) (s-cons h1 (s-union (tail s1) s2 ))) ((< h2 h1) (s-cons h2 (s-union s1 (tail s2)))) (else (s-cons h1 (s-union (tail s1) (tail s2)))))))</lang>
Tree-folding
Finds composites as a tree of unions of each prime's multiples.
<lang scheme> ;;;; all primes' multiples are removed, merged through a tree of unions
;;;; runs in ~ n^1.15 run time in producing n = 100K .. 1M primes (define (primes-stream) (define (mults p) (from-By (* p p) (* 2 p))) (define (no-mults-From from) (s-diff (from-By from 2) (s-tree-join (s-map mults odd-primes)))) (define odd-primes (s-cons 3 (no-mults-From 5))) ;; inner feedback loop (s-cons 2 (no-mults-From 3))) ;; result stream
;;;; join an ordered stream of streams (here, of primes' multiples) ;;;; into one ordered stream, via an infinite right-deepening tree (define (s-tree-join sts) ;; sts -> s (define (join-With of-Tail sts) ;; sts -> s (s-cons (head (head sts)) (s-union (tail (head sts)) (of-Tail (tail sts))))) (define (pairs sts) ;; sts -> sts (s-cons (join-With head sts) (pairs (tail (tail sts))))) (join-With (lambda (t) (s-tree-join (pairs t))) sts))</lang>
Print 10 last primes of the first thousand primes:
(display (take 10 (drop 990 (primes-stream)))) ; (7841 7853 7867 7873 7877 7879 7883 7901 7907 7919)
The above can be accomplished by the following self contained code which follows the format of the Richard Bird code below with the added "pairs" function integrated into the "mrgmltpls" function:
<lang scheme>(define (treemergePrimes)
(define (mltpls p) (define pm2 (* p 2)) (let nxtmltpl ((cmpst (* p p))) (cons cmpst (lambda () (nxtmltpl (+ cmpst pm2)))))) (define (allmltpls ps) (cons (mltpls (car ps)) (lambda () (allmltpls ((cdr ps)))))) (define (merge xs ys) (let ((x (car xs)) (xt (cdr xs)) (y (car ys)) (yt (cdr ys))) (cond ((< x y) (cons x (lambda () (merge (xt) ys)))) ((> x y) (cons y (lambda () (merge xs (yt))))) (else (cons x (lambda () (merge (xt) (yt)))))))) (define (pairs mltplss) (let ((tl ((cdr mltplss)))) (cons (merge (car mltplss) (car tl)) (lambda () (pairs ((cdr tl))))))) (define (mrgmltpls mltplss) (cons (car (car mltplss)) (lambda () (merge ((cdr (car mltplss))) (mrgmltpls (pairs ((cdr mltplss)))))))) (define (minusstrtat n cmps) (if (< n (car cmps)) (cons n (lambda () (minusstrtat (+ n 2) cmps))) (minusstrtat (+ n 2) ((cdr cmps))))) (define (cmpsts) (mrgmltpls (allmltpls (oddprms)))) ;; internal define's are mutually recursive (define (oddprms) (cons 3 (lambda () (minusstrtat 5 (cmpsts))))) (cons 2 (lambda () (oddprms))))</lang>
It can be tested with the same code as the Richard Bird sieve just by calling "treemergePrimes" instead of "birdPrimes".
Richard Bird's sieve
Archetypal, straightforward code. Uses s-linear-join
, i.e. right fold, which is less efficient and of worse time complexity than the tree-folding above. Does not attempt to conserve space by arranging for the additional inner feedback loop, as is done in the tree-folding variant above.
<lang scheme> (define (primes-stream-ala-Bird)
(define (mults p) (from-By (* p p) (* 2 p))) (define odd-primes ;; primes are (s-cons 3 (s-diff (from-By 5 2) ;; odds, without (s-linear-join (s-map mults odd-primes))))) ;; multiples of primes (s-cons 2 odd-primes))
;;;; join streams using linear structure (define (s-linear-join sts) (s-cons (head (head sts)) (s-union (tail (head sts)) (s-linear-join (tail sts)))))</lang>
Here is a version of the same Richard Bird sieve (mentioned in the link Melissa E. O'Neill article), which is self contained with all the requisite functions wrapped in the overall function; optimized further. It arranges for a separate primes feed for the base primes separate from the output stream, calculated recursively by the recursive call to "oddprms" in forming "cmpsts". It also "fuses" two functions, s-diff
and from-By
, into one, minusstrtat
:
<lang scheme>(define (birdPrimes)
(define (mltpls p) (define pm2 (* p 2)) (let nxtmltpl ((cmpst (* p p))) (cons cmpst (lambda () (nxtmltpl (+ cmpst pm2)))))) (define (allmltpls ps) (cons (mltpls (car ps)) (lambda () (allmltpls ((cdr ps)))))) (define (merge xs ys) (let ((x (car xs)) (xt (cdr xs)) (y (car ys)) (yt (cdr ys))) (cond ((< x y) (cons x (lambda () (merge (xt) ys)))) ((> x y) (cons y (lambda () (merge xs (yt))))) (else (cons x (lambda () (merge (xt) (yt)))))))) (define (mrgmltpls mltplss) (cons (car (car mltplss)) (lambda () (merge ((cdr (car mltplss))) (mrgmltpls ((cdr mltplss))))))) (define (minusstrtat n cmps) (if (< n (car cmps)) (cons n (lambda () (minusstrtat (+ n 2) cmps))) (minusstrtat (+ n 2) ((cdr cmps))))) (define (cmpsts) (mrgmltpls (allmltpls (oddprms)))) ;; internal define's are mutually recursive (define (oddprms) (cons 3 (lambda () (minusstrtat 5 (cmpsts))))) (cons 2 (lambda () (oddprms))))</lang>
It can be tested with the following code:
<lang scheme>(define (nthPrime n)
(let nxtprm ((cnt 0) (ps (birdPrimes))) (if (< cnt n) (nxtprm (+ cnt 1) ((cdr ps))) (car ps))))
(nthPrime 1000000)</lang>
- Output:
15485863
The same code can easily be modified to perform the folded tree case just by writing and integrating a "pairs" function to do the folding along with the merge, which has been done as an alternate tree folding case above.
Scilab
<lang scliab> clear
n=99 sieve=ones(1,n+2) for i=2:n if sieve(i) then for j=i*2:i:n sieve(j)=0 end end end for i=2:n if sieve(i) then disp(i); end end</lang>
Seed7
The program below computes the number of primes between 1 and 10000000: <lang seed7>$ include "seed7_05.s7i";
const func set of integer: eratosthenes (in integer: n) is func
result var set of integer: sieve is EMPTY_SET; local var integer: i is 0; var integer: j is 0; begin sieve := {2 .. n}; for i range 2 to sqrt(n) do if i in sieve then for j range i ** 2 to n step i do excl(sieve, j); end for; end if; end for; end func;
const proc: main is func
begin writeln(card(eratosthenes(10000000))); end func;</lang>
Original source: [1]
SNOBOL4
Using strings instead of arrays, and the square/sqrt optimizations.
<lang SNOBOL4> define('sieve(n)i,j,k,p,str,res') :(sieve_end) sieve i = lt(i,n - 1) i + 1 :f(sv1)
str = str (i + 1) ' ' :(sieve)
sv1 str break(' ') . j span(' ') = :f(return)
sieve = sieve j ' ' sieve = gt(j ^ 2,n) sieve str :s(return) ;* Opt1 res = str (arb ' ') @p ((j ^ 2) ' ') ;* Opt2 str len(p) . res = ;* Opt2
sv2 str break(' ') . k span(' ') = :f(sv3)
res = ne(remdr(k,j),0) res k ' ' :(sv2)
sv3 str = res :(sv1) sieve_end
- # Test and display
output = sieve(100)
end</lang>
Output:
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
Tcl
<lang tcl>package require Tcl 8.5
proc sieve n {
if {$n < 2} {return {}} # create a container to hold the sequence of numbers. # use a dictionary for its speedy access (like an associative array) # and for its insertion order preservation (like a list) set nums [dict create] for {set i 2} {$i <= $n} {incr i} { # the actual value is never used dict set nums $i "" } set primes [list] while {[set nextPrime [lindex [dict keys $nums] 0]] <= sqrt($n)} { dict unset nums $nextPrime for {set i [expr {$nextPrime ** 2}]} {$i <= $n} {incr i $nextPrime} { dict unset nums $i } lappend primes $nextPrime } return [concat $primes [dict keys $nums]]
}
puts [sieve 100] ;# 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</lang>
Summary :/* TI-83 BASIC */
TI-83 BASIC
<lang ti83b>Input "Limit:",N N→Dim(L1) For(I,2,N) 1→L1(I) End For(I,2,SQRT(N)) If L1(I)=1 Then For(J,I*I,N,I) 0→L1(J) End End End For(I,2,N) If L1(I)=1 Then Disp i End End ClrList L1</lang>
UNIX Shell
With array
<lang bash>#!/usr/bin/zsh
function primes() { typeset -a a typeset i j
a[1]="" for (( i = 2; i <= $1; i++ )); do a[$i]=$i done
for (( i = 2; i * i <= $1; i++ )); do if [[ ! -z $a[$i] ]]; then for (( j = i * i; j <= $1; j += i )); do a[$j]="" done fi done print $a }
primes 1000</lang>
<lang bash>function primes { typeset a i=2 j m=$1 # No for (( ... )) loop in pdksh. Use while loop. while (( i <= m )); do a[$i]=$i (( i++ )) done
i=2 while (( j = i * i, j <= m )); do if [[ -n ${a[$i]} ]]; then while (( j <= m )); do unset a[$j] (( j += i )) done fi (( i++ )) done # No print command in bash. Use echo command. echo ${a[*]} }
primes 1000</lang>
Both scripts output a single long line.
2 3 5 7 11 13 17 19 23 ... 971 977 983 991 997
Using variables as fake array
Bourne Shell and Almquist Shell have no arrays. This script works with bash or dash (standard shell in Ubuntu), but uses no specifics of the shells, so it works with plain Bourne Shell as well.
<lang bash>#! /bin/sh
LIMIT=1000
- As a workaround for missing arrays, we use variables p2, p3, ...,
- p$LIMIT, to represent the primes. Values are true or false.
- eval p$i=true # Set value.
- eval \$p$i # Run command: true or false.
- A previous version of this script created a temporary directory and
- used files named 2, 3, ..., $LIMIT to represent the primes. We now use
- variables so that a killed script does not leave extra files. About
- performance, variables are about as slow as files.
i=2 while [ $i -le $LIMIT ] do
eval p$i=true # was touch $i i=`expr $i + 1`
done
i=2 while
j=`expr $i '*' $i` [ $j -le $LIMIT ]
do
if eval \$p$i # was if [ -f $i ] then while [ $j -le $LIMIT ] do eval p$j=false # was rm -f $j j=`expr $j + $i` done fi i=`expr $i + 1`
done
- was echo `ls|sort -n`
echo `i=2
while [ $i -le $LIMIT ]; do eval \\$p$i && echo $i i=\`expr $i + 1\` done`</lang>
With piping
This version works by piping 1s and 0s through sed. The string of 1s and 0s represents the array of primes.
<lang bash># Fill $1 characters with $2. fill () { # This pipeline would begin # head -c $1 /dev/zero | ... # but some systems have no head -c. Use dd. dd if=/dev/zero bs=$1 count=1 2>/dev/null | tr '\0' $2 }
filter () { # Use sed to put an 'x' after each multiple of $1, remove # first 'x', and mark non-primes with '0'. sed -e s/$2/\&x/g -e s/x// -e s/.x/0/g | { if expr $1 '*' $1 '<' $3 > /dev/null; then filter `expr $1 + 1` .$2 $3 else cat fi } }
- Generate a sequence of 1s and 0s indicating primality.
oz () { fill $1 1 | sed s/1/0/ | filter 2 .. $1 }
- Echo prime numbers from 2 to $1.
prime () { # Escape backslash inside backquotes. sed sees one backslash. echo `oz $1 | sed 's/./&\\ /g' | grep -n 1 | sed s/:1//` }
prime 1000</lang>
C Shell
<lang csh># Sieve of Eratosthenes: Echoes all prime numbers through $limit. @ limit = 80
if ( ( $limit * $limit ) / $limit != $limit ) then echo limit is too large, would cause integer overflow. exit 1 endif
- Use $prime[2], $prime[3], ..., $prime[$limit] as array of booleans.
- Initialize values to 1 => yes it is prime.
set prime=( `repeat $limit echo 1` )
- Find and echo prime numbers.
@ i = 2 while ( $i <= $limit ) if ( $prime[$i] ) then echo $i
# For each multiple of i, set 0 => no it is not prime. # Optimization: start at i squared. @ m = $i * $i while ( $m <= $limit ) set prime[$m] = 0 @ m += $i end endif @ i += 1 end</lang>
Ursala
with no optimizations <lang Ursala>#import nat
sieve = ~<{0,1}&& iota; @NttPX ~&r->lx ^/~&rhPlC remainder@rlX~|@r</lang> test program: <lang Ursala>#cast %nL
example = sieve 50</lang>
- Output:
<2,3,5,7,11,13,17,19,23,29,31,37,41,43,47>
Vala
Without any optimizations:
<lang vala>using Gee;
ArrayList<int> primes(int limit){ var sieve = new ArrayList<bool>(); var prime_list = new ArrayList<int>();
for(int i = 0; i <= limit; i++){ sieve.add(true); }
sieve[0] = false; sieve[1] = false;
for (int i = 2; i <= limit/2; i++){ if (sieve[i] != false){ for (int j = 2; i*j <= limit; j++){ sieve[i*j] = false; } } }
for (int i = 0; i < sieve.size; i++){ if (sieve[i] != false){ prime_list.add(i); } }
return prime_list; } // end primes
public static void main(){ var prime_list = primes(50);
foreach(var prime in prime_list) stdout.printf("%s ", prime.to_string());
stdout.printf("\n"); }</lang>{{out}
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47
VBScript
To run in console mode with cscript. <lang vb>
Dim sieve()
If WScript.Arguments.Count>=1 Then n = WScript.Arguments(0) Else n = 99 End If
ReDim sieve(n) For i = 1 To n sieve(i) = True Next For i = 2 To n If sieve(i) Then For j = i * 2 To n Step i sieve(j) = False Next End If Next For i = 2 To n If sieve(i) Then WScript.Echo i Next
</lang>
Vedit macro language
This implementation uses an edit buffer as an array for flags. After the macro has been run, you can see how the primes are located in the array. Primes are marked with 'P' and non-primes with '-'. The first character position represents number 0.
#10 = Get_Num("Enter number to search to: ", STATLINE) Buf_Switch(Buf_Free) // Use edit buffer as flags array Ins_Text("--") // 0 and 1 are not primes Ins_Char('P', COUNT, #10-1) // init rest of the flags to "prime" for (#1 = 2; #1*#1 < #10; #1++) { Goto_Pos(#1) if (Cur_Char=='P') { // this is a prime number for (#2 = #1*#1; #2 <= #10; #2 += #1) { Goto_Pos(#2) Ins_Char('-', OVERWRITE) } } }
Sample output showing numbers in range 0 to 599.
--PP-P-P---P-P---P-P---P-----P-P-----P---P-P---P-----P-----P -P-----P---P-P-----P---P-----P-------P---P-P---P-P---P------ -------P---P-----P-P---------P-P-----P-----P---P-----P-----P -P---------P-P---P-P-----------P-----------P---P-P---P-----P -P---------P-----P-----P-----P-P-----P---P-P---------P------ -------P---P-P---P-------------P-----P---------P-P---P-----P -------P-----P-----P---P-----P-------P---P-------P---------P -P---------P-P-----P---P-----P-------P---P-P---P-----------P -------P---P-------P---P-----P-----------P-P---------------- -P-----P---------P-----P-----P-P-----P---------P-----P-----P
Visual Basic
Works with: VB6 <lang vb>Sub Eratost()
Dim sieve() As Boolean Dim n As Integer, i As Integer, j As Integer n = InputBox("limit:", n) ReDim sieve(n) For i = 1 To n sieve(i) = True Next i For i = 2 To n If sieve(i) Then For j = i * 2 To n Step i sieve(j) = False Next j End If Next i For i = 2 To n If sieve(i) Then Debug.Print i Next i
End Sub 'Eratost</lang>
Visual Basic .NET
<lang vbnet>Dim n As Integer, k As Integer, limit As Integer Console.WriteLine("Enter number to search to: ") limit = Console.ReadLine Dim flags(limit) As Integer For n = 2 To Math.Sqrt(limit)
If flags(n) = 0 Then For k = n * n To limit Step n flags(k) = 1 Next k End If
Next n
' Display the primes For n = 2 To limit
If flags(n) = 0 Then Console.WriteLine(n) End If
Next n</lang>
Vorpal
<lang vorpal>self.print_primes = method(m){
p = new() p.fill(0, m, 1, true)
count = 0 i = 2 while(i < m){ if(p[i] == true){ p.fill(i+i, m, i, false) count = count + 1 } i = i + 1 } ('primes: ' + count + ' in ' + m).print() for(i = 2, i < m, i = i + 1){ if(p[i] == true){ ( + i + ', ').put() } } .print()
}
self.print_primes(100)</lang>
- Result:
primes: 25 in 100 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,
XPL0
<lang XPL0>include c:\cxpl\codes; \intrinsic 'code' declarations int Size, Prime, I, Kill; char Flag; [Size:= IntIn(0); Flag:= Reserve(Size+1); for I:= 2 to Size do Flag(I):= true; for I:= 2 to Size do
if Flag(I) then \found a prime [Prime:= I; IntOut(0, Prime); CrLf(0); Kill:= Prime + Prime; \first multiple to kill while Kill <= Size do [Flag(Kill):= false; \zero a non-prime Kill:= Kill + Prime; \next multiple ]; ];
]</lang>
- Example output:
202 3 5 7 11 13 17
19
zkl
<lang zkl>fcn sieve(limit){
if (limit<2) return(T); composite:=(0).pump(limit+1,Data,1); // bucket of bytes set to 1 (prime) (2).filter(limit.toFloat().sqrt()+1,T(composite.get, // if prime, zero multiples 'wrap(n){ [n*n .. limit,n].pump(Void,composite.set.fp1(0)) }, //composite[n*p]=0 False)); // turn filter into a no-result loop (2).filter(limit-1,composite.get); // bytes still 1 are prime
} sieve(53).println();</lang> The filter method, when given multiple filters, acts like a conditional and. Here, the first filter checks the table for a prime, if so, the second filter does some side effects and the third filter ensures that no items make it through the filter (False(<anything>)-->False) so that the filter returns an empty list, minimizing garbage.
- Output:
L(2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53)
- Programming Tasks
- Prime Numbers
- Clarified and Needing Review
- 360 Assembly
- 68000 Assembly
- ABAP
- ACL2
- Ada
- ALGOL 60
- ALGOL 68
- ALGOL W
- AutoHotkey
- AutoIt
- AWK
- BASIC
- Applesoft BASIC
- Locomotive Basic
- ZX Spectrum Basic
- BBC BASIC
- Batch File
- Befunge
- Bracmat
- C
- C++
- C sharp
- Chapel
- Clojure
- CMake
- COBOL
- Common Lisp
- D
- Dart
- Delphi
- DWScript
- Dylan
- E
- EC
- EC examples needing attention
- Examples needing attention
- Eiffel
- Elixir
- Elixir examples needing attention
- Emacs Lisp
- Erlang
- ERRE
- Euphoria
- F Sharp
- Forth
- Fortran
- GAP
- GLBasic
- Go
- Groovy
- GW-BASIC
- Haskell
- HicEst
- Icon
- Unicon
- J
- Java
- JavaScript
- Jq
- Liberty BASIC
- Logo
- Lua
- Lucid
- Lucid examples needing attention
- M4
- Mathematica
- MATLAB
- Maxima
- MAXScript
- Modula-3
- Modula-3 examples needing attention
- MUMPS
- NetRexx
- Nial
- Nial examples needing attention
- Nim
- Niue
- Niue examples needing attention
- Oberon-2
- OCaml
- Oforth
- Oz
- PARI/GP
- Pascal
- Perl
- Perl 6
- PHP
- PicoLisp
- PL/I
- Pop11
- PowerShell
- Processing
- Prolog
- PureBasic
- Python
- Numpy
- R
- Racket
- REXX
- Ruby
- Run BASIC
- Rust
- Scala
- Scheme
- Scilab
- Seed7
- SNOBOL4
- Tcl
- TI-83 BASIC
- UNIX Shell
- C Shell
- Ursala
- Ursala examples needing attention
- Vala
- Gee
- VBScript
- Vedit macro language
- Visual Basic
- Visual Basic .NET
- Vorpal
- XPL0
- Zkl