Sieve of Eratosthenes: Difference between revisions

{{trans|Python}}

V is_prime = [0B]*2 [+] [1B]*(limit - 1)

L(n) 0 .< Int(limit ^ 0.5 + 1.5)

R enumerate(is_prime).filter((i, prime) -> prime).map((i, prime) -> i)

{{out}}

=={{header|360 Assembly}}==

For maximum compatibility, this program uses only the basic instruction set.

ERATOST CSECT

USING ERATOST,R12

CRIBLE DC 100000X'01'

YREGS

{{out}}

<pre style="height:20ex">

=={{header|6502 Assembly}}==

If this subroutine is called with the value of <i>n</i> in the accumulator, it will store an array of the primes less than <i>n</i> beginning at address 1000 hex and return the number of primes it has found in the accumulator.

LDA #$00

LDX #$00

JMP COPY

COPIED: TYA ; how many found

=={{header|68000 Assembly}}==

Some of the macro code is derived from the examples included with EASy68K.

See 68000 "100 Doors" listing for additional information.

* Title : BitSieve

* Written by : G. A. Tippery

Summary2: DC.B ' prime numbers found.',CR,LF,$00

=={{header|8086 Assembly}}==

cpu 8086

bits 16

numbuf: db 13,10,'$'

section .bss

{{out}}

=={{header|8th}}==

with: n

: sieve SED: n -- a

dup make-sieve swap sieve>a ;

{{Out}}

<pre>

=={{header|AArch64 Assembly}}==

{{works with|as|Raspberry Pi 3B version Buster 64 bits}}

/* ARM assembly AARCH64 Raspberry PI 3B */

/* program cribleEras64.s */

/* for this file see task include a file in language AArch64 assembly */

.include "../includeARM64.inc"

<pre>

Prime : 2

=={{header|ABAP}}==

PARAMETERS: p_limit TYPE i OBLIGATORY DEFAULT 100.

ENDIF.

ENDLOOP.

=={{header|ACL2}}==

(declare (xargs :measure (nfix (- n i))))

(if (zp (- n i))

(defun sieve (limit)

=={{header|Action!}}==

PROC Main()

i==+1

OD

{{out}}

[https://gitlab.com/amarok8bit/action-rosetta-code/-/raw/master/images/Sieve_of_Eratosthenes.png Screenshot from Atari 8-bit computer]

=={{header|ActionScript}}==

Works with ActionScript 3.0 (this is utilizing the actions panel, not a separated class file)

{

var primes:Array = new Array();

}

var e:Array = eratosthenes(1000);

Output:

{{out}}

=={{header|Ada}}==

procedure Eratos is

end if;

end loop;

{{out}}

=={{header|Agda}}==

-- imports

open import Data.Nat as ℕ using (ℕ; suc; zero; _+_; _∸_)

remove _ ys zero = nothing ∷ ys i

remove y ys (suc j) = y ∷ ys j

=={{header|Agena}}==

Tested with Agena 2.9.5 Win32

# generate and return a sequence containing the primes up to sieveSize

# test the sieve proc

{{out}}

<pre>

'''Works with:''' ALGOL 60 for OS/360

'INTEGER' 'ARRAY' CANDIDATES(/0..1000/);

'INTEGER' I,J,K;

'END'

'END'

=={{header|ALGOL 68}}==

PROC eratosthenes = (INT n)[]BOOL:

(

);

{{out}}

<pre>

=={{header|ALGOL W}}==

=== Standard, non-optimised sieve ===

% implements the sieve of Eratosthenes %

end

{{out}}

<pre>

Alternative version that only stores odd numbers greater than 1 in the sieve.

% implements the sieve of Eratosthenes %

% only odd numbers appear in the sieve, which starts at 3 %

for i := 1 until arrayMax do if s( i ) then writeon( ( i * 2 ) + 1 );

end

{{out}}

Same as the standard version.

=={{header|ALGOL-M}}==

BEGIN

END

{{out}}

<pre>

=== Non-Optimized Version ===

b←⍵⍴1

b[⍳2⌊⍵]←0

}

The required list of prime divisors obtains by recursion (<tt>{⍵/⍳⍴⍵}∇⌈⍵*0.5</tt>).

=== Optimized Version ===

b←⍵⍴{∧⌿↑(×/⍵)⍴¨~⍵↑¨1}2 3 5

b[⍳6⌊⍵]←(6⌊⍵)⍴0 0 1 1 0 1

}

The optimizations are as follows:

=== Examples ===

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

+/∘sieve¨ 10*⍳10

The last expression computes the number of primes < 1e0 1e1 ... 1e9.

=={{header|AppleScript}}==

script o

property numberList : {missing value}

end sieveOfEratosthenes

{{out}}

=={{header|ARM Assembly}}==

{{works with|as|Raspberry Pi}}

/* ARM assembly Raspberry PI */

/***************************************************/

.include "../affichage.inc"

<pre>

Prime : 2

=={{header|Arturo}}==

composites: array.of: inc upto false

loop 2..to :integer sqrt upto 'x [

]

{{out}}

=={{header|AutoHotkey}}==

{{AutoHotkey case}}Source: [http://www.autohotkey.com/forum/topic44657.html AutoHotkey forum] by Laszlo

Sieve(n) { ; Sieve of Eratosthenes => string of 0|1 chars, 1 at position k: k is prime

}

Return S

===Alternative Version===

arr := []

loop % n-1

}

return Arr

Arr := Sieve_of_Eratosthenes(n)

loop, % n-1

output .= (Arr[A_Index] ? A_Index : ".") . (!Mod(A_Index, 10) ? "`n" : "`t")

MsgBox % output

{{out}}

<pre>. 2 3 . 5 . 7 . . .

=={{header|AutoIt}}==

$M = InputBox("Integer", "Enter biggest Integer")

Global $a[$M], $r[$M], $c = 1

$r[0] = $c - 1

ReDim $r[$c]

=={{header|AWK}}==

input from commandline as well as stdin,

and input is checked for valid numbers:

# usage: gawk -v n=101 -f sieve.awk

END { print "Bye!" }

Here is an alternate version that uses an associative array to record composites with a prime dividing it. It can be considered a slow version, as it does not cross out composites until needed. This version assumes enough memory to hold all primes up to ULIMIT. It prints out noncomposites greater than 1.

BEGIN { ULIMIT=100

else print ( S[(n+n)] = n )

}

==Bash==

{{works with|FreeBASIC}}

{{works with|RapidQ}}

INPUT "Enter number to search to: "; limit

FOR n = 2 TO limit

IF flags(n) = 0 THEN PRINT n; ", ";

==={{header|Applesoft BASIC}}===

20 DIM FLAGS(LIMIT)

30 FOR N = 2 TO SQR (LIMIT)

100 FOR N = 2 TO LIMIT

110 IF FLAGS(N) = 0 THEN PRINT N;", ";

==={{header|Atari BASIC}}===

{{trans|Commodore BASIC}}

Auto-initialization of arrays is not reliable, so we have to do our own. Also, PRINTing with commas doesn't quite format as nicely as one might hope, so we do a little extra work to keep the columns lined up.

110 PRINT "LIMIT";:INPUT LI

120 DIM N(LI):FOR I=0 TO LI:N(I)=1:NEXT I

250 IF C=3 THEN PRINT:C=0

260 NEXT I

{{Out}}

<pre> Ready

==={{header|Commodore BASIC}}===

Since C= BASIC initializes arrays to all zeroes automatically, we avoid needing our own initialization loop by simply letting 0 mean prime and using 1 for composite.

110 INPUT "LIMIT";LI

120 DIM N(LI)

230 NEXT I

240 PRINT

{{Out}}

<pre>

==={{header|IS-BASIC}}===

110 LET LIMIT=100

120 NUMERIC T(1 TO LIMIT)

230 FOR I=2 TO LIMIT ! Display the primes

240 IF T(I)=0 THEN PRINT I;

==={{header|Locomotive Basic}}===

20 INPUT "Limit";limit

30 DIM f(limit)

100 FOR n=2 TO limit

110 IF f(n)=0 THEN PRINT n;",";

==={{header|MSX Basic}}===

6 REM Translated from Rosetta's ZX Spectrum implementation

10 INPUT "Enter number to search to: ";l

100 FOR n=2 TO l

110 IF p(n)=0 THEN PRINT n;", ";

==={{header|Sinclair ZX81 BASIC}}===

A note on <code>FAST</code> and <code>SLOW</code>: under normal circumstances the CPU spends about 3/4 of its time driving the display and only 1/4 doing everything else. Entering <code>FAST</code> mode blanks the screen (which we do not want to update anyway), resulting in substantially improved performance; we then return to <code>SLOW</code> mode when we have something to print out.

20 FAST

30 DIM N(L)

110 FOR I=2 TO L

120 IF NOT N(I) THEN PRINT I;" ";

==={{header|ZX Spectrum Basic}}===

20 DIM p(l)

30 FOR n=2 TO SQR l

100 FOR n=2 TO l

110 IF p(n)=0 THEN PRINT n;", ";

==={{header|QBasic}}===

{{works with|QBasic|1.1}}

{{works with|QuickBasic|4.5}}

DIM flags(limit)

FOR n = 1 TO limit

IF flags(n) THEN PRINT USING "####"; n;

{{out}}

<pre>Prime numbers less than or equal to 120 are:

==={{header|BASIC256}}===

limit = 120

for n = 1 to limit

if flags[n] then print rjust(n,4);

{{out}}

<pre>Prime numbers less than or equal to 120 are:

==={{header|True BASIC}}===

{{trans|QBasic}}

DIM flags(0)

MAT redim flags(limit)

IF flags(n)<>0 then PRINT using "####": n;

NEXT n

{{out}}

<pre>Same as QBasic entry.</pre>

{{trans|ZX Spectrum Basic}}

Sets h$ to 1 for higher multiples of 2 via <code>FILL$</code>, later on sets <code>STEP</code> to 2n; replaces Floating Pt array p(z) with string variable h$(z) to sieve out all primes < z=441 (<code>l</code>=21) in under 1K, so that h$ is fillable to its maximum (32766), even on a 48K ZX Spectrum if translated back.

10 INPUT "Enter Stopping Pt for squared factors: ";z

15 LET l=SQR(z)

90 REM Display the primes

100 FOR n=2 TO z: IF h$(n)=0: PRINT n;", ";

====2i wheel emulation of Sinclair ZX81 BASIC====

10 INPUT Z

15 LET L=SQR(Z)

120 IF H$(I)="0" THEN PRINT I!

130 NEXT I

====2i wheel emulation of Sinclair ZX80 BASIC====

. . . with 2:1 compression (of 16-bit integer variables on ZX80s) such that it obviates having to account for any multiple of 2; one has to input odd upper limits on factors to be squared, <code>L</code> (=21 at most on 1K ZX80s for all primes till 439).

10 INPUT L

15 LET Z=(L+1)*(L- 1)/2

120 IF NOT H(I) THEN PRINT 0^I+1+I*2!

130 NEXT I

====Sieve of Sundaram====

duplicate entries are omitted. Thus, a slightly transformed Sieve of Sundaram is what Eratosthenes' Sieve becomes upon applying all

optimisations incorporated into the prior entries for QL SuperBASIC, except for any equivalent to line 50 in them.

10 INPUT O

15 LET Z=2*O*O+O*2

120 IF NOT H(I) THEN PRINT 0^I+1+I*2!

130 NEXT I

====Eulerian optimisation====

While slower than the optimised Sieve of Eratosthenes before it, the Sieve of Sundaram above has a compatible compression scheme that's more convenient than the conventional one used beforehand. It is therefore applied below along with Euler's alternative optimisation in a reversed implementation that lacks backward-compatibility to ZX80 BASIC. This program is designed around features & limitations of the QL, yet can be rewritten more efficiently for 1K ZX80s, as they allow integer variables to be parameters of <code>FOR</code> statements (& as their 1K of static RAM is equivalent to L1 cache, even in <code>FAST</code> mode). That's left as an exercise for ZX80 enthusiasts, who for o%=14 should be able to generate all primes < 841, i.e. 3 orders of (base 2) magnitude above the limit for the program listed under Sinclair ZX81 BASIC. In QL SuperBASIC, o% may at most be 127--generating all primes < 65,025 (almost 2x the upper limit for indices & integer variables used to calculate them ~2x faster than for floating point as used in line 30, after which the integer code mimics an assembly algorithm for the QL's 68008.)

10 INPUT "Enter highest value of diagonal index q%: ";o%

15 LET z%=o%*(2+o%*2) : h$=FILL$(" ",z%+o%) : q%=1 : q=q% : m=z% DIV (2*q%+1)

110 GOTO 100

127 DEF FN N%(i)=0^i+1+i*2

=={{header|Batch File}}==

@echo off

setlocal ENABLEDELAYEDEXPANSION

if !crible.%%i! EQU 1 echo %%i

)

{{Out}}

<pre style="height:20ex">limit: 100

=={{header|BBC BASIC}}==

DIM sieve% limit%

FOR I% = 1 TO limit%

IF sieve%?I% = 0 PRINT I%;

=={{header|BCPL}}==

manifest $( LIMIT = 1000 $)

$)

wrch('*N')

{{out}}

<pre> 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71

=== Odds-only bit packed array version (64 bit) ===

This sieve also uses an iterator structure to enumerate the primes in the sieve. It's inspired by the golang bit packed sieve that returns a closure as an iterator. However, BCPL does not support closures, so the code uses an iterator object.

GET "libhdr"

RESULTIS 0

}

{{Out}}

<pre>

A more efficient sieve (primes below one billion in under a minute) is provided as <code>PrimesTo</code> in bqn-libs [https://github.com/mlochbaum/bqn-libs/blob/master/primes.bqn primes.bqn].

𝕩≤2 ? ↕0 ; # No primes below 2

p ← 𝕊⌈√n←𝕩 # Initial primes by recursion

E ← {↕∘⌈⌾((𝕩×𝕩+⊢)⁼)n} # Multiples of 𝕩 under n, starting at 𝕩×𝕩

/ b E⊸{0¨⌾(𝕨⊸⊏)𝕩}´ p # Cross them out

{{out}}

⟨ 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 ⟩

≠∘Primes¨ 10⋆↕7 # Number of primes below 1e0, 1e1, ... 1e6

=={{header|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.

= n j i

. !arg:?n

)

& eratosthenes$100

{{out}}

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

=={{header|C}}==

#include <math.h>

}

return sieve;

'''Another example:'''

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.

#include <malloc.h>

void sieve(int *, int);

}

printf("\n\n");

j:2

Before a[2*2]: 1

i:9

i:10

=={{header|C sharp|C#}}==

{{works with|C sharp|C#|2.0+}}

using System.Collections;

using System.Collections.Generic;

}

}

===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](http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf) in Haskell:

using System.Collections;

using System.Collections.Generic;

}

IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); }

===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:

using System.Collections;

using System.Collections.Generic;

}

IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); }

The above code runs over ten times faster than the original Richard Bird algorithm.

First, an implementation of a Min Heap Priority Queue is provided as extracted from the entry at [http://rosettacode.org/wiki/Priority_queue#C.23 RosettaCode], with only the necessary methods duplicated here:

using KeyT = System.UInt32;

using System;

public static MinHeapPQ<V> replaceMin(KeyT k, V v, MinHeapPQ<V> pq) {

pq.rplcmin(k, v); return pq; }

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:

using System.Collections;

using System.Collections.Generic;

public IEnumerator<PrimeT> GetEnumerator() { return nmrtr(); }

IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); }

The above code is at least about 2.5 times faster than the Tree Folding version.

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:

using System.Collections;

using System.Collections.Generic;

public IEnumerator<PrimeT> GetEnumerator() { return nmrtr(); }

IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); }

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".

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:

using System.Collections;

using System.Collections.Generic;

public IEnumerator<PrimeT> GetEnumerator() { return nmrtr(); }

IEnumerator IEnumerable.GetEnumerator() { return (IEnumerator)GetEnumerator(); }

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.

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

Console.WriteLine(PrimesXXX().ElementAt(1000000 - 1)); // zero based indexing...

To produce the following output for all tested versions (although some are considerably faster than others):

This implementation follows the standard library pattern of [http://en.cppreference.com/w/cpp/algorithm/iota std::iota]. The start and end iterators are provided for the container. The destination container is used for marking primes and then filled with the primes which are less than the container size. This method requires no memory allocation inside the function.

#include <algorithm>

#include <vector>

std::cout << std::endl;

return 1;

=== Boost ===

template<class UnaryFunction>

void primesupto(int limit, UnaryFunction yield)

if (is_prime[n])

yield(n);

Full program:

$ g++ -I/path/to/boost sieve.cpp -o sieve && sieve 10000000

*/

cout << "limit sum pi(n)\n"

<< limit << " " << sum << " " << count << endl;

=={{header|Chapel}}==

{{incorrect|Chapel|Doesn't compile since at least Chapel version 1.20 to 1.24.1.}}

This solution uses nested iterators to create new wheels at run time:

iter sieve(prime):int {

yield candidate;

}

for p in sieve(2) {

if p > N {

}

write(" ", p);

===Alternate Conventional Bit-Packed Implementation===

compile with the `--fast` option

use BitOps;

for p in 2 .. n do if cmpsts[p >> 3] & (1: uint(8) << (p & 7)) == 0 then yield p;

{{out}}

===Alternate Odds-Only Bit-Packed Implementation===

use BitOps;

if cmpsts[i >> 3] & (1: uint(8) << (i & 7)) == 0 then yield i + i + 3;

{{out}}

{{works with|Chapel|1.25.1}}

config const limit = 100000000;

write("Found ", count, " primes up to ", limit);

writeln(" in ", timer.elapsed(TimeUnits.milliseconds), " milliseconds.");

{{out}}

Compile with the `--fast` compiler command line option

config const limit = 100000000;

write("Found ", count, " primes up to ", limit);

writeln(" in ", timer.elapsed(TimeUnits.milliseconds), " milliseconds.");

{{out}}

Compile with the `--fast` compiler command line option

use Map;

write("Found ", count, " primes up to ", limit);

writeln(" in ", timer.elapsed(TimeUnits.milliseconds), " milliseconds.");

This works in about exactly the same time as the last previous code, but doesn't require special custom adaptations of the associative array so that the standard library Map can be used.

{{trans|Nim}} [https://rosettacode.org/wiki/Sieve_of_Eratosthenes#Nim_Unbounded_Versions code link]

{{works with|Chapel|1.22|- compile with the --fast compiler command line flag for full optimization}}

type Prime = uint(32);

for p in primes() { if p > limit then break; cnt += 1; }

timer.stop(); write("\nFound ", cnt, " primes up to ", limit);

{{out}}

<pre>The first 25 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

To take advantage of the features that make Chapel shine, we need to use it to do some parallel computations, and to efficiently do that for the Sieve of Eratosthenes, we need to divide the work into page segments where we can assign each largish segment to a separate thread/task; this also improves the speed due to better cache associativity with most memory accesses to values that are already in the cache(s). Once we have divided the work, Chapel offers lots of means to implement the parallelism but to be a true Sieve of Eratosthenes, we need to have the ability to output the results in order; with many of the convenience mechanisms not doing that, the best/simplest option is likely task parallelism with the output results assigned to an rotary indexed array containing the `sync` results. It turns out that, although the Chapel compiler can sometimes optimize the code so the overhead of creating tasks is not onerous, for this case where the actual tasks are somewhat complex, the compiler can't recognize that an automatically generated thread pool(s) are required so we need to generate the thread pool(s) manually. The code that implements the multi-threading of page segments using thread pools is as follows:

{{works with|Chapel|1.24.1|- compile with the --fast compiler command line flag for full optimization}}

type Prime = uint(64);

write("Found ", count, " primes up to ", LIMIT);

{{out}}

<pre>The first 25 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

=={{header|Clojure}}==

(remove (set (mapcat #(range (* % %) n %)

(range 2 (Math/sqrt n))))

The above is **not strictly a Sieve of Eratosthenes** as the composite culling ranges (in the ''mapcat'') include all of the multiples of all of the numbers and not just the multiples of primes. When tested with <code clojure>(println (time (count (primes< 1000000))))</code>, it takes about 5.5 seconds just to find the number of primes up to a million, partly because of the extra work due to the use of the non-primes, and partly because of the constant enumeration using sequences with multiple levels of function calls. Although very short, this code is likely only useful up to about this range of a million.

It may be written using the ''into #{}'' function to run slightly faster due to the ''set'' function being concerned with only distinct elements whereas the ''into #{}'' only does the conjunction, and even at that doesn't do that much as it does the conjunction to an empty sequence, the code as follows:

(remove (into #{}

(mapcat #(range (* % %) n %)

(range 2 (Math/sqrt n))))

The above code is slightly faster for the reasons given, but is still not strictly a Sieve of Eratosthenes due to sieving by all numbers and not just by the base primes.

The following code also uses the ''into #{}'' transducer but has been slightly wheel-factorized to sieve odds-only:

(if (< n 2) ()

(cons 2 (remove (into #{}

(mapcat #(range (* % %) n %)

(range 3 (Math/sqrt n) 2)))

The above code is a little over twice as fast as the non-odds-only due to the reduced number of operations. It still isn't strictly a Sieve of Eratosthenes as it sieves by all odd base numbers and not only by the base primes.

The following code calculates primes up to and including ''n'' using a mutable boolean array but otherwise entirely functional code; it is tens (to a hundred) times faster than the purely functional codes due to the use of mutability in the boolean array:

"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes"

[n]

(do (dorun (map #(cullp %) (filter #(not (aget cmpsts %))

(range 2 (inc root)))))

'''Alternative implementation using Clojure's side-effect oriented list comprehension.'''

"Returns a lazy sequence of prime numbers less than lim"

[lim]

(doseq [j (range (* i i) lim i)]

(aset refs j false)))

'''Alternative implementation using Clojure's side-effect oriented list comprehension. Odds only.'''

"Returns a lazy sequence of prime numbers less than lim"

[lim]

(doseq [j (range (* (+ i i) (inc i)) max-i (+ i i 1))]

(aset refs j false)))

This implemantation is about twice as fast as the previous one and uses only half the memory.

'''Alternative very slow entirely functional implementation using lazy sequences'''

"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes"

[n]

(cons p (lazy-seq (nxtprm (-> (range (* p p) (inc n) p)

set (remove cs) rest)))))))]

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.

Here is an immutable boolean vector based non-lazy sequence version other than for the lazy sequence operations to output the result:

"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes"

[max-prime]

(assoc 1 false)

(sieve 2))

The above code is still quite slow due to the cost of the immutable copy-on-modify operations.

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:

(defn primes-to

(recur (inc i)))))))))]

(if (< n 2) nil

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.

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

"Computes lazy sequence of prime numbers up to a given number using sieve of Eratosthenes"

[n]

(toString [this] (if (= cnt tcnt) "()"

(.toString (seq (map identity this))))))

'(time (count (primes-tox 10000000)))' takes about 40 milliseconds (compiled) to produce 664579.

'''A Clojure version of Richard Bird's Sieve using Lazy Sequences (sieves odds only)'''

"Computes the unbounded sequence of primes using a Sieve of Eratosthenes algorithm by Richard Bird."

[]

(do (def oddprms (cons 3 (lazy-seq (let [cmpsts (-> oddprms (allmtpls) (mrgmltpls))]

(minusStrtAt 5 cmpsts)))))

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.

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:

"Computes the unbounded sequence of primes using a Sieve of Eratosthenes algorithm modified from Bird."

[]

(do (def oddprms (cons 3 (lazy-seq (let [cmpsts (-> oddprms (allmtpls) (mrgmltpls))]

(minusStrtAt 5 cmpsts)))))

The above code is still slower than it should be due to the slowness of Clojure's sequence operations.

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:

clojure.lang.ISeq

(first [_] v)

(do (def oddprms (->CIS 3 (fn [] (let [cmpsts (-> oddprms (allmtpls) (mrgmltpls))]

(minusStrtAt 5 cmpsts)))))

'(time (count (take-while #(<= (long %) 10000000) (primes-treeFoldingx))))' takes about 3.4 seconds for a range of 10 million.

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:

"Infinite sequence of primes using an incremental Sieve or Eratosthenes with a Hashmap"

[]

(cons c (lazy-seq (nxtoddprm (+ c 2) q bsprms cmpsts))))))]

(do (def baseoddprms (cons 3 (lazy-seq (nxtoddprm 5 9 baseoddprms {}))))

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.

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:

Object

(toString [_] (str "<" k "," v ">")))

(if (<= kl kr)

(recur l #(cont (->PQNode kvl % r)))

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:

"Infinite sequence of primes using an incremental Sieve or Eratosthenes with a Priority Queue"

[]

(cons c (lazy-seq (nxtoddprm (+ c 2) q bsprms cmpsts)))))))]

(do (def baseoddprms (cons 3 (lazy-seq (nxtoddprm 5 9 baseoddprms (empty-pq)))))

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.

To show that Clojure does not need to be particularly slow, the following version runs about twice as fast as the non-segmented unbounded array based version above (extremely fast compared to the non-array based versions) and only a little slower than other equivalent versions running on virtual machines: C# or F# on DotNet or Java and Scala on the JVM:

(def PGSZ (bit-shift-left 1 14)) ;; size of CPU cache

(next pgseq)

(+ cnt (count-pg PGBTS pg))))))]

The above code runs just as fast as other virtual machine languages when run on a 64-bit JVM; however, when run on a 32-bit JVM it runs almost five times slower. This is likely due to Clojure only using 64-bit integers for integer operations and these operations getting JIT compiled to use library functions to simulate those operations using combined 32-bit operations under a 32-bit JVM whereas direct CPU operations can be used on a 64-bit JVM

=={{header|CLU}}==

eratosthenes = proc (n: int) returns (array[bool])

prime: array[bool] := array[bool]$fill(1, n, true)

end

end

{{out}}

<pre> 2 3 5 7 11 13 17 19 23 29

=={{header|CMake}}==

# Check for integer overflow. With CMake using 32-bit signed integer,

# this check fails when limit > 46340.

endforeach(i)

set(${var} ${list} PARENT_SCOPE)

# Print all prime numbers through 100.

eratosthenes(primes 100)

=={{header|COBOL}}==

*> which are kludges to get syntax highlighting to work.

IDENTIFICATION DIVISION.

GOBACK

=={{header|Comal}}==

{{trans|BASIC}}

input "Limit? ": limit

dim sieve(1:limit)

endif

endfor i

{{Out}}

=={{header|Common Lisp}}==

(loop

with sieve = (make-array (1+ maximum)

and do (loop for composite from (expt candidate 2)

to maximum by candidate

Working with odds only (above twice speedup), and marking composites only for primes up to the square root of the maximum:

"Prime numbers sieve for odd numbers.

Returns a list with all the primes that are less than or equal to maximum."

:do (loop :for j :from (* i (1+ i) 2) :to maxi :by odd-number

:do (setf (sbit sieve j) 1))

The indexation scheme used here interprets each index <code>i</code> as standing for the value <code>2i+1</code>. Bit <code>0</code> is unused, a small price to pay for the simpler index calculations compared with the <code>2i+3</code> indexation scheme. The multiples of a given odd prime <code>p</code> are enumerated in increments of <code>2p</code>, which corresponds to the index increment of <code>p</code> on the sieve array. The starting point <code>p*p = (2i+1)(2i+1) = 4i(i+1)+1</code> corresponds to the index <code>2i(i+1)</code>.

=={{header|Cowgol}}==

# To change the maximum prime, change the size of this array

end if;

cand := cand + 1;

{{out}}

This implementation uses a `BitArray` so it is automatically bit-packed to use just one bit per number representation:

require "bit_array"

print SoE.new(1_000_000_000).each.size

elpsd = (Time.monotonic - start_time).total_milliseconds

{{out}}

the non-odds-only version as per the above should never be used because in not using odds-only, it uses twice the memory and over two and a half times the CPU operations as the following odds-only code, which is very little more complex:

require "bit_array"

print SoE_Odds.new(1_000_000_000).each.size

elpsd = (Time.monotonic - start_time).total_milliseconds

{{out}}