Rosetta Code

Sieve of Eratosthenes: Difference between revisions

← Older edit
Content deleted Content added
→‎{{header|Mojo}}: updated to recent 24.4.0 version...
 
(35 intermediate revisions by 13 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 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 5,764 ⟶ 5,797:
 
=={{header|EasyLang}}==
<syntaxhighlight lang="text">
len is_divisiblesieve[] 100
max = sqrt len is_divisiblesieve[]
for di = 2 to max
if is_divisiblesieve[di] = 0
for ij = di * d step d to len is_divisible[]i
while j <= is_divisiblelen sieve[i] = 1
sieve[j] = 1
j += i
.
.
.
for i = 2 to len is_divisiblesieve[]
if is_divisiblesieve[i] = 0
print i
.
.
.</syntaxhighlight>
</syntaxhighlight>
 
=={{header|eC}}==
Line 6,418 ⟶ 6,454:
=={{header|Elm}}==
 
===Elm with immutable arrays===
The Elm language doesn't handle efficiently doing the Sieve of Eratosthenes (SoE) algorithm efficiently because it doesn't have directly accessible linear 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:
<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,426 ⟶ 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,468 ⟶ 6,773:
in CIS 2 <| \ () -> oddprms()
 
type alias Model = ( Int,List String, String )
 
type alias Msg = Start Int | Done ( Int, Int )Model
 
test : (() -> CIS Int) -> Int -> Cmd Msg
cLIMIT : Int
test primesf lmt =
cLIMIT = 1000000
let
 
to100 = primesf() |> uptoCIS2List 100
timemillis : () -> Task Never Int
|> List.map String.fromInt |> String.join ", "
timemillis() = now |> andThen (\ t -> succeed (posixToMillis t))
to100str = "The primes to 100 are: " ++ to100
 
timemillis() = now |> andThen (succeed << posixToMillis)
test : Int -> Cmd Msg
in timemillis() |> andThen (\ strt ->
test lmt =
let cnt = countCISTo lmt primesf(primesTreeFolding()) --|> (soeDict())countCISTo lmt
in timemillis() |> andThen (\ tstop -> succeed ( t, cnt )) |> perform Done
let answrstr = "Found " ++ (String.fromInt cnt) ++ " primes to "
 
++ (String.fromInt cLIMIT) ++ " in "
myupdate : Msg -> Model -> (Model, Cmd Msg)
++ (String.fromInt (stop - strt)) ++ " milliseconds."
myupdate msg mdl =
in succeed [to100str, answrstr] ) ) |> perform identity
let (strttm, oldstr, _) = mdl in
case msg of
Start tm -> ( ( tm, oldstr, "Running..." ), test cLIMIT )
Done (stptm, answr) ->
( ( stptm, oldstr, "Found " ++ String.fromInt answr ++
" primes to " ++ String.fromInt cLIMIT ++
" in " ++ String.fromInt (stptm - strttm) ++ " milliseconds." )
, Cmd.none )
 
myview : Model -> Html msg
myview mdl =
let (_, str1, str2) = mdl
in div [] [ div [] [text str1], div [] [text str2] ]
 
main : Program () Model Msg
main =
element { init = \ _ -> ( ([], 0test primesTreeFolding cLIMIT )
, update = \ msg _ -> (msg, "The primes up to 100 are: " ++Cmd.none)
(primesTreeFolding() |> uptoCIS2List 100
|> List.map String.fromInt
|> String.join ", ") ++ "."
, "" ), perform Start (timemillis()) )
, update = myupdate
, subscriptions = \ _ -> Sub.none
, view = myviewdiv [] << 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 295201 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).
 
===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">moduletype MainPriorityQ exposingcomparable (v main )=
 
import Task exposing ( Task, succeed, perform, andThen )
import Html exposing (Html, div, text)
import Browser exposing ( element )
import Time exposing ( now, posixToMillis )
 
cLIMIT : Int
cLIMIT = 1000000
 
-- an infinite non-empty non-memoizing Co-Inductive Stream (CIS)...
type CIS a = CIS a (() -> CIS a)
 
uptoCIS2List : comparable -> CIS comparable -> List comparable
uptoCIS2List n cis =
let loop (CIS hd tl) lst =
if hd > n then List.reverse lst
else loop (tl()) (hd :: lst)
in loop cis []
 
countCISTo : comparable -> CIS comparable -> Int
countCISTo lmt cis =
let cntem acc (CIS hd tlf) =
if hd > lmt then acc else cntem (acc + 1) (tlf())
in cntem 0 cis
 
type PriorityQ comparable v =
Mt
| Br comparable v (PriorityQ comparable v)
Line 6,602 ⟶ 6,866:
else CIS n <| \ () -> sieve (n + 2) pq q bps
oddprms() = CIS 3 <| \ () -> sieve 5 emptyPQ 9 <| oddprms()
in CIS 2 <| \ () -> oddprms()</syntaxhighlight>
 
type alias Model = ( Int, String, String )
 
type Msg = Start Int | Done ( Int, Int )
 
timemillis : () -> Task Never Int
timemillis() = now |> andThen (\ t -> succeed (posixToMillis t))
 
test : Int -> Cmd Msg
test lmt =
let cnt = countCISTo lmt (primesPQ())
in timemillis() |> andThen (\ t -> succeed ( t, cnt )) |> perform Done
 
myupdate : Msg -> Model -> (Model, Cmd Msg)
myupdate msg mdl =
let (strttm, oldstr, _) = mdl in
case msg of
Start tm -> ( ( tm, oldstr, "Running..." ), test cLIMIT )
Done (stptm, answr) ->
( ( stptm, oldstr, "Found " ++ String.fromInt answr ++
" primes to " ++ String.fromInt cLIMIT ++
" in " ++ String.fromInt (stptm - strttm) ++ " milliseconds." )
, Cmd.none )
 
myview : Model -> Html msg
myview mdl =
let (_, str1, str2) = mdl
in div [] [ div [] [text str1], div [] [text str2] ]
 
main : Program () Model Msg
main =
element { init = \ _ -> ( ( 0
, "The primes up to 100 are: " ++
(primesPQ() |> uptoCIS2List 100
|> List.map String.fromInt
|> String.join ", ") ++ "."
, "" ), perform Start (timemillis()) )
, update = myupdate
, subscriptions = \ _ -> Sub.none
, view = myview }</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 192124 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).
 
==={{header|Elm Sieve of Eratosthenes function with range of primes}}===
{{incorrect|Bash|This version uses remainder testing (modBy) and so is a trial division algorithm, not a sieve of Eratosthenes.}}
Note: The David Turner algorithm is incorrect as per [https://en.m.wikipedia.org/wiki/Sieve_of_Eratosthenes Sieve of Eratosthenes article on Wikipedia]
 
The range of primes output is defined by the lower and upper limit.
This module exposes the function sieve to be imported by the module Main.
The test case is performed in elm repl.
 
<syntaxhighlight lang="elm">module Sieve exposing (sieve)
 
sieve : Int -> Int -> List Int
sieve min limit =
let
numbers : List Int
numbers =
List.range 2 limit
 
last : Int
last =
limit
|> toFloat
|> sqrt
|> round
 
isMultiple : Int -> Int -> Bool
isMultiple n m =
n /= m && modBy m n == 0
in
List.range 2 last
|> List.foldl
(\current result ->
List.filter
(\elem -> not (isMultiple elem current))
result
)
numbers
 
|> List.filter (\el -> el >= min)
 
</syntaxhighlight>
{{out}}
<pre>
% elm repl
---- Elm 0.19.1 ----------------------------------------------------------------
Say :help for help and :exit to exit! More at <https://elm-lang.org/0.19.1/repl>
--------------------------------------------------------------------------------
> import Sieve exposing (sieve)
> a= sieve 1000200 1000600
[1000211,1000213,1000231,1000249,1000253,1000273,1000289,1000291,1000303,1000313,1000333,1000357,1000367,1000381,1000393,1000397,1000403,1000409,1000423,1000427,1000429,1000453,1000457,1000507,1000537,1000541,1000547,1000577,1000579,1000589]
: List Int
>
</pre>
 
=={{header|Emacs Lisp}}==
Line 6,951 ⟶ 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 11,571 ⟶ 11,778:
=={{header|langur}}==
{{trans|D}}
<syntaxhighlight lang="langur">val .sieve = f(.limit) {
val sieve = if .fn(limit < 2) {
if limit < 2: return []
}
 
var .composite = arr .limit, * [false]
.composite[1] = true
 
for .n in 2 to.. truncatetrunc(.limit ^/ 2) + 1 {
if not .composite[.n] {
for .k = .n^2 ; .k < .limit ; .k += .n {
.composite[.k] = true
}
}
}
 
filter f(.fn n) :not .composite[.n], series .(limit-1)
}
 
writeln .sieve(100)</syntaxhighlight>
</syntaxhighlight>
 
{{out}}
Line 12,466 ⟶ 12,673:
IO.Put("\n");
END Sieve.</syntaxhighlight>
 
=={{header|Mojo}}==
 
Tested with Mojo version 24.4.0:
 
<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):
var s = i * (i + 4) + 2
if s >= limit - 1: break
if self.cmpsts[i]: continue
var 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
var 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(prm, " ",end="")
print()
var strt0 = now()
var answr0 = len(SoEBasic(1_000_000))
var elpsd0 = (now() - strt0) / 1000000
print("Found", answr0, "primes up to 1,000,000 in", elpsd0, "milliseconds.")
var strt1 = now()
var answr1 = len(SoEBasic(cLIMIT))
var 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 iterator struct which seems to be the Mojo-preferred way to do this, and normally a might have been used the the culling array but here the raw unsafe 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 composite 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 optimized by modern compilers to use efficient auto-vectorization and the use of SIMD registers available to the architecture to reduce these culling operations to an average of a tiny fraction of a CPU clock cycle per cull.
 
Mojo version 24.4.0 was tested:
 
<syntaxhighlight lang="mojo">from memory import (memset_zero, memcpy)
from memory.unsafe import (DTypePointer)
from bit import pop_count
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)
var r1: Int = ((s + bp) >> 3) - (s >> 3)
var r2: Int = ((s + 2 * bp) >> 3) - (s >> 3)
var r3: Int = ((s + 3 * bp) >> 3) - (s >> 3)
var r4: Int = ((s + 4 * bp) >> 3) - (s >> 3)
var r5: Int = ((s + 5 * bp) >> 3) - (s >> 3)
var r6: Int = ((s + 6 * bp) >> 3) - (s >> 3)
var r7: Int = ((s + 7 * bp) >> 3) - (s >> 3)
var 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
var 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 mkExtremeFuncs[SIZE: Int]() -> UnsafePointer[UnrollFunc, 0]:
var jmptbl = UnsafePointer[UnrollFunc, 0].alloc(SIZE)
@parameter
for i in range(SIZE):
alias OFST = i >> 2
alias BP = ((i & 3) << 1) + 1
jmptbl[i] = extreme[OFST, BP]
return jmptbl
alias DenseFunc = fn(DTypePointer[DType.uint64], Int, Int) -> DTypePointer[DType.uint64]
 
@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)
var plmt = pcmps + (bufsz >> 3) - BP
while cp < plmt:
@parameter
for n in range(64):
alias MUL = n * BP
var cop = cp.offset(MUL >> 6)
cop.store(cop.load() | (1 << (MUL & 63)))
cp += BP
return cp
 
fn mkDenseFuncs[SIZE: Int]() -> UnsafePointer[DenseFunc, 0]:
var jmptbl = UnsafePointer[DenseFunc, 0].alloc(SIZE)
@parameter
for i in range(SIZE):
alias BP = (i << 1) + 3
jmptbl[i] = denseCullFunc[BP]
return jmptbl
 
@always_inline
fn cullPass(dfs: UnsafePointer[DenseFunc, 0], efs: UnsafePointer[UnrollFunc, 0],
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
var bcp = dfs[(bp - 3) >> 1](cmpsts.bitcast[DType.uint64](), bytesz, sm)
var ns = 0
var ncp = bcp
var 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
efs[((s & 7) << 2) + ((bp & 7) >> 1)](cmpsts, bytesz, s, bp)
# for c in range(s, bytesz * 8, bp): # slow bit twiddling way
# cmpsts[c >> 3] |= (1 << (c & 7))
 
fn countPagePrimes(ptr: DTypePointer[DType.uint8], bitsz: Int) -> Int:
var wordsz: Int = (bitsz + 63) // 64 # round up to nearest 64 bit boundary
var rslt: Int = wordsz * 64
var bigcmps = ptr.bitcast[DType.uint64]()
for i in range(wordsz - 1):
rslt -= int(pop_count(bigcmps[i]))
rslt -= int(pop_count(bigcmps[wordsz - 1] | (-2 << ((bitsz - 1) & 63))))
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
var 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)
var denseFuncs : UnsafePointer[DenseFunc, 0] = mkDenseFuncs[64]()
var extremeFuncs: UnsafePointer[UnrollFunc, 0] = mkExtremeFuncs[32]()
for i in range(self.len):
var s = (i + i) * (i + 3) + 3
if s >= self.len: break
if (self.cmpsts[i >> 3] >> (i & 7)) & 1 != 0: continue
var bp = i + i + 3
cullPass(denseFuncs, extremeFuncs, 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
var 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
var 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(prm, " ", end="")
print()
var strt0 = now()
var answr0 = len(SoEOdds(1_000_000))
var elpsd0 = (now() - strt0) / 1000000
print("Found", answr0, "primes up to 1,000,000 in", elpsd0, "milliseconds.")
var strt1 = now()
var answr1 = len(SoEOdds(cLIMIT))
var 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 24.4.0:
 
<syntaxhighlight lang="mojo">from memory import (memset_zero, memcpy)
from memory.unsafe import (DTypePointer)
from bit import pop_count
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
var 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)
var r1: Int = ((s + bp) >> 3) - (s >> 3)
var r2: Int = ((s + 2 * bp) >> 3) - (s >> 3)
var r3: Int = ((s + 3 * bp) >> 3) - (s >> 3)
var r4: Int = ((s + 4 * bp) >> 3) - (s >> 3)
var r5: Int = ((s + 5 * bp) >> 3) - (s >> 3)
var r6: Int = ((s + 6 * bp) >> 3) - (s >> 3)
var r7: Int = ((s + 7 * bp) >> 3) - (s >> 3)
var 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
var 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 mkExtremeFuncs[SIZE: Int]() -> UnsafePointer[UnrollFunc, 0]:
var jmptbl = UnsafePointer[UnrollFunc, 0].alloc(SIZE)
@parameter
for i in range(SIZE):
alias OFST = i >> 2
alias BP = ((i & 3) << 1) + 1
jmptbl[i] = extreme[OFST, BP]
return jmptbl
alias DenseFunc = fn(DTypePointer[DType.uint64], Int, Int) -> DTypePointer[DType.uint64]
 
@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)
var plmt = pcmps + (bufsz >> 3) - BP
while cp < plmt:
@parameter
for n in range(64):
alias MUL = n * BP
var cop = cp.offset(MUL >> 6)
cop.store(cop.load() | (1 << (MUL & 63)))
cp += BP
return cp
 
fn mkDenseFuncs[SIZE: Int]() -> UnsafePointer[DenseFunc, 0]:
var jmptbl = UnsafePointer[DenseFunc, 0].alloc(SIZE)
@parameter
for i in range(SIZE):
alias BP = (i << 1) + 3
jmptbl[i] = denseCullFunc[BP]
return jmptbl
 
@always_inline
fn cullPass(dfs: UnsafePointer[DenseFunc, 0], efs: UnsafePointer[UnrollFunc, 0],
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
var bcp = dfs[(bp - 3) >> 1](cmpsts.bitcast[DType.uint64](), bytesz, sm)
var ns = 0
var ncp = bcp
var 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
efs[((s & 7) << 2) + ((bp & 7) >> 1)](cmpsts, bytesz, s, bp)
# for c in range(s, bytesz * 8, bp): # slow bit twiddling way
# cmpsts[c >> 3] |= (1 << (c & 7))
 
fn cullPage(dfs: UnsafePointer[DenseFunc, 0], efs: UnsafePointer[UnrollFunc, 0],
lwi: Int, lmt: Int, cmpsts: DTypePointer[DType.uint8], bsprmrps: DTypePointer[DType.uint8]):
var bp = 1; var ndx = 0
while True:
bp += int(bsprmrps[ndx]) << 1
var 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(dfs, efs, cmpsts, cBufferSize, s, bp)
ndx += 1
 
fn countPagePrimes(ptr: DTypePointer[DType.uint8], bitsz: Int) -> Int:
var wordsz: Int = (bitsz + 63) // 64 # round up to nearest 64 bit boundary
var rslt: Int = wordsz * 64
var bigcmps = ptr.bitcast[DType.uint64]()
for i in range(wordsz - 1):
rslt -= int(pop_count(bigcmps[i]))
rslt -= int(pop_count(bigcmps[wordsz - 1] | (-2 << ((bitsz - 1) & 63))))
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
var 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)
var denseFuncs : UnsafePointer[DenseFunc, 0] = mkDenseFuncs[64]()
var extremeFuncs: UnsafePointer[UnrollFunc, 0] = mkExtremeFuncs[32]()
for i in range(self.len):
var s = (i + i) * (i + 3) + 3
if s >= self.len: break
if (self.cmpsts[i >> 3] >> (i & 7)) & 1 != 0: continue
var bp = i + i + 3
cullPass(denseFuncs, extremeFuncs, 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
var 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
var rslt = (self.ndx << 1) + 3; self.sz -= 1; self.ndx += 1; return rslt
 
struct SoEOddsPaged:
var denseFuncs : UnsafePointer[DenseFunc, 0]
var extremeFuncs: UnsafePointer[UnrollFunc, 0]
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.denseFuncs = mkDenseFuncs[64]()
self.extremeFuncs = mkExtremeFuncs[32]()
self.len = 0 if limit < 2 else int(((limit - 3) // 2 + 1))
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:
var 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)
var bsprmitr = SoEOdds(int(intsqrt(limit)))
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.denseFuncs = existing.denseFuncs
self.extremeFuncs = existing.extremeFuncs
self.len = existing.len
self.sz = existing.sz
var 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.denseFuncs = existing.denseFuncs
self.extremeFuncs = existing.extremeFuncs
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
var cmpsts = DTypePointer[DType.uint8].alloc(cBufferSize)
memset_zero(cmpsts, cBufferSize)
cullPage(self.denseFuncs, self.extremeFuncs, 0, cBufferBits, cmpsts, self.bsprmrps)
while lwi + cBufferBits <= self.len:
cnt += countPagePrimes(cmpsts, cBufferBits)
lwi += cBufferBits
memset_zero(cmpsts, cBufferSize)
var lmt = lwi + cBufferBits if lwi + cBufferBits <= self.len else self.len
cullPage(self.denseFuncs, self.extremeFuncs, 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:
var bytesz = ((self.len + 63) & -64) >> 3 # round up to nearest 64 bit boundary
memset_zero(self.cmpsts, bytesz)
for i in range(self.len):
var s = (i + i) * (i + 3) + 3
if s >= self.len: break
if (self.cmpsts[i >> 3] >> (i & 7)) & 1 != 0: continue
var bp = i + i + 3
cullPass(self.denseFuncs, self.extremeFuncs, self.cmpsts, bytesz, s, bp)
else:
memset_zero(self.cmpsts, cBufferSize)
cullPage(self.denseFuncs, self.extremeFuncs, 0, cBufferBits, self.cmpsts, self.bsprmrps)
return 2
var 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)
var lmt = self.lwi + cBufferBits if self.lwi + cBufferBits <= self.len else self.len
cullPage(self.denseFuncs, self.extremeFuncs, self.lwi, lmt, self.cmpsts, self.bsprmrps)
var 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(prm, " ", end="")
print()
var strt0 = now()
var answr0 = SoEOddsPaged(1_000_000).countPrimes()
var elpsd0 = (now() - strt0) / 1000000
print("Found", answr0, "primes up to 1,000,000 in", elpsd0, "milliseconds.")
var strt1 = now()
var answr1 = SoEOddsPaged(cLIMIT).countPrimes()
var 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++ when run single-threaded 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 100 million 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,441 ⟶ 14,258:
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]</syntaxhighlight>
 
===Yet Another Functional Version===
 
The problem with the above two "functional" versions using OCaml list's is that they aren't written to be Tail Call Recursive and will thus blow stack for more than trivial ranges; the following "incremental" (as in each new value depends on some previous state) "infinite series" (is a generator function that produces a continuous stream or actually a Co-Inductive Stream - CIS - of prime values) Sieve of Eratosthenes code based on adding "infinite tree folding" to a sieve by Richard Bird doesn't have the above limitations:
<syntaxhighlight lang="ocaml">let cLIMIT = 1_000_000
 
type 'a cis = CIS of 'a * (unit -> 'a cis)
 
let primesTF() =
let rec merge (CIS(x, xtl) as xs) (CIS(y, ytl) as ys) =
if x < y then CIS(x, fun() -> merge (xtl()) ys)
else if y < x then CIS(y, fun() -> merge xs (ytl()))
else CIS(x, fun() -> merge (xtl()) (ytl()))
in let bpmults bp =
let adv = bp + bp in
let rec pmlt vr = CIS(vr, fun() -> pmlt (vr + adv))
in pmlt (bp * bp)
in let rec allmlts = function
| CIS(bp, bptf) -> CIS(bpmults bp, fun() -> allmlts (bptf()))
in let rec pairs = function
| CIS(mcs, cstf) ->
let CIS(nmcs, ncstf) = cstf() in
CIS(merge mcs nmcs, fun() -> pairs (ncstf()))
in let rec cmpsts = function
| CIS(CIS(cr, ctfr), cstf) ->
CIS(cr, fun() -> merge (ctfr()) (cstf() |> pairs |> cmpsts))
in let rec testAt n csr =
match csr with
| CIS(cr, ctlr) ->
if n >= cr then testAt (n + 2) (ctlr())
else CIS(n, fun() -> testAt (n + 2) csr)
in let rec oddprms() = CIS(3, fun() -> oddprms() |> allmlts |> cmpsts |> testAt 5)
in CIS(2, fun() -> oddprms())
 
let showprmsTo lmt pcis =
let rec loop lst = function
| CIS(p, ptf) -> if p > lmt then List.rev lst |> List.map string_of_int
|> String.concat "; "|> print_endline
else loop (p :: lst) (ptf())
in loop [] pcis
 
let countprmsTo lmt pcis =
let rec loop cnt = function
| CIS(p, ptf) -> if p > lmt then cnt else loop (cnt + 1) (ptf())
in loop 0 pcis
 
let _ = showprmsTo 100 (primesTF())
let strt = Sys.time()
let answr = countprmsTo cLIMIT (primesTF())
let elpsd = (Sys.time() -. strt) *. 1000.
let _ = Printf.printf "Found %d primes to %d in %f milliseconds.\r\n" answr cLIMIT elpsd</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
Found 78498 primes to 1000000 in 71.197000 milliseconds.</pre>
 
This is actually an odds-only sieve since two is the only even prime, and works recursively in feeding a secondary version of the resulting odd primes back into the sieving mechanism to avoid the size of the "infinite folded tree" of merged base prime multiples from being as large as the range currently sieved but rather to only the square root of the current range. It takes the base primes stream and generates an infinite (as required) stream containing the streams (a stream of streams) of the multiples of each of the base primes generated by `allmults` and `bpmults, then merges this stream of streams using the `pairs` function (thereby doing the "infinite tree folding") and the `merge` function into a single stream of odd composite numbers (`cmpsts`). This works as written because the first value (the square of the base prime) is always less than the first value of the next multiples stream (the next base prime's square). Finally, the `testAt` function compares all of the odd numbers with the stream of composites (only needing to check the first value of the composites stream being consumed) and if not included, adds the current odd test value to the output stream, skipping over any numbers that are found in the composites stream. The odd primes stream as generated by the function `oddprms` has a first value set as three and the testing starts at five in order to avoid the data race due to the recursion so that there are always composite values starting at nine in place before the `testAt` comparisons start with five.
 
This algorithm using CIS's does have the overhead of needing to generate closures (functions that capture external values that are not function parameters) for every step of new prime generated, but this is quite well optimized for OCaml.
 
This function is of O(n (log n) log log n) complexity due to the (binary) folded merge tree comparisons, but the OCaml List versions above don't have the ideal O(n log log n) complexity either since they need to scan the lists from the start for every culling pass, as well as being limited by generating stack overflows...
 
=== Yet another (mostly) functional version ===
 
This sieve is just "almost"/mostly a purely functional algorithm because it still mutates the contents of the sieving array.
 
<syntaxhighlight lang="ocaml">
let sieve limit =
let p = Array.make (limit + 1) true in
let rec sieve_outer d =
if d * d > limit then p
else if p.(d) then
let rec sieve_inner m =
if m > limit then sieve_outer (d + 1)
else ((p.(m) <- false); sieve_inner (m + d))
in sieve_inner (d * d)
else sieve_outer (d + 1)
in sieve_outer 2
 
let primes limit =
let s = (sieve limit) in
List.init (limit - 1) (fun i -> i + 2)
|> List.filter (fun x -> s.(x))</syntaxhighlight>
 
In the top-level:
<syntaxhighlight lang="ocaml"># 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]
</syntaxhighlight>
 
=={{header|Oforth}}==
Line 13,881 ⟶ 14,787:
 
</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,954 ⟶ 16,886:
== [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 16,213 ⟶ 17,174:
2 n √ '''FOR''' ii
'''IF''' A ii GET '''THEN'''
ii 2 ^SQ n '''FOR''' j
'A' j 0 PUT ii '''STEP'''
'''END NEXT'''
{ } 2 n '''FOR''' ii '''IF''' A ii GET '''THEN''' ii + '''END NEXT'''
{ }
≫ ≫ ‘'''SIEVE'''’ STO
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.
Line 16,226 ⟶ 17,190:
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 19,389 ⟶ 20,368:
 
=={{header|Wren}}==
<syntaxhighlight lang="ecmascriptwren">var sieveOfE = Fn.new { |n|
if (n < 2) return []
var comp = List.filled(n-1, false)
Line 19,501 ⟶ 20,480:
print "time: ", t, "\n"
print peek("millisrunning")</syntaxhighlight>
 
=={{header|YAMLScript}}==
<syntaxhighlight lang="yaml">
!yamlscript/v0
 
defn main(n):
primes =: primes-up-to(n)
say: |
The $(primes.count()) prime numbers less than $n are:
$(primes.join("\n"))
 
defn primes-up-to(limit):
"Returns a lazy sequence of prime numbers less than limit":
max-i =: int((limit - 1) / 2)
refs =: boolean-array(max-i true)
root =: dec(int(Math/sqrt(limit))) / 2
doseq [i (1 .. root) :when aget(refs i)]:
doseq [j range((2 * i * inc(i)) max-i (i + i + 1))]:
aset refs: j false
cons 2:
map \(+ %1 %1 1):
filter \(aget refs %1):
range: 1 max-i
</syntaxhighlight>
 
=={{header|Zig}}==
Retrieved from "https://rosettacode.org/wiki/Sieve_of_Eratosthenes"