Sieve of Pritchard
You are encouraged to solve this task according to the task description, using any language you may know.
The Sieve of Pritchard is a modern algorithm for finding prime numbers. It takes many fewer operations than the Sieve of Eratosthenes (better time complexity), at the cost of greater storage requirements (worse space complexity).
Conceptually, it works by constructing a series of "wheels" marked along their circumference with the pattern of primes up to the value of successive primorial numbers (where the Nth primorial is the product of the first N primes). Those wheels are then rolled along the number line, and only the numbers touched by the marks are considered as candidate primes, in contrast to Eratosthenes' sieve in which all the integers in the range start out as candidates. (The Sieve of Pritchard is an example of the "wheel-based optimizations" mentioned in the Eratosthenes task.)
For example, the second-order wheel has size 6 (the product of the first two primes, 2 and 3) and is marked only at the numbers between 1 and 6 that are not multiples of 2 or 3, namely 1 and 5. As this wheel is rolled along the number line, it will pick up only numbers of the form 6k+1 or 6k+5 (that is, n where n mod 6 is in {1,5}). By the time it stops at 30 (2x3x5) it has added only 8 of the numbers between 6 and 30 as candidates for primality, only one of which is actually composite and must be removed (25). In the process it has constructed the next wheel, which will add only nine out of every 30 numbers as it rolls up to 210.
This YouTube video tells a story to help motivate the algorithm's design;this one presents the execution of the algorithm for N=150 in a format that permits single-stepping forward and backward through the run. In that implementation, the list of primes is populated into a sparse global array s such that s[p] contains the next prime after p iff p is itself a prime in the target range; this allows numbers to be removed from consideration quickly without any the copying/shifting that would be required from a normally-packed array.
- Task
Write a program/subprogram that uses the Sieve of Pritchard algorithm to find all primes up to a specified limit. Show the result of running it with a limit of 150.
- Related tasks
- Sieve of Eratosthenes
- Emirp primes
- count in factors
- prime decomposition
- factors of an integer
- extensible prime generator
- primality by trial division
- factors of a Mersenne number
- trial factoring of a Mersenne number
- partition an integer X into N primes
- sequence of primes by Trial Division
AppleScript
on sieveOfPritchard(limit)
if (limit < 2) then return {}
script o
property primes : {}
property wheel : {1, 2}
property oldWheel : missing value
end script
set {oldCircumference, circumference} to {missing value, 2}
repeat until (oldCircumference = limit)
set o's oldWheel to o's wheel's numbers
set prime to o's oldWheel's second item
set end of o's primes to prime
set oldCircumference to circumference
set circumference to oldCircumference * prime
if (circumference > limit) then set circumference to limit
repeat with n from (oldCircumference + 1) to circumference
if (o's wheel's item ((n - 1) mod oldCircumference + 1) is missing value) then
set end of o's wheel to missing value
else
set end of o's wheel to n
end if
end repeat
repeat with this in o's oldWheel
set n to this * prime
if (n > circumference) then exit repeat
set o's wheel's item n to missing value
end repeat
end repeat
return o's primes & rest of o's wheel's numbers
end sieveOfPritchard
sieveOfPritchard(150)
- Output:
{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149}
While the above's fine for the current task, if it were needed to return primes up to the hundreds of thousands and beyond, it would be much faster to prefabricate the 'wheel' list to its final length by means of concatenation than to grow it on the fly by appending items.
on sieveOfPritchard(limit)
if (limit < 2) then return {}
script o
property primes : {}
property wheel : makeList(limit, missing value)
property oldWheel : missing value
end script
set {o's wheel's 1st item, o's wheel's 2nd item} to {1, 2}
set {oldCircumference, circumference} to {missing value, 2}
repeat until (oldCircumference = limit)
set o's oldWheel to o's wheel's numbers
set prime to o's oldWheel's second item
set end of o's primes to prime
set oldCircumference to circumference
set circumference to oldCircumference * prime
if (circumference > limit) then set circumference to limit
repeat with n from (oldCircumference + 1) to circumference
if (o's wheel's item ((n - 1) mod oldCircumference + 1) is not missing value) then
set o's wheel's item n to n
end if
end repeat
repeat with this in o's oldWheel
set n to this * prime
if (n > circumference) then exit repeat
set o's wheel's item n to missing value
end repeat
end repeat
return o's primes & rest of o's wheel's numbers
end sieveOfPritchard
on makeList(limit, filler)
if (limit < 1) then return {}
script o
property lst : {filler}
end script
set counter to 1
repeat until (counter + counter > limit)
set o's lst to o's lst & o's lst
set counter to counter + counter
end repeat
if (counter < limit) then set o's lst to o's lst & o's lst's items 1 thru (limit - counter)
return o's lst
end makeList
sieveOfPritchard(1000000)
C#
Loosely based on the Python version. I cut a couple of things out and it still worked. Not too crazy about having to create temporary lists to add or remove from the SortedSet, seems inefficient. But that is the work-around I employed, since SortedSets can't be accessed by indexing, and are non-mutable in a foreach loop. I haven't yet directly tested this against a Sieve of Eratosthenes to compare performance. The Wikipedia article suggests using a doubly linked list, so this C# incarnation is a kludge at best.
Compared to the prototype algorithm, it appears there isn't any code to do the follow-up end-of-wheel additions when necessary. But the main loop limit has been changed to go to the next prime, and the existing code handles the additions.
Updated to include "numbers added / removed (to / from members)" and performance statistics. The "removed" figure includes both composite numbers and prime numbers less than the square root of limit. The Wikipedia article indicates only eight removals (for limit = 150) because it doesn't count the removed primes and the initial 1 that the members array is initialized with.
using System;
using System.Collections.Generic;
class Program {
// Returns list of primes up to limit using Pritchard (wheel) sieve
static List<int> PrimesUpTo(int limit, bool verbose = false) {
var sw = System.Diagnostics.Stopwatch.StartNew();
var members = new SortedSet<int>{ 1 };
int stp = 1, prime = 2, n, nxtpr, rtlim = 1 + (int)Math.Sqrt(limit), nl, ac = 2, rc = 1;
List<int> primes = new List<int>(), tl = new List<int>();
while (prime < rtlim) {
nl = Math.Min(prime * stp, limit);
if (stp < limit) {
tl.Clear();
foreach (var w in members)
for (n = w + stp; n <= nl; n += stp) tl.Add(n);
members.UnionWith(tl); ac += tl.Count;
}
stp = nl; // update wheel size to wheel limit
nxtpr = 5; // for obtaining the next prime
tl.Clear();
foreach (var w in members) {
if (nxtpr == 5 && w > prime) nxtpr = w;
if ((n = prime * w) > nl) break; else tl.Add(n);
}
foreach (var itm in tl) members.Remove(itm); rc += tl.Count;
primes.Add(prime);
prime = prime == 2 ? 3 : nxtpr;
}
members.Remove(1); primes.AddRange(members); sw.Stop();
if (verbose) Console.WriteLine("Up to {0}, added:{1}, removed:{2}, primes counted:{3}, time:{4} ms", limit, ac, rc, primes.Count, sw.Elapsed.TotalMilliseconds);
return primes;
}
static void Main(string[] args) {
Console.WriteLine("[{0}]", string.Join(", ", PrimesUpTo(150, true)));
PrimesUpTo(1000000, true);
}
}
- Output:
Timing from Tio.run
Up to 150, added:45, removed:14, primes counted:35, time:13.2842 ms [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] Up to 1000000, added:186825, removed:108494, primes counted:78498, time:139.4323 ms
C++
The starting idea is to represent W as simply as possible, with an array w[] containing the members in order; i.e. w[i] is the i'th member (indexing from 0). When the current wheel is extended by rolling it, the code simply iterates through the array w, adding the length of the wheel to each member w[i] and appending the result. The other step is to delete the composites formed by multiplying the values in the current wheel by the current prime p. However, this presents problems, firstly because each multiple cannot be found in w in O(1) time. Accordingly, a bit array d[] (for "deleted") is introduced such that d[x] is initialized to false when a value x is appended to W, and is set to true should x be deleted as a multiple p*w[i] of p. Deletions are now fast, but the array is left containing deleted elements.
So if the new W will be extended in the next iteration, because its length < N, then the array w is compressed by eliminating the deleted values. But once the length reaches N (which happens very quickly), it would be way too costly to compress w at the end of each iteration. However, only the values in W up to N/p will be used as factors in the next lot of deletions. So it suffices to compress only this initial section of w. When the remaining primes are gathered on completion, it is necessary to skip zero and deleted values in w.
Each low-level operation in the resulting algorithm can be associated with an abstract operation so that each of the latter gets O(1) operations. So the resulting program still runs in O(N/log log N) time, and the implicit constant factor is quite small.
/* Sieve of Pritchard in C++, as described at https://en.wikipedia.org/wiki/Sieve_of_Pritchard */
/* simple but more efficient implementation than the standard one that uses an indexed linked-list */
/* 2 <= N <= 1000000000 */
/* (like the standard Sieve of Eratosthenes, this algorithm is not suitable for very large N due to memory requirements) */
#include <cstring>
#include <string>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <ctime>
#define Append(x) {\
w[++w_end] = x; /* Append x to the ordered set W */\
d[x] = false;\
}
void Extend (uint32_t w[], uint32_t &w_end, uint32_t &length, uint32_t n, uint32_t N, bool d[], uint32_t &w_end_max) {
/* Rolls full wheel W up to n, and sets length=n */
uint32_t i, x;
i = 0;
x = length + 1; /* length+w[0] */
while (x <= n) {
Append(x);
x = length + w[++i];
}
length = n;
if (length == N) Append(N+1); /* sentinel */
if (w_end > w_end_max) w_end_max = w_end;
}
void Delete (uint32_t w[], uint32_t length, uint32_t p, bool d[], uint32_t &imaxf) {
/* Deletes multiples p*w[i] of p from W, and sets imaxf to last i deleted */
uint32_t i, x;
i = 0;
x = p; /* p*w[0]=p*1 */
while (x <= length) {
d[x] = true; /* Remove x from W; */
x = p*w[++i];
}
imaxf = i-1;
}
void Compress(uint32_t w[], uint32_t length, uint32_t N, bool d[], uint32_t imax, uint32_t &w_end) {
/* Removes deleted values in w[0..imax], and if length < N, updates w_end, otherwise pads with zeros on right */
uint32_t i, j;
j = 0;
for (i=1; i <= imax; i++) {
if (!d[w[i]]) {
w[++j] = w[i];
}
}
if (length < N) {
w_end = j;
} else {
for (uint32_t k=j+1; k <= imax; k++) w[k] = 0;
}
}
void Sift(uint32_t N, bool printPrimes, uint32_t &nrPrimes, uint32_t &vBound) {
/* finds the nrPrimes primes up to N, printing them if printPrimes */
uint32_t *w = new uint32_t[N/4+5];
bool *d = new bool[N+1];
uint32_t w_end, length;
/* representation invariant (for the main loop): */
/* if length < N (so W is a complete wheel), w[0..w_end] is the ordered set W; */
/* otherwise, w[0..w_end], omitting zeros and values w with d[w] true, is the ordered set W, */
/* and no values <= N/p are omitted */
uint32_t w_end_max, p, imaxf, i;
/* W,k,length = {1},1,2: */
w_end = 0; w[0] = 1;
w_end_max = 0;
length = 2;
/* Pr = {2}: */
nrPrimes = 1;
if (printPrimes) printf("%d", 2);
p = 3;
/* invariant: p = p_(k+1) and W = W_k inter {1,...,N} and length = min(P_k,N) and Pr = the first k primes */
/* (where p_i denotes the i'th prime, W_i denotes the i'th wheel, P_i denotes the product of the first i primes) */
while (p*p <= N) {
/* Append p to Pr: */
nrPrimes++;
if (printPrimes) printf(" %d", p);
if (length < N) {
/* Extend W,length to minimum of p*length,N: */
Extend (w, w_end, length, std::min(p*length,N), N, d, w_end_max);
}
Delete(w, length, p, d, imaxf);
Compress(w, length, N, d, (length < N ? w_end : imaxf), w_end); /* (can be inefficient for last full wheel) */
/* p = next(W, 1): */
i = 1;
while (w[i] == 0) i++;
p = w[i];
/* k++ */
}
if (length < N) {
/* Extend full wheel W,length to N: */
Extend (w, w_end, length, N, N, d, w_end_max);
}
if (N == 2) { Append(N+1); w_end_max = w_end; } /* sentinel */
/* gather remaining primes: */
i = 0;
while (w[++i] <= N) {
if (w[i] == 0 || d[w[i]]) continue;
if (printPrimes) printf(" %d", w[i]);
nrPrimes++;
}
vBound = w_end_max+1;
}
int main (int argc, char *argw[]) {
uint32_t N, nrPrimes, vBound;
N = 150;
int start_s = clock();
Sift(N, true, nrPrimes, vBound);
int stop_s=clock();
printf("\n%d primes up to %lu found in %.3f ms using array w[%d]\n", nrPrimes, (unsigned long)N, (stop_s-start_s)*1000.0/double(CLOCKS_PER_SEC), vBound);
}
- Output:
Output
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 35 primes up to 150 found in 0.093 ms using array w[41]
J
Implementation model:
pritchard=: {{N=. y
root=. >.@%: N
spokes=. 1
primes=. ''
p=. 0
while. p<:root do.
primes=. primes, p=. 2+(}.spokes) i.1 NB. find next prime
rim=. #spokes NB. "length" of "circumference" of wheel
spokes=. (N<.p*rim)$spokes NB. roll next larger wheel
NB. remove multiples of this next prime:
spokes=. 0 ((#spokes) (>#]) _1+p*1+i.rim)} spokes NB. remove newly recognized prime from wheel
end.
N (>:#]) primes,1+}.,I.spokes
}}
Task example:
pritchard 150
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
However, this approach exhibits performance problems when N is large.
A faster approach recognizes when the wheel is large enough and treats all subsequent "next primes" specially:
pr =: {{N=.y
root=. <.%:N NB. performance optimzation
circumference=. 1
spokes=. ,1
primes=. ''
while. N > L=. circumference do.
primes=. primes, p =. 1{ spokes,L+1 NB. next prime from sieve
circumference=. N <. p * L NB. next larger wheel:
spokes=. circumference (>:#]), spokes +/~ L * i.circumference >.@% L
NB. remove multiples of this next prime:
spokes=. spokes -. p * spokes ( [{.~ >:@:(I.-(-.@e.)~))circumference<.@%p
end.
NB. set up for optimized version of above code
comb=. root (>:#]) }. spokes NB. candidate next primes to consider
discardp=. discard=. '' NB. what we'll be eliminating
for_p. comb do.
if. p e. comb =. comb (-. }.) discardp do.
NB. remove multiples of this next prime:
discardp=. p * spokes ( [{.~ >:@:(I.-(-.@e.)~))circumference<.@%p
discard =. discard, discardp
end.
end.
primes,comb,}.spokes-.discard
}}
Here, pr 150 gives the same result as pritchard 150 but pr 1e7 takes well under a second.
Julia
Added add/remove statistics. "Removed" figure is a combination of composites and primes under sqrt of limit. Getting a >20x speedup from using a BitArray instead of a Set.
""" Rosetta Code task rosettacode.org/wiki/Sieve_of_Pritchard """
""" Pritchard sieve of primes up to limit. Uses type of `limit` arg for type of primes """
function pritchard(limit::T, verbose=false) where {T<:Integer}
members = falses(limit)
members[1] = true
steplength = 1 # wheel size
prime = T(2)
primes = T[]
nlimit = prime * steplength # wheel limit
ac = 2 # added count, since adding 1 & 2 during initialization
rc = 1 # removed count, since 1 will be removed at the end
rtlim = T(isqrt(limit)) # this allows the main loop to go
while prime < rtlim # one extra time, eliminating the follow-up for
# the last partial wheel (if present)
if steplength < limit
for w in 1:steplength
if members[w]
n = w + steplength
while n <= nlimit
members[n] = true
ac += 1
n += steplength
end
end
end
steplength = nlimit # advance wheel size
end
np = 5
mcopy = copy(members)
for w in 1:nlimit
if mcopy[w]
np == 5 && w > prime && (np = w)
n = prime * w
n > nlimit && break
rc += 1
members[n] = false
end
end
np < prime && break
push!(primes, prime)
prime = prime == 2 ? 3 : np
nlimit = min(steplength * prime, limit) # advance wheel limit
end
members[1] = false
newprimes = [i for i in eachindex(members) if members[i]]
verbose && println(
"up to $limit, added $ac, removed $rc, prime count ",
length(primes) + length(newprimes),
)
return append!(primes, newprimes)
end
println(pritchard(150))
pritchard(1000000, true)
- Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149]up to 1000000, added 186825, removed 108494, prime count 78498
Perl
use v5.36;
use List::Util 'min';
my($limit, $maxS, $length, $p, @s) = (150, 1, 2, 3);
sub next_($w) { $s[$w-1] }
sub prev_($w) { $s[$w-2] }
sub append($w) {
$s[$maxS-1] = $w;
$s[$w-2] = $maxS;
$maxS = $w;
}
sub delete_multiples_of($p) {
my $f = $p;
while ($p*$f <= $length) { $f = next_ $f }
while ( $f > 1 ) { $f = prev_ $f; delete_pf($p*$f) }
}
sub delete_pf($pf) {
my($temp1, $temp2) = ($s[$pf-2], $s[$pf-1]);
$s[ $temp1-1 ] = $temp2;
$s[($temp2-2)%@s] = $temp1;
}
sub extend_to($n) {
my($w, $x) = (1, $length+1);
while ($x <= $n) {
append $x;
$w = next_ $w;
$x = $length + $w;
}
$length = $n;
append $limit+2 if $length == $limit
}
do {
extend_to min $p*$length, $limit if $length < $limit;
delete_multiples_of $p;
$p = next_ 1;
extend_to $limit if $length < $limit
} until $p*$p > $limit;
my @primes = 2;
for (my $p = 3; $p <= $limit; $p = next_ $p) { push @primes, $p }
say "@primes";
- Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149
Phix
with javascript_semantics function pritchard(integer limit) sequence members = {1}, primes = {} integer steplength = 1, prime = 2 while prime * prime <= limit do if steplength < limit then integer mpsll = min(prime * steplength, limit) for w in members do integer n = w + steplength while n <= mpsll do members &= n n += steplength end while end for steplength = mpsll end if members = sort(filter(members,"out",sq_mul(members,prime))) primes &= prime prime = iff(prime=2?3:members[2]) end while primes &= members[2..$] return primes end function printf(1,"%s\n",{join_by(pritchard(150),1,7," ",fmt:="%3d")})
- Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149
Python
""" Rosetta Code task rosettacode.org/wiki/Sieve_of_Pritchard """
from numpy import ndarray
from math import isqrt
def pritchard(limit):
""" Pritchard sieve of primes up to limit """
members = ndarray(limit + 1, dtype=bool)
members.fill(False)
members[1] = True
steplength, prime, rtlim, nlimit = 1, 2, isqrt(limit), 2
primes = []
while prime <= rtlim:
if steplength < limit:
for w in range(1, len(members)):
if members[w]:
n = w + steplength
while n <= nlimit:
members[n] = True
n += steplength
steplength = nlimit
np = 5
mcpy = members.copy()
for w in range(1, len(members)):
if mcpy[w]:
if np == 5 and w > prime:
np = w
n = prime * w
if n > nlimit:
break # no use trying to remove items that can't even be there
members[n] = False # no checking necessary now
if np < prime:
break
primes.append(prime)
prime = 3 if prime == 2 else np
nlimit = min(steplength * prime, limit) # advance wheel limit
newprimes = [i for i in range(2, len(members)) if members[i]]
return sorted(primes + newprimes)
print(pritchard(150))
print('Number of primes up to 1,000,000:', len(pritchard(1000000)))
- Output:
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149] Number of primes up to 1,000,000: 78498
Raku
First, a direct translation of the implementation in the YouTube video:
unit sub MAIN($limit = 150);
my $maxS = 1;
my $length = 2;
my $p = 3;
my @s = ();
while $p*$p <= $limit {
if $length < $limit {
extend-to [$p*$length, $limit].min;
}
delete-multiples-of($p);
$p = next(1);
}
if $length < $limit {
extend-to $limit;
}
# Done, build the list of actual primes from the array
$p = 3;
my @primes = 2, |gather while $p <= $limit {
take $p;
$p = next($p);
};
say @primes;
exit;
sub extend-to($n) {
my $w = 1;
my $x = $length + 1;
while $x <= $n {
append $x;
$w = next($w);
$x = $length + $w;
}
$length = $n;
if $length == $limit {
append $limit+2;
}
}
sub delete-multiples-of($p) {
my $f = $p;
while $p*$f <= $length {
$f = next($f);
}
while $f > 1 {
$f = prev($f);
delete($p*$f);
}
}
sub append($w) {
@s[$maxS-1] = $w;
@s[$w-2] = $maxS;
$maxS = $w;
}
sub next($w) { @s[$w-1]; }
sub prev($w) { @s[$w-2]; }
sub delete($pf) {
my $temp1 = @s[$pf-2];
my $temp2 = @s[$pf-1];
@s[$temp1-1] = $temp2;
@s[($temp2-2)%@s] = $temp1;
}
- Output:
[2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149
Then a slightly more Raku-ish implementation based on the description in the Wikipedia article:
unit sub MAIN($limit = 150);
class Wheel {
has $.members is rw;
has $.length is rw;
method extend(*@limits) {
my @members = $.members.keys;
for @members -> $w {
my $n = $w + $.length;
while $n <= @limits.all {
$.members.set($n);
$n += $.length;
}
}
$.length = @limits.min;
}
}
# start with W₀=({1},1)
my $wheel = Wheel.new: :members(SetHash(1)), :length(1);
my $prime = 2;
my @primes = ();
while $prime * $prime <= $limit {
if $wheel.length < $limit {
$wheel.extend($prime*$wheel.length, $limit);
}
for $wheel.members.keys.sort(-*) -> $w {
$wheel.members.unset($prime * $w);
}
@primes.push: $prime;
$prime = $prime == 2 ?? 3 !! $wheel.members.keys.grep(*>1).sort[0];
}
if $wheel.length < $limit {
$wheel.extend($limit);
}
@primes.append: $wheel.members.keys.grep: * != 1;
say @primes.sort;
The only difference in the output is that the result of `.sort` is a list rather than an array, so it's printed in parentheses instead of square brackets:
- Output:
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149)
Wren
import "./sort" for SortedList
import "./fmt" for Fmt
var extend = Fn.new { |W, length, n|
var w = 1
var x = length + 1
while (x <= n) {
W.add(x)
var ix = W.indexOf(w)
w = W[ix+1]
x = length + w
}
}
var deleteMultiples = Fn.new { |W, length, p|
var w = p
while (p * w <= length) {
var ix = W.indexOf(w)
w = W[ix+1]
}
while (w > 1) {
var ix = W.indexOf(w)
w = W[ix-1]
W.remove(p*w)
}
}
var sieveOfPritchard = Fn.new { |N|
if (N < 2) return []
var W = SortedList.fromOne(1)
var Pr = SortedList.fromOne(2)
var k = 1
var length = 2
var p = 3
while (p * p <= N) {
if (length < N) {
var n = N.min(p*length)
extend.call(W, length, n)
length = n
}
deleteMultiples.call(W, length, p)
Pr.add(p)
k = k + 1
p = W[1]
}
if (length < N) extend.call(W, length, N)
return (Pr + W)[1..-1]
}
Fmt.tprint("$3d", sieveOfPritchard.call(150), 7)- Output:
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149