Rosetta Code

Sieve of Eratosthenes: Difference between revisions

Page Discussion

Sieve of Eratosthenes (view source)

Revision as of 19:34, 1 May 2024

38,112 bytes added ,  17 days ago
→‎{{header|langur}}
 
(29 intermediate revisions by 15 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 4,702 ⟶ 4,735:
endif
 
next n</syntaxhighlight>
{{out| Output}}<pre>
 
prime numbers less than or equal to 120 are:
end</syntaxhighlight>
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,523 ⟶ 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 6,407 ⟶ 6,451:
=={{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,415 ⟶ 6,725:
import Html exposing (Html, div, text)
import Time exposing (now, posixToMillis)
 
cLIMIT : Int
cLIMIT = 1000000
 
type CIS a = CIS a (() -> CIS a)
Line 6,457 ⟶ 6,770:
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,591 ⟶ 6,863:
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|Emacs Lisp}}==
Line 8,256 ⟶ 8,490:
I: DWord;
begin
for I := 0 toif aList.Size -<> 0 2then dobegin
Write(if aList[I],.Size > ',1 ');then
WriteLn(aList[ for I := 0 to aList.Size - 1]);2 do
Write(aList[I], ', ');
WriteLn(aList[aList.Size - 1]);
end;
aList.Free;
end;
Line 8,397 ⟶ 8,634:
for I := 0 to Length(aList) - 2 do
Write(aList[I], ', ');
WriteLn(if aList[High(aList)]); <> nil then
WriteLn(aList[High(aList)]);
end;
begin
Line 8,561 ⟶ 8,799:
=={{header|Fōrmulæ}}==
 
{{FormulaeEntry|page=https://formulae.org/?script=examples/Sieve_of_Eratosthenes}}
Fōrmulæ programs are not textual, visualization/edition of programs is done showing/manipulating structures but not text. Moreover, there can be multiple visual representations of the same program. Even though it is possible to have textual representation &mdash;i.e. XML, JSON&mdash; they are intended for storage and transfer purposes more than visualization and edition.
 
'''Solution'''
Programs in Fōrmulæ are created/edited online in its [https://formulae.org website], However they run on execution servers. By default remote servers are used, but they are limited in memory and processing power, since they are intended for demonstration and casual use. A local server can be downloaded and installed, it has no limitations (it runs in your own computer). Because of that, example programs can be fully visualized and edited, but some of them will not run if they require a moderate or heavy computation/memory resources, and no local server is being used.
 
[[File:Fōrmulæ - Sieve of Eratosthenes 01.png]]
In '''[https://formulae.org/?example=Sieve_of_Eratosthenes this]''' page you can see the program(s) related to this task and their results.
 
'''Test case'''
 
[[File:Fōrmulæ - Sieve of Eratosthenes 02.png]]
 
[[File:Fōrmulæ - Sieve of Eratosthenes 03.png]]
 
=={{header|GAP}}==
Line 8,964 ⟶ 9,208:
 
var p int
p = <-primes
p = <-primes
 
Line 10,539 ⟶ 10,783:
Alternate version using <code>findall</code> to get all primes at once in the end
 
<syntaxhighlight lang="julia">function sieve(n :: IntInteger)
isprimeprimes = truesfill(true, n)
isprimeprimes[1] = false
for p in 2:n
if isprimeprimes[p] || continue
primes[p .* (2:n÷p)] j .= p * pfalse
if j > n
return findall(isprime)
else
for k in j:p:n
isprime[k] = false
end
end
end
end
findall(primes)
end</syntaxhighlight>
 
Line 11,503 ⟶ 11,740:
=={{header|langur}}==
{{trans|D}}
<syntaxhighlight lang="langur">val .sieve = ffn(.limit) {
if .limit < 2: {return []
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
}
Line 11,519 ⟶ 11,754:
}
 
filter ffn(.n) { not .composite[.n] }, series .limit-1
}
 
writeln .sieve(100)</syntaxhighlight>
</syntaxhighlight>
 
{{out}}
Line 12,398 ⟶ 12,634:
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 15,886 ⟶ 16,759:
== [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,145 ⟶ 17,047:
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,158 ⟶ 17,063:
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,449 ⟶ 18,369:
 
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 19,303 ⟶ 20,241:
 
=={{header|Wren}}==
<syntaxhighlight lang="ecmascriptwren">var sieveOfE = Fn.new { |n|
if (n < 2) return []
var comp = List.filled(n-1, false)
Langurmonkey
885

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/343386...364154"