# Sieve of Pritchard

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 an algorithm for finding the prime numbers up to a given limit N, published in 1981. It removes many fewer composite numbers than the Sieve of Eratosthenes (and has a better asymptotic time complexity). However, unlike the latter, it cannot be modified to greatly reduce its space requirement, making it unsuitable for very large limits.

Conceptually, it works by building a wheel representing the repeating pattern of numbers not divisible by one of the first k primes, increasing k until the square of the k'th prime is at least N. Since wheels grow very quickly, the algorithm restricts attention to the initial portions of wheels up to N. (Small examples of the wheels constructed by the Sieve of Pritchard are used in the "wheel-based optimizations" mentioned in the Eratosthenes task.)

For example, the second-order wheel has circumference 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. Those that are multiples of 5 (only 2: 1*5 and 5*5) are obtained by multiplying the members of the second-order wheel. Removing them gives the next wheel, and so on.

this YouTube video 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 wheel is represented by a sparse global array s such that for each member w of the wheel, s[w] contains the next member of the wheel; along with a similar "previous member" value, this allows numbers to be removed in a constant number of operations. But the simple abstract algorithm is based on an ordered set, and there is plenty of scope for different implementations.

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.

## AppleScript

```on sieveOfPritchard(limit)
if (limit < 2) then return {}
script o
property primes : {}
property wheel : {1, 2}
end script
set {oldCircumference, circumference} to {missing value, count o's wheel}

set prime to 1
repeat until (prime * prime > limit)
-- Get the lowest prime currently in the wheel.
set prime to o's wheel's second item
set end of o's primes to prime
-- Expand the wheel's circumference to prime times the current value. Populate with the existing
-- numbers added to multiples of the current circumference, omitting multiples of the prime.
set oldCircumference to circumference
set searchLimit to (count o's wheel)
set circumference to oldCircumference * prime
if (circumference > limit) then set circumference to limit
repeat with augend from oldCircumference to (circumference - 1) by oldCircumference
repeat with i from 1 to searchLimit
set n to augend + (o's wheel's item i)
if (n > circumference) then exit repeat
if (n mod prime > 0) then set end of o's wheel to n
end repeat
end repeat
-- Find and delete multiples of the prime occurring in the old part of the wheel.
set maxCompFactor to oldCircumference div prime
set i to 2
repeat while ((o's wheel's item i) < maxCompFactor)
set i to i + 1
end repeat
repeat with i from i to 1 by -1
set j to binarySearch((o's wheel's item i) * prime, o's wheel, i, searchLimit)
if (j > 0) then
set o's wheel's item j to missing value
set searchLimit to j - 1
end if
end repeat
set o's wheel to o's wheel's numbers
end repeat

return o's primes & rest of o's wheel
end sieveOfPritchard

on binarySearch(n, theList, l, r)
script o
property lst : theList
end script
repeat until (l = r)
set m to (l + r) div 2
if (item m of o's lst < n) then
set l to m + 1
else
set r to m
end if
end repeat

if (item l of o's lst is n) then return l
return 0
end binarySearch

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}
```

## 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;
prime = prime == 2 ? 3 : nxtpr;
}
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++

We obtain a simple but high-performance implementation.

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 high-performance implementation using a simple array of integers and a bit array (using bytes for speed). */
/* 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>

void Extend (uint32_t w[], uint32_t &w_end, uint32_t &length, uint32_t n, bool d[], uint32_t &w_end_max) {
/* Rolls full wheel W up to n, and sets length=n */
uint32_t i, j, x;
i = 0; j = w_end;
x = length + 1; /* length+w[0] */
while (x <= n) {
w[++j] = x; /* Append x to the ordered set W */
d[x] = false;
x = length + w[++i];
}
length = n; w_end = j;
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 for deletion */
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[], bool d[], uint32_t to, uint32_t &w_end) {
/* Removes deleted values in w[0..to], and if to=w_end, updates w_end, otherwise pads with zeros on right */
uint32_t i, j;
j = 0;
for (i=1; i <= to; i++) {
if (!d[w[i]]) {
w[++j] = w[i];
}
}
if (to == w_end) {
w_end = j;
} else {
for (uint32_t k=j+1; k <= to; 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;
/* 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 with length to minimum of p*length and N: */
Extend (w, w_end, length, std::min(p*length,N), d, w_end_max);
}
Delete(w, length, p, d, imaxf);
Compress(w, d, (length < N ? w_end : imaxf), w_end);
/* p = next(W, 1): */
p = w[1];
if (p == 0) break; /* next p is after zeroed section so is too big */
/* k++ */
}
if (length < N) {
/* Extend full wheel W,length to N: */
Extend (w, w_end, length, N, d, w_end_max);
}
/* gather remaining primes: */
for (uint32_t i=1; i <= w_end; i++) {
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[]) {
bool error = false;
bool printPrimes = false;
uint32_t N, nrPrimes, vBound;
if (argc == 3) {
if (strcmp(argw[2], "-p") == 0) {
printPrimes = true;
argc--;
} else {
error = true;
}
}
if (argc == 2) {
N = atoi(argw[1]);
if (N < 2 || N > 1000000000) error = true;
} else {
error = true;
}
if (error) {
printf("call with: %s N -p where 2 <= N <= 1000000000 and -p to print the primes is optional \n", argw[0]);
exit(1);
}
int start_s = clock();
Sift(N, printPrimes, 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)*1E3/double(CLOCKS_PER_SEC), vBound);
}
```
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.038 ms using array w[40]

50847534 primes up to 1000000000 found in 2339.566 ms using array w[163588196]
```

## 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
}}
```
```   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
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
end.
end.
}}
```

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

Translation of: Raku
```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

Translation of: Julia
```""" 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;
}
}

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

Library: Wren-sort
Library: Wren-fmt
```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) {
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)
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
```