Sieve of Eratosthenes: Difference between revisions
Content deleted Content added
Improve C++ style for stdlib version. |
No edit summary |
||
(58 intermediate revisions by 25 users not shown) | |||
Line 761:
ENDLOOP.
</syntaxhighlight>
=={{header|ABC}}==
<syntaxhighlight lang="ABC">HOW TO SIEVE UP TO n:
SHARE sieve
PUT {} IN sieve
FOR cand IN {2..n}: PUT 1 IN sieve[cand]
FOR cand IN {2..floor root n}:
IF sieve[cand] = 1:
PUT cand*cand IN comp
WHILE comp <= n:
PUT 0 IN sieve[comp]
PUT comp+cand IN comp
HOW TO REPORT prime n:
SHARE sieve
IF n<2: FAIL
REPORT sieve[n] = 1
SIEVE UP TO 100
FOR n IN {1..100}:
IF prime n: WRITE n</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|ACL2}}==
Line 1,287 ⟶ 1,310:
The required list of prime divisors obtains by recursion (<tt>{⍵/⍳⍴⍵}∇⌈⍵*0.5</tt>).
=== Optimized Version ===|
<syntaxhighlight lang="apl">sieve←{
Line 2,248 ⟶ 2,271:
@ ^ p3\" ":<
2 234567890123456789012345678901234567890123456789012345678901234567890123456789
=={{header|Binary Lambda Calculus}}==
The BLC sieve of Eratosthenes as documented at https://github.com/tromp/AIT/blob/master/characteristic_sequences/primes.lam is the 167 bit program
<pre>00010001100110010100011010000000010110000010010001010111110111101001000110100001110011010000000000101101110011100111111101111000000001111100110111000000101100000110110</pre>
The infinitely long output is
<pre>001101010001010001010001000001010000010001010001000001000001010000010001010000010001000001000000010001010001010001000000000000010001000001010000000001010000010000010001000001000001010000000001010001010000000000010000000000010001010001000001010000000001000001000001000001010000010001010000000001000000000000010001010001000000000000010000010000000001010001000001000...</pre>
=={{header|BQN}}==
Line 2,794 ⟶ 2,827:
{
auto prime = std::distance(start, it) + 1;
// mark all multiples of this prime number as composite.
while (std::distance(multiple_it, end) > prime)
{
std::advance(
*
}
}
Line 4,666 ⟶ 4,699:
4987
4999</pre>
=={{header|Craft Basic}}==
<syntaxhighlight lang="basic">define limit = 120
dim flags[limit]
for n = 2 to limit
let flags[n] = 1
next n
print "prime numbers less than or equal to ", limit ," are:"
for n = 2 to sqrt(limit)
if flags[n] = 1 then
for i = n * n to limit step n
let flags[i] = 0
next i
endif
next n
for n = 1 to limit
if flags[n] then
print n
endif
next n</syntaxhighlight>
{{out| Output}}<pre>
prime numbers less than or equal to 120 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 </pre>
=={{header|Crystal}}==
Line 5,484 ⟶ 5,557:
The algorithm can be sped up by a factor of four by extreme wheel factorization and (likely) about a factor of the effective number of CPU cores by using multi-processing isolates, but there isn't much point if one is to use the prime generator for output. For most purposes, it is better to use custom functions that directly manipulate the culled bit-packed page segments as `countPrimesTo` does here.
=={{header|dc}}==
<syntaxhighlight lang="dc">[dn[,]n dsx [d 1 r :a lx + d ln!<.] ds.x lx] ds@
[sn 2 [d;a 0=@ 1 + d ln!<#] ds#x] se
100 lex</syntaxhighlight>
{{out}}
<pre>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,</pre>
=={{header|Delphi}}==
Line 5,714 ⟶ 5,797:
=={{header|EasyLang}}==
<syntaxhighlight
len
max = sqrt len
for
if
sieve[j] = 1
j += i
.
.
.
for i = 2 to len
if
print i
.
.
</syntaxhighlight>
Line 6,369 ⟶ 6,454:
=={{header|Elm}}==
===Elm with immutable arrays===
<syntaxhighlight lang="elm">
module PrimeArray exposing (main)
import Array exposing (Array, foldr, map, set)
import Html exposing (div, h1, p, text)
import Html.Attributes exposing (style)
{-
The Eratosthenes sieve task in Rosetta Code does not accept the use of modulo function (allthough Elm functions modBy and remainderBy work always correctly as they require type Int excluding type Float). Thus the solution needs an indexed work array as Elm has no indexes for lists.
In this method we need no division remainder calculations, as we just set the markings of non-primes into the array. We need the indexes that we know, where the marking of the non-primes shall be set.
Because everything is immutable in Elm, every change of array values will create a new array save the original array unchanged. That makes the program running slower or consuming more space of memory than with non-functional imperative languages. All conventional loops (for, while, until) are excluded in Elm because immutability requirement.
Live: https://ellie-app.com/pTHJyqXcHtpa1
-}
alist =
List.range 2 150
-- Work array contains integers 2 ... 149
workArray =
Array.fromList alist
n : Int
n =
List.length alist
-- The max index of integers used in search for primes
-- limit * limit < n
-- equal: limit <= √n
limit : Int
limit =
round (0.5 + sqrt (toFloat n))
-- Remove zero cells of the array
findZero : Int -> Bool
findZero =
\el -> el > 0
zeroFree : Array Int
zeroFree =
Array.filter findZero workResult
nrFoundPrimes =
Array.length zeroFree
workResult : Array Int
workResult =
loopI 2 limit workArray
{- As Elm has no loops (for, while, until)
we must use recursion instead!
The search of prime starts allways saving the
first found value (not setting zero) and continues setting the multiples of prime to zero.
Zero is no integer and may thus be used as marking of non-prime numbers. At the end, only the primes remain in the array and the zeroes are removed from the resulted array to be shown in Html.
-}
-- The recursion increasing variable i follows:
loopI : Int -> Int -> Array Int -> Array Int
loopI i imax arr =
if i > imax then
arr
else
let
arr2 =
phase i arr
in
loopI (i + 1) imax arr2
-- The helper function
phase : Int -> Array Int -> Array Int
phase i =
arrayMarker i (2 * i - 2) n
lastPrime =
Maybe.withDefault 0 <| Array.get (nrFoundPrimes - 1) zeroFree
outputArrayInt : Array Int -> String
outputArrayInt arr =
decorateString <|
foldr (++) "" <|
Array.map (\x -> String.fromInt x ++ " ") arr
decorateString : String -> String
decorateString str =
"[ " ++ str ++ "]"
-- Recursively marking the multiples of p with zero
-- This loop operates with constant p
arrayMarker : Int -> Int -> Int -> Array Int -> Array Int
arrayMarker p min max arr =
let
arr2 =
set min 0 arr
min2 =
min + p
in
if min < max then
arrayMarker p min2 max arr2
else
arr
main =
div [ style "margin" "2%" ]
[ h1 [] [ text "Sieve of Eratosthenes" ]
, text ("List of integers [2, ... ," ++ String.fromInt n ++ "]")
, p [] [ text ("Total integers " ++ String.fromInt n) ]
, p [] [ text ("Max prime of search " ++ String.fromInt limit) ]
, p [] [ text ("The largest found prime " ++ String.fromInt lastPrime) ]
, p [ style "color" "blue", style "font-size" "1.5em" ]
[ text (outputArrayInt zeroFree) ]
, p [] [ text ("Found " ++ String.fromInt nrFoundPrimes ++ " primes") ]
]
</syntaxhighlight>
{{out}}
<pre>
List of integers [2, ... ,149]
Total integers 149
Max prime of search 13
The largest found prime 149
[ 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 ]
Found 35 primes </pre>
===Concise Elm Immutable Array Version===
Although functional, the above code is written in quite an imperative style, so the following code is written in a more concise functional style and includes timing information for counting the number of primes to a million:
<syntaxhighlight lang="elm">module Main exposing (main)
import Browser exposing (element)
import Task exposing (Task, succeed, perform, andThen)
import Html exposing (Html, div, text)
import Time exposing (now, posixToMillis)
import Array exposing (repeat, get, set)
cLIMIT : Int
cLIMIT = 1000000
primesArray : Int -> List Int
primesArray n =
if n < 2 then [] else
let
sz = n + 1
loopbp bp arr =
let s = bp * bp in
if s >= sz then arr else
let tst = get bp arr |> Maybe.withDefault True in
if tst then loopbp (bp + 1) arr else
let
cullc c iarr =
if c >= sz then iarr else
cullc (c + bp) (set c True iarr)
in loopbp (bp + 1) (cullc s arr)
cmpsts = loopbp 2 (repeat sz False)
cnvt (i, t) = if t then Nothing else Just i
in cmpsts |> Array.toIndexedList
|> List.drop 2 -- skip the values for zero and one
|> List.filterMap cnvt -- primes are indexes of not composites
type alias Model = List String
type alias Msg = Model
test : (Int -> List Int) -> Int -> Cmd Msg
test primesf lmt =
let
to100 = primesf 100 |> List.map String.fromInt |> String.join ", "
to100str = "The primes to 100 are: " ++ to100
timemillis() = now |> andThen (succeed << posixToMillis)
in timemillis() |> andThen (\ strt ->
let cnt = primesf lmt |> List.length
in timemillis() |> andThen (\ stop ->
let answrstr = "Found " ++ (String.fromInt cnt) ++ " primes to "
++ (String.fromInt cLIMIT) ++ " in "
++ (String.fromInt (stop - strt)) ++ " milliseconds."
in succeed [to100str, answrstr] ) ) |> perform identity
main : Program () Model Msg
main =
element { init = \ _ -> ( [], test primesArray cLIMIT )
, update = \ msg _ -> (msg, Cmd.none)
, subscriptions = \ _ -> Sub.none
, view = div [] << List.map (div [] << List.singleton << text) }</syntaxhighlight>
{{out}}
<pre>The primes up to 100 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.
Found 78498 primes to 1000000 in 958 milliseconds.</pre>
The above output is the contents of the HTML web page as shown with Google Chrome version 1.23 running on an AMD 7840HS CPU at 5.1 GHz (single thread boosted).
===Concise Elm Immutable Array Odds-Only Version===
The following code can replace the `primesArray` function in the above program and called from the testing and display code (two places):
<syntaxhighlight lang="elm">primesArrayOdds : Int -> List Int
primesArrayOdds n =
if n < 2 then [] else
let
sz = (n - 1) // 2
loopi i arr =
let s = (i + i) * (i + 3) + 3 in
if s >= sz then arr else
let tst = get i arr |> Maybe.withDefault True in
if tst then loopi (i + 1) arr else
let
bp = i + i + 3
cullc c iarr =
if c >= sz then iarr else
cullc (c + bp) (set c True iarr)
in loopi (i + 1) (cullc s arr)
cmpsts = loopi 0 (repeat sz False)
cnvt (i, t) = if t then Nothing else Just <| i + i + 3
oddprms = cmpsts |> Array.toIndexedList |> List.filterMap cnvt
in 2 :: oddprms</syntaxhighlight>
{{out}}
<pre>The primes up to 100 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.
Found 78498 primes to 1000000 in 371 milliseconds.</pre>
The above output is the contents of the HTML web page as shown with Google Chrome version 1.23 running on an AMD 7840HS CPU at 5.1 GHz (single thread boosted).
===Richard Bird Tree Folding Elm Version===
The Elm language doesn't efficiently handle the Sieve of Eratosthenes (SoE) algorithm because it doesn't have directly accessible linear arrays (the Array module used above is based on a persistent tree of sub arrays) and also does Copy On Write (COW) for every write to every location as well as a logarithmic process of updating as a "tree" to minimize the COW operations. Thus, there is better performance implementing the Richard Bird Tree Folding functional algorithm, as follows:
{{trans|Haskell}}
<syntaxhighlight lang="elm">module Main exposing (main)
Line 6,377 ⟶ 6,728:
import Html exposing (Html, div, text)
import Time exposing (now, posixToMillis)
cLIMIT : Int
cLIMIT = 1000000
type CIS a = CIS a (() -> CIS a)
Line 6,419 ⟶ 6,773:
in CIS 2 <| \ () -> oddprms()
type alias Model =
type alias Msg =
test : (() -> CIS Int) -> Int -> Cmd Msg
test primesf lmt =
let
to100 = primesf() |> uptoCIS2List 100
|> List.map String.fromInt |> String.join ", "
to100str = "The primes to 100 are: " ++ to100
timemillis() = now |> andThen (succeed << posixToMillis)
in timemillis() |> andThen (\ strt ->
let cnt =
in timemillis() |> andThen (\
let answrstr = "Found " ++ (String.fromInt cnt) ++ " primes to "
++ (String.fromInt cLIMIT) ++ " in "
++ (String.fromInt (stop - strt)) ++ " milliseconds."
in succeed [to100str, answrstr] ) ) |> perform identity
main : Program () Model Msg
main =
element { init = \ _ -> (
, update = \ msg _ ->
, subscriptions = \ _ -> Sub.none
, view =
{{out}}
<pre>The primes up to 100 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.
Found 78498 primes to 1000000 in
The above output is the contents of the HTML web page as shown with Google Chrome version 1.23 running on an AMD 7840HS CPU at 5.1 GHz (single thread boosted).
===Elm Priority Queue Version===
Using a Binary Minimum Heap Priority Queue is a constant factor faster than the above code as the data structure is balanced rather than "heavy to the right" and requires less memory allocations/deallocation in the following code, which implements enough of the Priority Queue for the purpose. Just substitute the following code for `primesTreeFolding` and pass `primesPQ` as an argument to `test` rather than `primesTreeFolding`:
<syntaxhighlight lang="elm">type PriorityQ comparable v =
Mt
| Br comparable v (PriorityQ comparable v)
Line 6,553 ⟶ 6,866:
else CIS n <| \ () -> sieve (n + 2) pq q bps
oddprms() = CIS 3 <| \ () -> sieve 5 emptyPQ 9 <| oddprms()
in CIS 2 <| \ () -> oddprms()</syntaxhighlight>
{{out}}
<pre>The primes up to 100 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.
Found 78498 primes to 1000000 in
The above output is the contents of the HTML web page as shown with Google Chrome version 1.23 running on an AMD 7840HS CPU at 5.1 GHz (single thread boosted).
=={{header|Emacs Lisp}}==
Line 6,848 ⟶ 7,123:
another several erlang implementation: http://mijkenator.github.io/2015/11/29/project-euler-problem-10/
===Erlang using lists:seq\3 for the initial list and the lists of multiples to be removed ===
<syntaxhighlight lang="erlang">
-module(primesieve).
-export([primes/1]).
mult(N, Limit) ->
case Limit > N * N of
true -> lists:seq(N * N, Limit, N);
false -> []
end.
primes(Limit) ->
case Limit > 1 of
true -> sieve(Limit, 3, [2] ++ lists:seq(3, Limit, 2), mult(3, Limit));
false -> []
end.
sieve(Limit, D, S, M) ->
case Limit < D * D of
true -> S;
false -> sieve(Limit, D + 2, S -- M, mult(D + 2, Limit))
end.
</syntaxhighlight>
{{out}}
<pre>
2> primesieve:primes(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]
3> timer:tc(primesieve, primes, [100]).
{20,
[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]}
</pre>
=={{header|ERRE}}==
Line 6,891 ⟶ 7,201:
16301 16319 16333 16339 16349 16361 16363 16369 16381
1899 PRIMES
</pre>
=={{header|Euler}}==
The original Euler doesn't have loops built-in. Loops can easily be added by defining and calling suitable procedures with literal procedures as parameters. In this sample, a C-style "for" loop procedure is defined and used to sieve and print the primes.<br>
'''begin'''
'''new''' sieve; '''new''' for; '''new''' prime; '''new''' i;
for <- ` '''formal''' init; '''formal''' test; '''formal''' incr; '''formal''' body;
'''begin'''
'''label''' again;
init;
again: '''if''' test '''then''' '''begin''' body; incr; '''goto''' again '''end''' '''else''' 0
'''end'''
'
;
sieve <- ` '''formal''' n;
'''begin'''
'''new''' primes; '''new''' i; '''new''' i2; '''new''' j;
primes <- '''list''' n;
for( ` i <- 1 ', ` i <= n ', ` i <- i + 1 '
, ` primes[ i ] <- '''true''' '
);
primes[ 1 ] <- '''false''';
for( ` i <- 2 '
, ` [ i2 <- i * i ] <= n '
, ` i <- i + 1 '
, ` '''if''' primes[ i ] '''then'''
for( ` j <- i2 ', ` j <= n ', ` j <- j + i '
, ` primes[ j ] <- '''false''' '
)
'''else''' 0
'
);
primes
'''end'''
'
;
prime <- sieve( 30 );
for( ` i <- 1 ', ` i <= '''length''' prime ', ` i <- i + 1 '
, ` '''if''' prime[ i ] '''then''' '''out''' i '''else''' 0 '
)
'''end''' $
{{out}}
<pre>
NUMBER 2
NUMBER 3
NUMBER 5
NUMBER 7
NUMBER 11
NUMBER 13
NUMBER 17
NUMBER 19
NUMBER 23
NUMBER 29
</pre>
Line 8,160 ⟶ 8,528:
I: DWord;
begin
Write(aList[I], ', ');
WriteLn(aList[aList.Size - 1]);
end;
aList.Free;
end;
Line 8,301 ⟶ 8,672:
for I := 0 to Length(aList) - 2 do
Write(aList[I], ', ');
WriteLn(aList[High(aList)]);
end;
begin
Line 8,410 ⟶ 8,782:
</syntaxhighlight>
=={{header|Peri}}==
''Note: With benchmark function''
<syntaxhighlight lang="peri">
###sysinclude standard.uh
tick sto startingtick
#g 100000 sto MAX
@MAX mem !maximize sto primeNumbers
one count
2 0 sto#s primeNumbers
2 @MAX külső: {{ ,
@count {{
{{}}§külső primeNumbers[{{}}] !/ else {{<}}§külső
}} // @count vége
//{{}} gprintnl // A talált prímszám kiiratásához kommentezzük ki e sort
{{}} @count++ sto#s primeNumbers
}} // @MAX vége
@primeNumbers inv mem
//."Time : " tick @startingtick - print ." tick\n"
."Prímek száma = " @count printnl
end
{ „MAX” } { „startingtick” } { „primeNumbers” } { „count” }
</syntaxhighlight>
=={{header|FutureBasic}}==
Line 8,442 ⟶ 8,837:
=={{header|Fōrmulæ}}==
{{FormulaeEntry|page=https://formulae.org/?script=examples/Sieve_of_Eratosthenes}}
'''Solution'''
[[File:Fōrmulæ - Sieve of Eratosthenes 01.png]]
'''Test case'''
[[File:Fōrmulæ - Sieve of Eratosthenes 02.png]]
[[File:Fōrmulæ - Sieve of Eratosthenes 03.png]]
=={{header|GAP}}==
Line 8,845 ⟶ 9,246:
var p int
p = <-primes
Line 9,066 ⟶ 9,467:
a = accumArray (\ b c -> False) True (1,q-1)
[(i,()) | (s,y) <- fs, i <- [y+s, y+s+s..q]]</syntaxhighlight>
====As list comprehension====
<syntaxhighlight lang="haskell">import Data.Array.Unboxed
import Data.List (tails, inits)
primes = 2 : [ n |
(r:q:_, px) <- zip (tails (2 : [p*p | p <- primes]))
(inits primes),
(n, True) <- assocs ( accumArray (\_ _ -> False) True
(r+1,q-1)
[ (m,()) | p <- px
, s <- [ div (r+p) p * p]
, m <- [s,s+p..q-1] ] :: UArray Int Bool
) ]</syntaxhighlight>
===Basic list-based sieve===
Line 10,312 ⟶ 10,728:
def eratosthenes:
def erase(i):
if .[i] then
reduce (range(2*i; length; i)) as $j (.; .[$j] = false)
else .
end;
Line 10,405 ⟶ 10,821:
Alternate version using <code>findall</code> to get all primes at once in the end
<syntaxhighlight lang="julia">function sieve(n
for p in 2:n
primes[p .* (2:n÷p)]
end
findall(primes)
end</syntaxhighlight>
Line 11,279 ⟶ 11,688:
=={{header|Lambdatalk}}==
<syntaxhighlight lang="scheme">
• 1) create an array of natural numbers, [0,1,2,3, ... ,n-1]
• 2) the 3rd number is 2, we set to dots all its composites by steps of 2,
• 3) the 4th number is 3, we set to dots all its composites by steps of 3,
• 4) the 6th number is 5, we set to dots all its composites by steps of 5,
• 5) the remaining numbers are primes and we clean all dots.
For instance:
1: 0 0 0 0 0 0 0 0 0 9 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3
0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0
2: 0 1 2 3 . 5 . 7 . 9 . 1 . 3 . 5 . 7 . 9 . 1 . 3 . 5 . 7 . 9 .
3: 0 1 2 3 . 5 . 7 . . . 1 . 3 . . . 7 . 9 . . . 3 . 5 . . . 9 .
4: 0 1 2 3 . 5 . 7 . .
| | | | | | | | | |
5: 0 0 0 0 1 1 1 1 2 2
2 3 5 7 1 3 7 9 3 9
1) recursive version {rsieve n}
{def rsieve
{def rsieve.mark
{lambda {:n :a :i :j}
{if {< :j :n}
then {rsieve.mark :n
{A.set! :j . :a}
:i
{+ :i :j}}
else :a}}}
{def rsieve.loop
{lambda {:n :a :i}
{if {< {* :i :i} :n}
then {rsieve.loop :n
{if {W.equal? {A.get :i :a} .}
then :a
else {rsieve.mark :n :a :i {* :i :i}}}
{+ :i 1}}
else :a}}}
{lambda {:n}
{S.replace \s by space in
{S.replace (\[|\]|\.|,) by space in
{A.disp
{A.slice 2 -1
{rsieve.loop :n
{A.new {S.serie 0 :n}} 2}}}}}}}
-> rsieve
{rsieve 1000}
-> 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
Note: this version doesn't avoid stackoverflow.
2) iterative version {isieve n}
{def isieve
{def isieve.mark
{lambda {:n :a :i}
{S.map {{lambda {:a :j} {A.set! :j . :a}
} :a}
{S.serie {* :i :i} :n :i} }}}
{lambda {:n}
{S.replace \s by space in
{S.replace (\[|\]|\.|,) by space in
{A.disp
{A.slice 2 -1
{S.last
{S.map {{lambda {:n :a :i} {if {W.equal? {A.get :i :a} .}
then
else {isieve.mark :n :a :i}}
} :n {A.new {S.serie 0 :n}}}
{S.serie 2 {sqrt :n} 1}}}}}}}}}
-> isieve
{isieve 1000}
-> 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
Notes:
- this version avoids stackoverflow.
- From 1 to 1000000 there are 78500 primes (computed in ~15000ms) and the last is 999983.
</syntaxhighlight>
=={{header|langur}}==
{{trans|D}}
<syntaxhighlight lang="langur">
val sieve =
if limit < 2: return []
var
for
if not
for
}
}
}
filter
}
writeln
</syntaxhighlight>
{{out}}
<pre>[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]</pre>
=={{header|LFE}}==
<syntaxhighlight lang="lisp">
(defmodule eratosthenes
(export (sieve 1)))
(defun sieve (limit)
(sieve limit (lists:seq 2 limit)))
(defun sieve
((limit (= l (cons p _))) (when (> (* p p) limit))
l)
((limit (cons p ns))
(cons p (sieve limit (remove-multiples p (* p p) ns)))))
(defun remove-multiples (p q l)
(lists:reverse (remove-multiples p q l '())))
(defun remove-multiples
((_ _ '() s) s)
((p q (cons q ns) s)
(remove-multiples p q ns s))
((p q (= r (cons a _)) s) (when (> a q))
(remove-multiples p (+ q p) r s))
((p q (cons n ns) s)
(remove-multiples p q ns (cons n s))))
</syntaxhighlight>
{{Out}}
<pre>
lfe> (slurp "sieve.lfe")
#(ok eratosthenes)
lfe> (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)
</pre>
=={{header|Liberty BASIC}}==
Line 12,175 ⟶ 12,673:
IO.Put("\n");
END Sieve.</syntaxhighlight>
=={{header|Mojo}}==
Tested with Mojo version 0.7:
<syntaxhighlight lang="mojo">from memory import memset_zero
from memory.unsafe import (DTypePointer)
from time import (now)
alias cLIMIT: Int = 1_000_000_000
struct SoEBasic(Sized):
var len: Int
var cmpsts: DTypePointer[DType.bool] # because DynamicVector has deep copy bug in mojo version 0.7
var sz: Int
var ndx: Int
fn __init__(inout self, limit: Int):
self.len = limit - 1
self.sz = limit - 1
self.ndx = 0
self.cmpsts = DTypePointer[DType.bool].alloc(limit - 1)
memset_zero(self.cmpsts, limit - 1)
for i in range(limit - 1):
let s = i * (i + 4) + 2
if s >= limit - 1: break
if self.cmpsts[i]: continue
let bp = i + 2
for c in range(s, limit - 1, bp):
self.cmpsts[c] = True
for i in range(limit - 1):
if self.cmpsts[i]: self.sz -= 1
fn __del__(owned self):
self.cmpsts.free()
fn __copyinit__(inout self, existing: Self):
self.len = existing.len
self.cmpsts = DTypePointer[DType.bool].alloc(self.len)
for i in range(self.len):
self.cmpsts[i] = existing.cmpsts[i]
self.sz = existing.sz
self.ndx = existing.ndx
fn __moveinit__(inout self, owned existing: Self):
self.len = existing.len
self.cmpsts = existing.cmpsts
self.sz = existing.sz
self.ndx = existing.ndx
fn __len__(self: Self) -> Int: return self.sz
fn __iter__(self: Self) -> Self: return self
fn __next__(inout self: Self) -> Int:
if self.ndx >= self.len: return 0
while (self.ndx < self.len) and (self.cmpsts[self.ndx]):
self.ndx += 1
let rslt = self.ndx + 2; self.sz -= 1; self.ndx += 1
return rslt
fn main():
print("The primes to 100 are:")
for prm in SoEBasic(100): print_no_newline(prm, " ")
print()
let strt0 = now()
let answr0 = len(SoEBasic(1_000_000))
let elpsd0 = (now() - strt0) / 1000000
print("Found", answr0, "primes up to 1,000,000 in", elpsd0, "milliseconds.")
let strt1 = now()
let answr1 = len(SoEBasic(cLIMIT))
let elpsd1 = (now() - strt1) / 1000000
print("Found", answr1, "primes up to", cLIMIT, "in", elpsd1, "milliseconds.")</syntaxhighlight>
{{out}}
<pre>The primes to 100 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
Found 78498 primes up to 1,000,000 in 1.2642770000000001 milliseconds.
Found 50847534 primes up to 1000000000 in 6034.328751 milliseconds.</pre>
as run on an AMD 7840HS CPU at 5.1 GHz.
Note that due to the huge memory array used, when large ranges are selected, the speed is disproportional in speed slow down by about four times.
This solution uses an interator struct which seems to be the Mojo-preferred way to do this, and normally a DynamicVector would have been used the the culling array except that there is a bug in this version of DynamicVector where the array is not properly deep copied when copied to a new location, so the raw pointer type is used.
===Odds-Only with Optimizations===
This version does three significant improvements to the above code as follows:
1) It is trivial to skip the processing to store representations for and cull the even comosite numbers other than the prime number two, saving half the storage space and reducing the culling time to about 40 percent.
2) There is a repeating pattern of culling composite representations over a bit-packed byte array (which reduces the storage requirement by another eight times) that repeats every eight culling operations, which can be encapsulated by a extreme loop unrolling technique with compiler generated constants as done here.
3) Further, there is a further extreme optimization technique of dense culling for small base prime values whose culling span is less than one register in size where the loaded register is repeatedly culled for different base prime strides before being written out (with such optimization done by the compiler), again using compiler generated modification constants. This technique is usually further optimizated by modern compilers to use efficient autovectorization and the use of SIMD registers available to the architecture to reduce these culling operations to an avererage of a tiny fraction of a CPU clock cycle per cull.
Mojo version 0.7 was tested:
<syntaxhighlight lang="mojo">from memory import (memset_zero, memcpy)
from memory.unsafe import (DTypePointer)
from math.bit import ctpop
from time import (now)
alias cLIMIT: Int = 1_000_000_000
alias cBufferSize: Int = 262144 # bytes
alias cBufferBits: Int = cBufferSize * 8
alias UnrollFunc = fn(DTypePointer[DType.uint8], Int, Int, Int) -> None
@always_inline
fn extreme[OFST: Int, BP: Int](pcmps: DTypePointer[DType.uint8], bufsz: Int, s: Int, bp: Int):
var cp = pcmps + (s >> 3)
let r1: Int = ((s + bp) >> 3) - (s >> 3)
let r2: Int = ((s + 2 * bp) >> 3) - (s >> 3)
let r3: Int = ((s + 3 * bp) >> 3) - (s >> 3)
let r4: Int = ((s + 4 * bp) >> 3) - (s >> 3)
let r5: Int = ((s + 5 * bp) >> 3) - (s >> 3)
let r6: Int = ((s + 6 * bp) >> 3) - (s >> 3)
let r7: Int = ((s + 7 * bp) >> 3) - (s >> 3)
let plmt: DTypePointer[DType.uint8] = pcmps + bufsz - r7
while cp < plmt:
cp.store(cp.load() | (1 << OFST))
(cp + r1).store((cp + r1).load() | (1 << ((OFST + BP) & 7)))
(cp + r2).store((cp + r2).load() | (1 << ((OFST + 2 * BP) & 7)))
(cp + r3).store((cp + r3).load() | (1 << ((OFST + 3 * BP) & 7)))
(cp + r4).store((cp + r4).load() | (1 << ((OFST + 4 * BP) & 7)))
(cp + r5).store((cp + r5).load() | (1 << ((OFST + 5 * BP) & 7)))
(cp + r6).store((cp + r6).load() | (1 << ((OFST + 6 * BP) & 7)))
(cp + r7).store((cp + r7).load() | (1 << ((OFST + 7 * BP) & 7)))
cp += bp
let eplmt: DTypePointer[DType.uint8] = plmt + r7
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << OFST))
cp += r1
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + BP) & 7)))
cp += r2 - r1
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 2 * BP) & 7)))
cp += r3 - r2
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 3 * BP) & 7)))
cp += r4 - r3
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 4 * BP) & 7)))
cp += r5 - r4
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 5 * BP) & 7)))
cp += r6 - r5
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 6 * BP) & 7)))
cp += r7 - r6
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 7 * BP) & 7)))
fn mkExtrm[CNT: Int](pntr: Pointer[UnrollFunc]):
@parameter
if CNT >= 32:
return
alias OFST = CNT >> 2
alias BP = ((CNT & 3) << 1) + 1
pntr.offset(CNT).store(extreme[OFST, BP])
mkExtrm[CNT + 1](pntr)
@always_inline
fn mkExtremeFuncs() -> Pointer[UnrollFunc]:
let jmptbl: Pointer[UnrollFunc] = Pointer[UnrollFunc].alloc(32)
mkExtrm[0](jmptbl)
return jmptbl
let extremeFuncs = mkExtremeFuncs()
alias DenseFunc = fn(DTypePointer[DType.uint64], Int, Int) -> DTypePointer[DType.uint64]
fn mkDenseCull[N: Int, BP: Int](cp: DTypePointer[DType.uint64]):
@parameter
if N >= 64:
return
alias MUL = N * BP
var cop = cp.offset(MUL >> 6)
cop.store(cop.load() | (1 << (MUL & 63)))
mkDenseCull[N + 1, BP](cp)
@always_inline
fn denseCullFunc[BP: Int](pcmps: DTypePointer[DType.uint64], bufsz: Int, s: Int) -> DTypePointer[DType.uint64]:
var cp: DTypePointer[DType.uint64] = pcmps + (s >> 6)
let plmt = pcmps + (bufsz >> 3) - BP
while cp < plmt:
mkDenseCull[0, BP](cp)
cp += BP
return cp
fn mkDenseFunc[CNT: Int](pntr: Pointer[DenseFunc]):
@parameter
if CNT >= 64:
return
alias BP = (CNT << 1) + 3
pntr.offset(CNT).store(denseCullFunc[BP])
mkDenseFunc[CNT + 1](pntr)
@always_inline
fn mkDenseFuncs() -> Pointer[DenseFunc]:
let jmptbl : Pointer[DenseFunc] = Pointer[DenseFunc].alloc(64)
mkDenseFunc[0](jmptbl)
return jmptbl
let denseFuncs : Pointer[DenseFunc] = mkDenseFuncs()
@always_inline
fn cullPass(cmpsts: DTypePointer[DType.uint8], bytesz: Int, s: Int, bp: Int):
if bp <= 129: # dense culling
var sm = s
while (sm >> 3) < bytesz and (sm & 63) != 0:
cmpsts[sm >> 3] |= (1 << (sm & 7))
sm += bp
let bcp = denseFuncs[(bp - 3) >> 1](cmpsts.bitcast[DType.uint64](), bytesz, sm)
var ns = 0
var ncp = bcp
let cmpstslmtp = (cmpsts + bytesz).bitcast[DType.uint64]()
while ncp < cmpstslmtp:
ncp[0] |= (1 << (ns & 63))
ns += bp
ncp = bcp + (ns >> 6)
else: # extreme loop unrolling culling
extremeFuncs[((s & 7) << 2) + ((bp & 7) >> 1)](cmpsts, bytesz, s, bp)
# for c in range(s, self.len, bp): # slow bit twiddling way
# self.cmpsts[c >> 3] |= (1 << (c & 7))
fn countPagePrimes(ptr: DTypePointer[DType.uint8], bitsz: Int) -> Int:
let wordsz: Int = (bitsz + 63) // 64 # round up to nearest 64 bit boundary
var rslt: Int = wordsz * 64
let bigcmps = ptr.bitcast[DType.uint64]()
for i in range(wordsz - 1):
rslt -= ctpop(bigcmps[i]).to_int()
rslt -= ctpop(bigcmps[wordsz - 1] | (-2 << ((bitsz - 1) & 63))).to_int()
return rslt
struct SoEOdds(Sized):
var len: Int
var cmpsts: DTypePointer[DType.uint8] # because DynamicVector has deep copy bug in Mojo version 0.7
var sz: Int
var ndx: Int
fn __init__(inout self, limit: Int):
self.len = 0 if limit < 2 else (limit - 3) // 2 + 1
self.sz = 0 if limit < 2 else self.len + 1 # for the unprocessed only even prime of two
self.ndx = -1
let bytesz = 0 if limit < 2 else ((self.len + 63) & -64) >> 3 # round up to nearest 64 bit boundary
self.cmpsts = DTypePointer[DType.uint8].alloc(bytesz)
memset_zero(self.cmpsts, bytesz)
for i in range(self.len):
let s = (i + i) * (i + 3) + 3
if s >= self.len: break
if (self.cmpsts[i >> 3] >> (i & 7)) & 1 != 0: continue
let bp = i + i + 3
cullPass(self.cmpsts, bytesz, s, bp)
self.sz = countPagePrimes(self.cmpsts, self.len) + 1 # add one for only even prime of two
fn __del__(owned self):
self.cmpsts.free()
fn __copyinit__(inout self, existing: Self):
self.len = existing.len
let bytesz = (self.len + 7) // 8
self.cmpsts = DTypePointer[DType.uint8].alloc(bytesz)
memcpy(self.cmpsts, existing.cmpsts, bytesz)
self.sz = existing.sz
self.ndx = existing.ndx
fn __moveinit__(inout self, owned existing: Self):
self.len = existing.len
self.cmpsts = existing.cmpsts
self.sz = existing.sz
self.ndx = existing.ndx
fn __len__(self: Self) -> Int: return self.sz
fn __iter__(self: Self) -> Self: return self
@always_inline
fn __next__(inout self: Self) -> Int:
if self.ndx < 0:
self.ndx = 0; self.sz -= 1; return 2
while (self.ndx < self.len) and ((self.cmpsts[self.ndx >> 3] >> (self.ndx & 7)) & 1 != 0):
self.ndx += 1
let rslt = (self.ndx << 1) + 3; self.sz -= 1; self.ndx += 1; return rslt
fn main():
print("The primes to 100 are:")
for prm in SoEOdds(100): print_no_newline(prm, " ")
print()
let strt0 = now()
let answr0 = len(SoEOdds(1_000_000))
let elpsd0 = (now() - strt0) / 1000000
print("Found", answr0, "primes up to 1,000,000 in", elpsd0, "milliseconds.")
let strt1 = now()
let answr1 = len(SoEOdds(cLIMIT))
let elpsd1 = (now() - strt1) / 1000000
print("Found", answr1, "primes up to", cLIMIT, "in", elpsd1, "milliseconds.")</syntaxhighlight>
{{out}}
<pre>The primes to 100 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
Found 78498 primes up to 1,000,000 in 0.085067000000000004 milliseconds.
Found 50847534 primes up to 1000000000 in 1204.866606 milliseconds.</pre>
This was run on the same computer as the above example; notice that while this is much faster than that version, it is still very slow as the sieving range gets large such that the relative processing time for a range that is 1000 times as large is about ten times slower than as might be expected by simple scaling. This is due to the "one huge sieving buffer" algorithm that gets very large with increasing range (and in fact will eventually limit the sieving range that can be used) to exceed the size of CPU cache buffers and thus greatly slow average memory access times.
===Page-Segmented Odds-Only with Optimizations===
While the above version performs reasonably well for small sieving ranges that fit within the CPU caches of a few tens of millions, as one can see it gets much slower with larger ranges and as well its huge RAM memory consumption limits the maximum range over which it can be used. This version solves these problems be breaking the huge sieving array into "pages" that each fit within the CPU cache size and processing each "page" sequentially until the target range is reached. This technique also greatly reduces memory requirements to only that required to store the base prime value representations up to the square root of the range limit (about O(n/log n) storage plus a fixed size page buffer. In this case, the storage for the base primes has been reduced by a constant factor by storing them as single byte deltas from the previous value, which works for ranges up to the 64-bit number range where the biggest gap is two times 192 and since we store only for odd base primes, the gap values are all half values to fit in a single byte.
Currently, Mojo has problems with some functions in the standard libraries such as the integer square root function is not accurate nor does it work for the required integer types so a custom integer square root function is supplied. As well, current Mojo does not support recursion for hardly any useful cases (other than compile time global function recursion), so the `SoeOdds` structure from the previous answer had to be kept to generate the base prime representation table (or this would have had to be generated from scratch within the new `SoEOddsPaged` structure). Finally, it didn't seem to be worth using the `Sized` trait for the new structure as this would seem to sometimes require processing the pages twice, one to obtain the size and once if iteration across the prime values is required.
Tested with Mojo version 0.7:
<syntaxhighlight lang="mojo">from memory import (memset_zero, memcpy)
from memory.unsafe import (DTypePointer)
from math.bit import ctpop
from time import (now)
alias cLIMIT: Int = 1_000_000_000
alias cBufferSize: Int = 262144 # bytes
alias cBufferBits: Int = cBufferSize * 8
fn intsqrt(n: UInt64) -> UInt64:
if n < 4:
if n < 1: return 0 else: return 1
var x: UInt64 = n; var qn: UInt64 = 0; var r: UInt64 = 0
while qn < 64 and (1 << qn) <= n:
qn += 2
var q: UInt64 = 1 << qn
while q > 1:
if qn >= 64:
q = 1 << (qn - 2); qn = 0
else:
q >>= 2
let t: UInt64 = r + q
r >>= 1
if x >= t:
x -= t; r += q
return r
alias UnrollFunc = fn(DTypePointer[DType.uint8], Int, Int, Int) -> None
@always_inline
fn extreme[OFST: Int, BP: Int](pcmps: DTypePointer[DType.uint8], bufsz: Int, s: Int, bp: Int):
var cp = pcmps + (s >> 3)
let r1: Int = ((s + bp) >> 3) - (s >> 3)
let r2: Int = ((s + 2 * bp) >> 3) - (s >> 3)
let r3: Int = ((s + 3 * bp) >> 3) - (s >> 3)
let r4: Int = ((s + 4 * bp) >> 3) - (s >> 3)
let r5: Int = ((s + 5 * bp) >> 3) - (s >> 3)
let r6: Int = ((s + 6 * bp) >> 3) - (s >> 3)
let r7: Int = ((s + 7 * bp) >> 3) - (s >> 3)
let plmt: DTypePointer[DType.uint8] = pcmps + bufsz - r7
while cp < plmt:
cp.store(cp.load() | (1 << OFST))
(cp + r1).store((cp + r1).load() | (1 << ((OFST + BP) & 7)))
(cp + r2).store((cp + r2).load() | (1 << ((OFST + 2 * BP) & 7)))
(cp + r3).store((cp + r3).load() | (1 << ((OFST + 3 * BP) & 7)))
(cp + r4).store((cp + r4).load() | (1 << ((OFST + 4 * BP) & 7)))
(cp + r5).store((cp + r5).load() | (1 << ((OFST + 5 * BP) & 7)))
(cp + r6).store((cp + r6).load() | (1 << ((OFST + 6 * BP) & 7)))
(cp + r7).store((cp + r7).load() | (1 << ((OFST + 7 * BP) & 7)))
cp += bp
let eplmt: DTypePointer[DType.uint8] = plmt + r7
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << OFST))
cp += r1
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + BP) & 7)))
cp += r2 - r1
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 2 * BP) & 7)))
cp += r3 - r2
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 3 * BP) & 7)))
cp += r4 - r3
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 4 * BP) & 7)))
cp += r5 - r4
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 5 * BP) & 7)))
cp += r6 - r5
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 6 * BP) & 7)))
cp += r7 - r6
if eplmt == cp or eplmt < cp: return
cp.store(cp.load() | (1 << ((OFST + 7 * BP) & 7)))
fn mkExtrm[CNT: Int](pntr: Pointer[UnrollFunc]):
@parameter
if CNT >= 32:
return
alias OFST = CNT >> 2
alias BP = ((CNT & 3) << 1) + 1
pntr.offset(CNT).store(extreme[OFST, BP])
mkExtrm[CNT + 1](pntr)
@always_inline
fn mkExtremeFuncs() -> Pointer[UnrollFunc]:
let jmptbl: Pointer[UnrollFunc] = Pointer[UnrollFunc].alloc(32)
mkExtrm[0](jmptbl)
return jmptbl
let extremeFuncs = mkExtremeFuncs()
alias DenseFunc = fn(DTypePointer[DType.uint64], Int, Int) -> DTypePointer[DType.uint64]
fn mkDenseCull[N: Int, BP: Int](cp: DTypePointer[DType.uint64]):
@parameter
if N >= 64:
return
alias MUL = N * BP
var cop = cp.offset(MUL >> 6)
cop.store(cop.load() | (1 << (MUL & 63)))
mkDenseCull[N + 1, BP](cp)
@always_inline
fn denseCullFunc[BP: Int](pcmps: DTypePointer[DType.uint64], bufsz: Int, s: Int) -> DTypePointer[DType.uint64]:
var cp: DTypePointer[DType.uint64] = pcmps + (s >> 6)
let plmt = pcmps + (bufsz >> 3) - BP
while cp < plmt:
mkDenseCull[0, BP](cp)
cp += BP
return cp
fn mkDenseFunc[CNT: Int](pntr: Pointer[DenseFunc]):
@parameter
if CNT >= 64:
return
alias BP = (CNT << 1) + 3
pntr.offset(CNT).store(denseCullFunc[BP])
mkDenseFunc[CNT + 1](pntr)
@always_inline
fn mkDenseFuncs() -> Pointer[DenseFunc]:
let jmptbl : Pointer[DenseFunc] = Pointer[DenseFunc].alloc(64)
mkDenseFunc[0](jmptbl)
return jmptbl
let denseFuncs : Pointer[DenseFunc] = mkDenseFuncs()
@always_inline
fn cullPass(cmpsts: DTypePointer[DType.uint8], bytesz: Int, s: Int, bp: Int):
if bp <= 129: # dense culling
var sm = s
while (sm >> 3) < bytesz and (sm & 63) != 0:
cmpsts[sm >> 3] |= (1 << (sm & 7))
sm += bp
let bcp = denseFuncs[(bp - 3) >> 1](cmpsts.bitcast[DType.uint64](), bytesz, sm)
var ns = 0
var ncp = bcp
let cmpstslmtp = (cmpsts + bytesz).bitcast[DType.uint64]()
while ncp < cmpstslmtp:
ncp[0] |= (1 << (ns & 63))
ns += bp
ncp = bcp + (ns >> 6)
else: # extreme loop unrolling culling
extremeFuncs[((s & 7) << 2) + ((bp & 7) >> 1)](cmpsts, bytesz, s, bp)
# for c in range(s, self.len, bp): # slow bit twiddling way
# self.cmpsts[c >> 3] |= (1 << (c & 7))
fn cullPage(lwi: Int, lmt: Int, cmpsts: DTypePointer[DType.uint8], bsprmrps: DTypePointer[DType.uint8]):
var bp = 1; var ndx = 0
while True:
bp += bsprmrps[ndx].to_int() << 1
let i = (bp - 3) >> 1
var s = (i + i) * (i + 3) + 3
if s >= lmt: break
if s >= lwi: s -= lwi
else:
s = (lwi - s) % bp
if s != 0: s = bp - s
cullPass(cmpsts, cBufferSize, s, bp)
ndx += 1
fn countPagePrimes(ptr: DTypePointer[DType.uint8], bitsz: Int) -> Int:
let wordsz: Int = (bitsz + 63) // 64 # round up to nearest 64 bit boundary
var rslt: Int = wordsz * 64
let bigcmps = ptr.bitcast[DType.uint64]()
for i in range(wordsz - 1):
rslt -= ctpop(bigcmps[i]).to_int()
rslt -= ctpop(bigcmps[wordsz - 1] | (-2 << ((bitsz - 1) & 63))).to_int()
return rslt
struct SoEOdds(Sized):
var len: Int
var cmpsts: DTypePointer[DType.uint8] # because DynamicVector has deep copy bug in Mojo version 0.7
var sz: Int
var ndx: Int
fn __init__(inout self, limit: Int):
self.len = 0 if limit < 2 else (limit - 3) // 2 + 1
self.sz = 0 if limit < 2 else self.len + 1 # for the unprocessed only even prime of two
self.ndx = -1
let bytesz = 0 if limit < 2 else ((self.len + 63) & -64) >> 3 # round up to nearest 64 bit boundary
self.cmpsts = DTypePointer[DType.uint8].alloc(bytesz)
memset_zero(self.cmpsts, bytesz)
for i in range(self.len):
let s = (i + i) * (i + 3) + 3
if s >= self.len: break
if (self.cmpsts[i >> 3] >> (i & 7)) & 1 != 0: continue
let bp = i + i + 3
cullPass(self.cmpsts, bytesz, s, bp)
self.sz = countPagePrimes(self.cmpsts, self.len) + 1 # add one for only even prime of two
fn __del__(owned self):
self.cmpsts.free()
fn __copyinit__(inout self, existing: Self):
self.len = existing.len
let bytesz = (self.len + 7) // 8
self.cmpsts = DTypePointer[DType.uint8].alloc(bytesz)
memcpy(self.cmpsts, existing.cmpsts, bytesz)
self.sz = existing.sz
self.ndx = existing.ndx
fn __moveinit__(inout self, owned existing: Self):
self.len = existing.len
self.cmpsts = existing.cmpsts
self.sz = existing.sz
self.ndx = existing.ndx
fn __len__(self: Self) -> Int: return self.sz
fn __iter__(self: Self) -> Self: return self
@always_inline
fn __next__(inout self: Self) -> Int:
if self.ndx < 0:
self.ndx = 0; self.sz -= 1; return 2
while (self.ndx < self.len) and ((self.cmpsts[self.ndx >> 3] >> (self.ndx & 7)) & 1 != 0):
self.ndx += 1
let rslt = (self.ndx << 1) + 3; self.sz -= 1; self.ndx += 1; return rslt
struct SoEOddsPaged:
var len: Int
var cmpsts: DTypePointer[DType.uint8] # because DynamicVector has deep copy bug in Mojo version 0.7
var sz: Int # 0 means finished; otherwise contains number of odd base primes
var ndx: Int
var lwi: Int
var bsprmrps: DTypePointer[DType.uint8] # contains deltas between odd base primes starting from zero
fn __init__(inout self, limit: UInt64):
self.len = 0 if limit < 2 else ((limit - 3) // 2 + 1).to_int()
self.sz = 0 if limit < 2 else 1 # means iterate until this is set to zero
self.ndx = -1 # for unprocessed only even prime of two
self.lwi = 0
if self.len < cBufferBits:
let bytesz = ((self.len + 63) & -64) >> 3 # round up to nearest 64 bit boundary
self.cmpsts = DTypePointer[DType.uint8].alloc(bytesz)
self.bsprmrps = DTypePointer[DType.uint8].alloc(self.sz)
else:
self.cmpsts = DTypePointer[DType.uint8].alloc(cBufferSize)
let bsprmitr = SoEOdds(intsqrt(limit).to_int())
self.sz = len(bsprmitr)
self.bsprmrps = DTypePointer[DType.uint8].alloc(self.sz)
var ndx = -1; var oldbp = 1
for bsprm in bsprmitr:
if ndx < 0: ndx += 1; continue # skip over the 2 prime
self.bsprmrps[ndx] = (bsprm - oldbp) >> 1
oldbp = bsprm; ndx += 1
self.bsprmrps[ndx] = 255 # one extra value to go beyond the necessary cull space
fn __del__(owned self):
self.cmpsts.free(); self.bsprmrps.free()
fn __copyinit__(inout self, existing: Self):
self.len = existing.len
self.sz = existing.sz
let bytesz = cBufferSize if self.len >= cBufferBits
else ((self.len + 63) & -64) >> 3 # round up to nearest 64 bit boundary
self.cmpsts = DTypePointer[DType.uint8].alloc(bytesz)
memcpy(self.cmpsts, existing.cmpsts, bytesz)
self.ndx = existing.ndx
self.lwi = existing.lwi
self.bsprmrps = DTypePointer[DType.uint8].alloc(self.sz)
memcpy(self.bsprmrps, existing.bsprmrps, self.sz)
fn __moveinit__(inout self, owned existing: Self):
self.len = existing.len
self.cmpsts = existing.cmpsts
self.sz = existing.sz
self.ndx = existing.ndx
self.lwi = existing.lwi
self.bsprmrps = existing.bsprmrps
fn countPrimes(self) -> Int:
if self.len <= cBufferBits: return len(SoEOdds(2 * self.len + 1))
var cnt = 1; var lwi = 0
let cmpsts = DTypePointer[DType.uint8].alloc(cBufferSize)
memset_zero(cmpsts, cBufferSize)
cullPage(0, cBufferBits, cmpsts, self.bsprmrps)
while lwi + cBufferBits <= self.len:
cnt += countPagePrimes(cmpsts, cBufferBits)
lwi += cBufferBits
memset_zero(cmpsts, cBufferSize)
let lmt = lwi + cBufferBits if lwi + cBufferBits <= self.len else self.len
cullPage(lwi, lmt, cmpsts, self.bsprmrps)
cnt += countPagePrimes(cmpsts, self.len - lwi)
return cnt
fn __len__(self: Self) -> Int: return self.sz
fn __iter__(self: Self) -> Self: return self
@always_inline
fn __next__(inout self: Self) -> Int: # don't count number of primes by interating - slooow
if self.ndx < 0:
self.ndx = 0; self.lwi = 0
if self.len < 2: self.sz = 0
elif self.len <= cBufferBits:
let bytesz = ((self.len + 63) & -64) >> 3 # round up to nearest 64 bit boundary
memset_zero(self.cmpsts, bytesz)
for i in range(self.len):
let s = (i + i) * (i + 3) + 3
if s >= self.len: break
if (self.cmpsts[i >> 3] >> (i & 7)) & 1 != 0: continue
let bp = i + i + 3
cullPass(self.cmpsts, bytesz, s, bp)
else:
memset_zero(self.cmpsts, cBufferSize)
cullPage(0, cBufferBits, self.cmpsts, self.bsprmrps)
return 2
let rslt = ((self.lwi + self.ndx) << 1) + 3; self.ndx += 1
if self.lwi + cBufferBits >= self.len:
while (self.lwi + self.ndx < self.len) and ((self.cmpsts[self.ndx >> 3] >> (self.ndx & 7)) & 1 != 0):
self.ndx += 1
else:
while (self.ndx < cBufferBits) and ((self.cmpsts[self.ndx >> 3] >> (self.ndx & 7)) & 1 != 0):
self.ndx += 1
while (self.ndx >= cBufferBits) and (self.lwi + cBufferBits <= self.len):
self.ndx = 0; self.lwi += cBufferBits; memset_zero(self.cmpsts, cBufferSize)
let lmt = self.lwi + cBufferBits if self.lwi + cBufferBits <= self.len else self.len
cullPage(self.lwi, lmt, self.cmpsts, self.bsprmrps)
let buflmt = cBufferBits if self.lwi + cBufferBits <= self.len else self.len - self.lwi
while (self.ndx < buflmt) and ((self.cmpsts[self.ndx >> 3] >> (self.ndx & 7)) & 1 != 0):
self.ndx += 1
if self.lwi + self.ndx >= self.len: self.sz = 0
return rslt
fn main():
print("The primes to 100 are:")
for prm in SoEOddsPaged(100): print_no_newline(prm, " ")
print()
let strt0 = now()
let answr0 = SoEOddsPaged(1_000_000).countPrimes()
let elpsd0 = (now() - strt0) / 1000000
print("Found", answr0, "primes up to 1,000,000 in", elpsd0, "milliseconds.")
let strt1 = now()
let answr1 = SoEOddsPaged(cLIMIT).countPrimes()
let elpsd1 = (now() - strt1) / 1000000
print("Found", answr1, "primes up to", cLIMIT, "in", elpsd1, "milliseconds.")</syntaxhighlight>
{{out}}
<pre>The primes to 100 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
Found 78498 primes up to 1,000,000 in 0.084122000000000002 milliseconds.
Found 50847534 primes up to 1000000000 in 139.509275 milliseconds.</pre>
This was tested on the same computer as the previous Mojo versions. Note that the time now scales quite well with range since there are no longer the huge RAM access time bottleneck's. This version is only about 2.25 times slower than Kim Walich's primesieve program written in C++ and the mostly constant factor difference will be made up if one adds wheel factorization to the same level as he uses (basic wheel factorization ratio of 48/105 plus some other more minor optimizations). This version can count the number of primes to 1e11 in about 21.85 seconds on this machine. It will work reasonably efficiently up to a range of about 1e14 before other optimization techniques such as "bucket sieving" should be used.
For counting the number of primes to a billion (1e9), this version has reduced the time by about a factor of 40 from the original version and over eight times from the odds-only version above. Adding wheel factorization will make it almost two and a half times faster yet for a gain in speed of about a hundred times over the original version.
=={{header|MUMPS}}==
Line 13,590 ⟶ 14,725:
</pre>
=={{header|PascalABC.NET}}==
<syntaxhighlight lang="delphi">
function Eratosthenes(N: integer): List<integer>;
type primetype = (nonprime,prime);
begin
var sieve := |nonprime|*2 + |prime|*(N-1);
for var i:=2 to N.Sqrt.Round do
if sieve[i] = prime then
for var j := i*i to N step i do
sieve[j] := nonprime;
Result := new List<integer>;
for var i:=2 to N do
if sieve[i] = prime then
Result.Add(i);
end;
begin
Eratosthenes(1000).Println
end.
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
=={{header|Perl}}==
Line 15,663 ⟶ 16,824:
== [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]
</syntaxhighlight>
=={{header|Refal}}==
<syntaxhighlight lang="refal">$ENTRY Go {
= <Print <Primes 100>>;
};
Primes {
s.N = <Sieve <Iota 2 s.N>>;
};
Iota {
s.End s.End = s.End;
s.Start s.End = s.Start <Iota <+ 1 s.Start> s.End>;
};
Cross {
s.Step e.List = <Cross (s.Step 1) s.Step e.List>;
(s.Step s.Skip) = ;
(s.Step 1) s.Item e.List = X <Cross (s.Step s.Step) e.List>;
(s.Step s.N) s.Item e.List = s.Item <Cross (s.Step <- s.N 1>) e.List>;
};
Sieve {
= ;
X e.List = <Sieve e.List>;
s.N e.List = s.N <Sieve <Cross s.N e.List>>;
};</syntaxhighlight>
{{out}}
<pre>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</pre>
=={{header|REXX}}==
Line 15,908 ⟶ 17,098:
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
</pre>
=={{header|RPL}}==
This is a direct translation from Wikipedia. The variable <code>i</code> has been renamed <code>ii</code> to avoid confusion with the language constant <code>i</code>=√ -1
{{works with|Halcyon Calc|4.2.8}}
{| class="wikitable"
! RPL code
! Comment
|-
|
≪ → n
≪ { } n + 1 CON 'A' STO
2 n √ '''FOR''' ii
'''IF''' A ii GET '''THEN'''
ii SQ n '''FOR''' j
'A' j 0 PUT ii '''STEP'''
'''END'''
'''NEXT'''
{ }
2 n '''FOR''' ii '''IF''' A ii GET '''THEN''' ii + '''END NEXT'''
'A' PURGE
≫ ≫ '<span style="color:blue">SIEVE</span>' STO
|
<span style="color:blue">SIEVE</span> ''( n -- { prime_numbers } )''
let A be an array of Boolean values, indexed by 2 to n,
initially all set to true.
for i = 2, 3, 4, ..., not exceeding √n do
if A[i] is true
for j = i^2, i^2+i,... not exceeding n do
set A[j] := false
return all i such that A[i] is true.
|}
100 <span style="color:blue">SIEVE</span>
{{out}}
<pre>
1: { 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 }
</pre>
Latest RPL versions allow to remove some slow <code>FOR..NEXT</code> loops and use local variables only.
{{works with|HP|49}}
« 'X' DUP 1 4 PICK 1 SEQ DUP → n a seq123
« 2 n √ '''FOR''' ii
'''IF''' a ii GET '''THEN'''
ii SQ n '''FOR''' j
'a' j 0 PUT ii '''STEP'''
'''END'''
'''NEXT'''
a seq123 IFT TAIL
» » '<span style="color:blue">SIEVE</span>' STO
{{works with|HP|49}}
=={{header|Ruby}}==
Line 17,190 ⟶ 18,434:
Original source: [http://seed7.sourceforge.net/algorith/math.htm#sieve_of_eratosthenes]
=={{header|SETL}}==
<syntaxhighlight lang="setl">program eratosthenes;
print(sieve 100);
op sieve(n);
numbers := [1..n];
numbers(1) := om;
loop for i in [2..floor sqrt n] do
loop for j in [i*i, i*i+i..n] do
numbers(j) := om;
end loop;
end loop;
return [n : n in numbers | n /= om];
end op;
end program;</syntaxhighlight>
{{out}}
<pre>[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]</pre>
=={{header|Sidef}}==
Line 17,477 ⟶ 18,739:
Output:
<pre>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</pre>
=={{header|SparForte}}==
As a structured script.
<syntaxhighlight lang="ada">#!/usr/local/bin/spar
pragma annotate( summary, "sieve" );
pragma annotate( description, "The Sieve of Eratosthenes is a simple algorithm that" );
pragma annotate( description, "finds the prime numbers up to a given integer. Implement ");
pragma annotate( description, "this algorithm, with the only allowed optimization that" );
pragma annotate( description, "the outer loop can stop at the square root of the limit," );
pragma annotate( description, "and the inner loop may start at the square of the prime" );
pragma annotate( description, "just found. That means especially that you shouldn't" );
pragma annotate( description, "optimize by using pre-computed wheels, i.e. don't assume" );
pragma annotate( description, "you need only to cross out odd numbers (wheel based on" );
pragma annotate( description, "2), numbers equal to 1 or 5 modulo 6 (wheel based on 2" );
pragma annotate( description, "and 3), or similar wheels based on low primes." );
pragma annotate( see_also, "http://rosettacode.org/wiki/Sieve_of_Eratosthenes" );
pragma annotate( author, "Ken O. Burtch" );
pragma license( unrestricted );
pragma restriction( no_external_commands );
procedure sieve is
last_bool : constant positive := 20;
type bool_array is array(2..last_bool) of boolean;
a : bool_array;
test_num : positive;
-- limit : positive := positive(numerics.sqrt(float(arrays.last(a))));
-- n : positive := 2;
begin
for i in arrays.first(a)..last_bool loop
a(i) := true;
end loop;
for num in arrays.first(a)..last_bool loop
if a(num) then
test_num := num * num;
while test_num <= last_bool loop
a(test_num) := false;
test_num := @ + num;
end loop;
end if;
end loop;
for i in arrays.first(a)..last_bool loop
if a(i) then
put_line(i);
end if;
end loop;
end sieve;</syntaxhighlight>
=={{header|Standard ML}}==
Line 18,993 ⟶ 20,306:
=={{header|Wren}}==
<syntaxhighlight lang="
if (n < 2) return []
var comp = List.filled(n-1, false)
|