Sieve of Pritchard: Difference between revisions
(Added Easylang) |
|||
(27 intermediate revisions by 13 users not shown) | |||
Line 2: | Line 2: | ||
[[Category:Simple]] |
[[Category:Simple]] |
||
The [[wp:Sieve_of_Pritchard|Sieve of Pritchard]] is |
The [[wp:Sieve_of_Pritchard|Sieve of Pritchard]] is an algorithm for finding the prime numbers up to a given limit N, published in 1981. It considers many fewer composite numbers than the [[Sieve of Eratosthenes|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 |
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 |
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. |
||
[https://www.youtube.com/watch?v= |
[https://www.youtube.com/watch?v=GxgGMwLfTjE 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 <tt>s</tt> such that for each member w of the wheel, <tt>s[w]</tt> 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. |
||
;Task: |
;Task: |
||
Line 31: | Line 31: | ||
if (limit < 2) then return {} |
if (limit < 2) then return {} |
||
script o |
script o |
||
property |
property wheel : {2} |
||
property |
property extension : missing value |
||
property oldWheel : missing value |
|||
end script |
end script |
||
set {oldCircumference, circumference} to {missing value, 2} |
|||
set {x, newCircumference, currentPrime, mv} to {0, 2, 1, missing value} |
|||
repeat until (oldCircumference = limit) |
|||
repeat until (currentPrime * currentPrime > limit) |
|||
set o's oldWheel to o's wheel's numbers |
|||
-- Get the next confirmed prime from the wheel. |
|||
set |
set x to x + 1 |
||
set |
set currentPrime to o's wheel's item x |
||
-- Get an extension list nominally expanding the wheel to the lesser of |
|||
set circumference to oldCircumference * prime |
|||
-- its current circumference * currentPrime and the limit parameter. |
|||
-- It'll be far longer than needed, but hey. |
|||
set oldCircumference to newCircumference |
|||
set newCircumference to oldCircumference * currentPrime |
|||
if (newCircumference > limit) then set newCircumference to limit |
|||
else |
|||
set o's extension to makeList(newCircumference - oldCircumference, mv) |
|||
-- Insert numbers that are oldCircumference added to 1 and to each number currently in the |
|||
-- unpartitioned part of the wheel, except where the results are multiples of currentPrime. |
|||
set k to 0 |
|||
set listLen to (count o's wheel) |
|||
repeat with augend from oldCircumference to (newCircumference - 1) by oldCircumference |
|||
set n to augend + 1 |
|||
if (n mod currentPrime > 0) then |
|||
set k to k + 1 |
|||
set o's extension's item k to n |
|||
end if |
end if |
||
repeat with i from x to listLen |
|||
set n to augend + (o's wheel's item i) |
|||
if (n > newCircumference) then exit repeat |
|||
if (n mod currentPrime > 0) then |
|||
set k to k + 1 |
|||
set o's extension's item k to n |
|||
end if |
|||
end repeat |
|||
end repeat |
end repeat |
||
-- Find and delete multiples of the current prime which occur in the old part of the wheel. |
|||
repeat with this in o's oldWheel |
|||
set maxMultiple to oldCircumference div currentPrime |
|||
set i to x |
|||
if (n > circumference) then exit repeat |
|||
repeat while ((o's wheel's item i) < maxMultiple) |
|||
set i to i + 1 |
|||
end repeat |
end repeat |
||
repeat with i from i to x by -1 |
|||
set j to binarySearch((o's wheel's item i) * currentPrime, o's wheel, i, listLen) |
|||
if (j > 0) then |
|||
set o's wheel's item j to mv |
|||
end sieveOfPritchard |
|||
set listLen to j - 1 |
|||
sieveOfPritchard(150)</syntaxhighlight> |
|||
{{output}} |
|||
<syntaxhighlight lang="applescript">{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}</syntaxhighlight> |
|||
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. |
|||
<syntaxhighlight lang="applescript">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 if |
||
end repeat |
end repeat |
||
-- Keep the undeleted numbers and any in the extension list. |
|||
repeat with this in o's oldWheel |
|||
set o's wheel to o's wheel's numbers |
|||
if (k > 0) then set o's wheel to o's wheel & o's extension's items 1 thru k |
|||
set o's wheel's item n to missing value |
|||
end repeat |
|||
end repeat |
end repeat |
||
return |
return o's wheel |
||
end sieveOfPritchard |
end sieveOfPritchard |
||
Line 115: | Line 102: | ||
end makeList |
end makeList |
||
on binarySearch(n, theList, l, r) |
|||
sieveOfPritchard(1000000)</syntaxhighlight> |
|||
script o |
|||
property lst : theList |
|||
end script |
|||
repeat until (l = r) |
|||
set m to (l + r) div 2 |
|||
if (o's lst's item m < n) then |
|||
set l to m + 1 |
|||
else |
|||
set r to m |
|||
end if |
|||
end repeat |
|||
if (o's lst's item l is n) then return l |
|||
return 0 |
|||
end binarySearch |
|||
sieveOfPritchard(150)</syntaxhighlight> |
|||
{{output}} |
|||
<syntaxhighlight lang="applescript">{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}</syntaxhighlight> |
|||
=={{header|BASIC}}== |
|||
==={{header|BASIC256}}=== |
|||
<syntaxhighlight lang="vb">arraybase 1 |
|||
call sieveOfPritchard(150, true) |
|||
call sieveOfPritchard(1e6, false) |
|||
end |
|||
function min(a, b) |
|||
if a < b then return a else return b |
|||
end function |
|||
subroutine sieveOfPritchard(limit, imprime) |
|||
dim members[limit + 1] fill false |
|||
members[1] = true |
|||
ub = members[?] |
|||
stepLength = 1 |
|||
prime = 2 |
|||
rtlim = sqr(limit) |
|||
nlimit = 2 |
|||
dim primes[1] |
|||
cont = 0 |
|||
while prime <= rtlim |
|||
if stepLength < limit then |
|||
for w = 1 to ub |
|||
if members[w] then |
|||
dim n = w + stepLength |
|||
while n <= nlimit |
|||
members[n] = true |
|||
n += stepLength |
|||
end while |
|||
end if |
|||
next |
|||
stepLength = nlimit |
|||
end if |
|||
np = 5 |
|||
dim mcpy[ub] |
|||
for i = 1 to ub |
|||
mcpy[i] = members[i] |
|||
next |
|||
for i = 1 to ub |
|||
if mcpy[i] then |
|||
if np = 5 and i > prime then np = i |
|||
dim n = prime * i |
|||
if n > limit then exit for |
|||
members[n] = false |
|||
end if |
|||
next |
|||
if np < prime then exit while |
|||
cont += 1 |
|||
redim primes(cont) |
|||
primes[cont] = prime |
|||
if prime = 2 then prime = 3 else prime = np |
|||
nlimit = min(stepLength * prime, limit) |
|||
end while |
|||
dim newPrimes(ub) |
|||
for i = 2 to ub |
|||
if members[i] then newPrimes[i] = i |
|||
next |
|||
cont = 0 |
|||
for i = 1 to primes[?] |
|||
if imprime then print " "; primes[i]; |
|||
cont += 1 |
|||
next |
|||
for i = 1 to ub |
|||
if newPrimes[i] then |
|||
cont += 1 |
|||
if imprime then print " "; i; |
|||
end if |
|||
next |
|||
if not imprime then print : print "Number of primes up to "; limit; ": "; cont |
|||
end subroutine</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Similar to FreeBASIC entry.</pre> |
|||
==={{header|FreeBASIC}}=== |
|||
{{trans|Wren}} |
|||
<syntaxhighlight lang="vb">#define min(a, b) iif((a) < (b), (a), (b)) |
|||
Sub sieveOfPritchard(limit As Uinteger, imprime As Boolean) |
|||
Dim As Boolean members(1 To limit + 1) |
|||
members(1) = True |
|||
Dim As Uinteger ub = Ubound(members) |
|||
Dim As Uinteger stepLength = 1 |
|||
Dim As Uinteger prime = 2 |
|||
Dim As Uinteger rtlim = Sqr(limit) |
|||
Dim As Uinteger nlimit = 2 |
|||
Dim As Integer primes() |
|||
Dim As Integer i, cont = 0 |
|||
While prime <= rtlim |
|||
If stepLength < limit Then |
|||
For w As Integer = 1 To ub |
|||
If members(w) Then |
|||
Dim As Integer n = w + stepLength |
|||
While n <= nlimit |
|||
members(n) = True |
|||
n += stepLength |
|||
Wend |
|||
End If |
|||
Next |
|||
stepLength = nlimit |
|||
End If |
|||
Dim As Uinteger np = 5 |
|||
Dim As Boolean mcpy(ub) |
|||
For i = 1 To ub |
|||
mcpy(i) = members(i) |
|||
Next |
|||
For i = 1 To ub |
|||
If mcpy(i) Then |
|||
If np = 5 And i > prime Then np = i |
|||
Dim As Uinteger n = prime * i |
|||
If n > limit Then Exit For there. |
|||
members(n) = False |
|||
End If |
|||
Next |
|||
If np < prime Then Exit While |
|||
cont += 1 |
|||
Redim Preserve primes(cont) |
|||
primes(cont) = prime |
|||
prime = Iif(prime = 2, 3, np) |
|||
nlimit = min(stepLength * prime, limit) |
|||
Wend |
|||
Dim As Integer newPrimes(ub) |
|||
For i = 2 To ub |
|||
If members(i) Then newPrimes(i) = i |
|||
Next |
|||
cont = 0 |
|||
For i = 1 To Ubound(primes) |
|||
If imprime Then Print primes(i); |
|||
cont += 1 |
|||
Next |
|||
For i = 1 To ub |
|||
If newPrimes(i) Then |
|||
cont += 1 |
|||
If imprime Then Print i; |
|||
End If |
|||
Next |
|||
If Not imprime Then Print !"\nNumber of primes up to "; limit; ":"; cont |
|||
End Sub |
|||
sieveOfPritchard(150, True) |
|||
sieveOfPritchard(1e6, False) |
|||
Sleep</syntaxhighlight> |
|||
{{out}} |
|||
<pre> 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 |
|||
Number of primes up to 1000000: 78498</pre> |
|||
==={{header|True BASIC}}=== |
|||
<syntaxhighlight lang="qbasic">FUNCTION min(a, b) |
|||
IF a < b THEN LET min = a ELSE LET min = b |
|||
END FUNCTION |
|||
SUB sieveofpritchard (limit) |
|||
DIM members(0) |
|||
MAT REDIM members(limit+1) |
|||
LET members(1) = 1 |
|||
LET ub = UBOUND(members) |
|||
LET steplength = 1 |
|||
LET prime = 2 |
|||
LET rtlim = SQR(limit) |
|||
LET nlimit = 2 |
|||
DIM primes(10) |
|||
LET cnt = 0 |
|||
DIM mcpy(0) |
|||
MAT REDIM mcpy(ub) |
|||
DO WHILE prime <= rtlim |
|||
IF steplength < limit THEN |
|||
FOR w = 1 TO ub |
|||
IF members(w)<>0 THEN |
|||
LET n = w+steplength |
|||
DO WHILE n <= nlimit |
|||
LET members(n) = 1 |
|||
LET n = n+steplength |
|||
LOOP |
|||
END IF |
|||
NEXT w |
|||
LET steplength = nlimit |
|||
END IF |
|||
LET np = 5 |
|||
FOR i = 1 TO ub |
|||
LET mcpy(i) = members(i) |
|||
NEXT i |
|||
FOR i = 1 TO ub |
|||
IF mcpy(i)<>0 THEN |
|||
IF np = 5 AND i > prime THEN LET np = i |
|||
LET n = prime*i |
|||
IF n > limit THEN EXIT FOR |
|||
LET members(n) = 0 |
|||
END IF |
|||
NEXT i |
|||
IF np < prime THEN EXIT DO |
|||
LET cnt = cnt+1 |
|||
MAT REDIM primes(cnt) |
|||
LET primes(cnt) = prime |
|||
IF prime = 2 THEN LET prime = 3 ELSE LET prime = np |
|||
LET nlimit = min(steplength*prime, limit) |
|||
LOOP |
|||
DIM newprimes(0) |
|||
MAT REDIM newprimes(ub) |
|||
FOR i = 2 TO ub |
|||
IF members(i)<>0 THEN LET newprimes(i) = i |
|||
NEXT i |
|||
LET cnt = 0 |
|||
FOR i = 1 TO UBOUND(primes) |
|||
PRINT primes(i); |
|||
LET cnt = cnt+1 |
|||
NEXT i |
|||
FOR i = 1 TO ub |
|||
IF newprimes(i)<>0 THEN |
|||
PRINT i; |
|||
LET cnt = cnt+1 |
|||
END IF |
|||
NEXT i |
|||
END SUB |
|||
CLEAR |
|||
CALL sieveofpritchard (150) |
|||
END</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Similar to FreeBASIC entry.</pre> |
|||
==={{header|Yabasic}}=== |
|||
<syntaxhighlight lang="vb">sieveOfPritchard(150, true) |
|||
sieveOfPritchard(1e6, false) |
|||
end |
|||
sub sieveOfPritchard(limit, imprime) |
|||
dim members(limit + 1) |
|||
members(1) = true |
|||
ub = arraysize(members(),1) |
|||
stepLength = 1 |
|||
prime = 2 |
|||
rtlim = sqrt(limit) |
|||
nlimit = 2 |
|||
dim primes(1) |
|||
cont = 0 |
|||
while prime <= rtlim |
|||
if stepLength < limit then |
|||
for w = 1 to ub |
|||
if members(w) then |
|||
n = w + stepLength |
|||
while n <= nlimit |
|||
members(n) = true |
|||
n = n + stepLength |
|||
wend |
|||
fi |
|||
next |
|||
stepLength = nlimit |
|||
fi |
|||
np = 5 |
|||
dim mcpy(ub) |
|||
for i = 1 to ub |
|||
mcpy(i) = members(i) |
|||
next |
|||
for i = 1 to ub |
|||
if mcpy(i) then |
|||
if np = 5 and i > prime np = i |
|||
n = prime * i |
|||
if n > limit break |
|||
members(n) = false |
|||
fi |
|||
next |
|||
if np < prime break |
|||
cont = cont + 1 |
|||
redim primes(cont) |
|||
primes(cont) = prime |
|||
if prime = 2 then prime = 3 else prime = np : fi |
|||
nlimit = min(stepLength * prime, limit) |
|||
wend |
|||
dim newPrimes(ub) |
|||
for i = 2 to ub |
|||
if members(i) newPrimes(i) = i |
|||
next |
|||
cont = 0 |
|||
for i = 1 to arraysize(primes(),1) |
|||
if imprime print " ", primes(i); |
|||
cont = cont + 1 |
|||
next |
|||
for i = 1 to ub |
|||
if newPrimes(i) then |
|||
cont = cont + 1 |
|||
if imprime print " ", i; |
|||
fi |
|||
next |
|||
if not imprime then print : print "Number of primes up to ", limit, ": ", cont : fi |
|||
end sub</syntaxhighlight> |
|||
{{out}} |
|||
<pre>Similar to FreeBASIC entry.</pre> |
|||
=={{header|C#|CSharp}}== |
=={{header|C#|CSharp}}== |
||
Line 121: | Line 442: | ||
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. |
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. |
|||
<syntaxhighlight lang="csharp">using System; |
<syntaxhighlight lang="csharp">using System; |
||
using System.Collections.Generic; |
using System.Collections.Generic; |
||
Line 127: | Line 450: | ||
// Returns list of primes up to limit using Pritchard (wheel) sieve |
// Returns list of primes up to limit using Pritchard (wheel) sieve |
||
static List<int> PrimesUpTo(int limit) { |
static List<int> PrimesUpTo(int limit, bool verbose = false) { |
||
var sw = System.Diagnostics.Stopwatch.StartNew(); |
|||
var members = new SortedSet<int>{ 1 }; |
var members = new SortedSet<int>{ 1 }; |
||
int stp = 1, prime = 2, n, nxtpr, rtlim = 1 + (int)Math.Sqrt(limit), nl = 2; |
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 < |
while (prime < rtlim) { |
||
nl = Math.Min(prime * stp, limit); |
|||
if (stp < limit) { |
if (stp < limit) { |
||
tl.Clear(); |
|||
foreach (var w in members) |
foreach (var w in members) |
||
for (n = w + stp; n <= nl; n += stp) |
for (n = w + stp; n <= nl; n += stp) tl.Add(n); |
||
members.UnionWith( |
members.UnionWith(tl); ac += tl.Count; |
||
} |
} |
||
stp = nl; // update wheel size to wheel limit |
stp = nl; // update wheel size to wheel limit |
||
nxtpr = |
nxtpr = 5; // for obtaining the next prime |
||
tl.Clear(); |
|||
foreach (var w in members) { |
foreach (var w in members) { |
||
if (nxtpr == |
if (nxtpr == 5 && w > prime) nxtpr = w; |
||
if ( |
if ((n = prime * w) > nl) break; else tl.Add(n); |
||
} |
} |
||
foreach (var itm in |
foreach (var itm in tl) members.Remove(itm); rc += tl.Count; |
||
primes.Add(prime); |
primes.Add(prime); |
||
prime = prime == 2 ? 3 : nxtpr; |
prime = prime == 2 ? 3 : nxtpr; |
||
nl = Math.Min(limit, prime * stp); |
|||
} |
} |
||
members.Remove(1); |
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); |
|||
primes.AddRange(members); |
|||
return primes; |
return primes; |
||
} |
} |
||
static void Main(string[] args) { |
static void Main(string[] args) { |
||
Console.WriteLine("[{0}]", string.Join(", ", PrimesUpTo(150))); |
Console.WriteLine("[{0}]", string.Join(", ", PrimesUpTo(150, true))); |
||
PrimesUpTo(1000000, true); |
|||
} |
} |
||
}</syntaxhighlight> |
|||
{{out}}Timing from Tio.run: |
|||
<pre>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 |
|||
</pre> |
|||
=={{header|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. |
|||
<syntaxhighlight lang="cpp">/* 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); |
|||
}</syntaxhighlight> |
}</syntaxhighlight> |
||
{{out}} |
{{out}} |
||
<pre> |
<pre>2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 |
||
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]</pre> |
|||
=={{header|EasyLang}}== |
|||
{{trans|Julia}} |
|||
<syntaxhighlight> |
|||
proc pritchard limit . primes[] . |
|||
len members[] limit |
|||
members[1] = 1 |
|||
steplength = 1 |
|||
prime = 2 |
|||
rtlimit = sqrt limit |
|||
nlimit = 2 |
|||
primes[] = [ ] |
|||
while prime <= rtlimit |
|||
if steplength < limit |
|||
for w to len members[] |
|||
if members[w] = 1 |
|||
n = w + steplength |
|||
while n <= nlimit |
|||
members[n] = 1 |
|||
n += steplength |
|||
. |
|||
. |
|||
. |
|||
steplength = nlimit |
|||
. |
|||
np = 5 |
|||
mcpy[] = members[] |
|||
for w to nlimit |
|||
if mcpy[w] = 1 |
|||
if np = 5 and w > prime |
|||
np = w |
|||
. |
|||
n = prime * w |
|||
if n > nlimit |
|||
break 1 |
|||
. |
|||
members[n] = 0 |
|||
. |
|||
. |
|||
if np < prime |
|||
break 1 |
|||
. |
|||
primes[] &= prime |
|||
if prime = 2 |
|||
prime = 3 |
|||
else |
|||
prime = np |
|||
. |
|||
nlimit = lower (steplength * prime) limit |
|||
. |
|||
for i = 2 to len members[] |
|||
if members[i] = 1 |
|||
primes[] &= i |
|||
. |
|||
. |
|||
. |
|||
pritchard 150 p[] |
|||
print p[] |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
[ 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149 ] |
|||
</pre> |
|||
=={{header|J}}== |
=={{header|J}}== |
||
Implementation |
Implementation model: |
||
<syntaxhighlight lang="j">pritchard=: {{N=. y |
|||
spokes=. $.6$4{.1 |
|||
root=. >.@%: N |
|||
spokes=. 1 |
|||
while. y > #spokes do. |
|||
primes=. '' |
|||
p=. 0 |
|||
while. p<:root do. |
|||
primes=. primes, p=. 2+(}.spokes) i.1 NB. find next prime |
primes=. primes, p=. 2+(}.spokes) i.1 NB. find next prime |
||
rim=. #spokes NB. "length" of "circumference" of wheel |
rim=. #spokes NB. "length" of "circumference" of wheel |
||
spokes=. ( |
spokes=. (N<.p*rim)$spokes NB. roll next larger wheel |
||
NB. remove multiples of this next prime: |
|||
spokes=. 8 $.0 ((#~ y>])_1+p*1+i.rim)} spokes NB. remove newly recognized prime from wheel |
|||
spokes=. 0 ((#spokes) (>#]) _1+p*1+i.rim)} spokes NB. remove newly recognized prime from wheel |
|||
end. |
end. |
||
N (>:#]) primes,1+}.,I.spokes |
|||
while. y > p*p do. |
|||
primes=. primes, p=. 2+(}.spokes) i.1 NB. find next prime |
|||
spokes=. 0 ((#~ y>])_1+p*1+i.rim)} spokes NB. scrub it out of wheel |
|||
end. |
|||
primes,1+}.,I.spokes |
|||
}}</syntaxhighlight> |
}}</syntaxhighlight> |
||
Line 182: | Line 732: | ||
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</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 127 131 137 139 149</syntaxhighlight> |
||
However, this approach exhibits performance problems when N is large. |
|||
This task exposed a bug in J's implementation of X i. Y when X is an empty sparse array. To avoid/work-around this bug, we initialize <code>spokes</code> and <code>primes</code> with the values that they would have had after the first two iterations of this algorithm if <code>i.</code> would have worked properly. (The bug has been fixed, but the fix is not yet widely deployed. (We could also have used a dense array instead of a sparse array, but it's probably more "in spirit with the design of this task" to "hard code" the first two rounds and use sparse arrays than it would be to use dense arrays for the wheels.)) |
|||
A faster approach recognizes when the wheel is large enough and treats all subsequent "next primes" specially: |
|||
<syntaxhighlight lang="j">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 |
|||
}}</syntaxhighlight> |
|||
Here, <code>pr 150</code> gives the same result as <code>pritchard 150</code> but <code>pr 1e7</code> takes well under a second. |
|||
=={{header|Java}}== |
|||
<syntaxhighlight lang="java"> |
|||
import java.util.ArrayList; |
|||
import java.util.BitSet; |
|||
import java.util.List; |
|||
public final class SieveOfPritchard { |
|||
public static void main(String[] args) { |
|||
System.out.println(sieveOfPritchard(150) + System.lineSeparator()); |
|||
System.out.println("Number of primes up to 1,000,000 is " + sieveOfPritchard(1_000_000).size() + "."); |
|||
System.out.println(); |
|||
final long start = System.currentTimeMillis(); |
|||
System.out.print("Number of primes up to 100,000,000 is " + sieveOfPritchard(100_000_000).size()); |
|||
final long finish = System.currentTimeMillis(); |
|||
System.out.println(". Obtained in a time of " + ( (double) finish - start ) / 1_000 + " seconds."); |
|||
} |
|||
private static List<Integer> sieveOfPritchard(int limit) { |
|||
List<Integer> primes = new ArrayList<Integer>(); |
|||
BitSet members = new BitSet(limit + 1); |
|||
members.set(1); |
|||
List<Integer> deletions = new ArrayList<Integer>(); |
|||
final int rootLimit = (int) Math.sqrt(limit); |
|||
int nLimit = 2; |
|||
int stepLength = 1; |
|||
int prime = 2; |
|||
while ( prime <= rootLimit ) { |
|||
if ( stepLength < limit ) { |
|||
for ( int w = 1; w >= 0; w = members.nextSetBit(w + 1) ) { |
|||
int n = w + stepLength; |
|||
while ( n <= nLimit ) { |
|||
members.set(n); |
|||
n += stepLength; |
|||
} |
|||
} |
|||
stepLength = nLimit; |
|||
} |
|||
deletions.clear(); |
|||
int nextPrime = 5; |
|||
for ( int w = 1; w < nLimit; w = members.nextSetBit(w + 1) ) { |
|||
if ( nextPrime == 5 && w > prime ) { |
|||
nextPrime = w; |
|||
} |
|||
final int n = prime * w; |
|||
if ( n > nLimit ) { |
|||
break; |
|||
} |
|||
deletions.add(n); |
|||
} |
|||
for ( int deletion : deletions ) { |
|||
members.clear(deletion); |
|||
} |
|||
if ( nextPrime < prime ) { |
|||
break; |
|||
} |
|||
primes.add(prime); |
|||
prime = ( prime == 2 ) ? 3 : nextPrime; |
|||
nLimit = (int) Math.min((long) stepLength * prime, limit); |
|||
} |
|||
if ( stepLength < limit ) { |
|||
for ( int w = 1; w >= 0; w = members.nextSetBit(w + 1) ) { |
|||
int n = w + stepLength; |
|||
while ( n <= limit ) { |
|||
members.set(n); |
|||
n += stepLength; |
|||
} |
|||
} |
|||
} |
|||
members.clear(1); |
|||
for ( int i = members.nextSetBit(0); i >= 0; i = members.nextSetBit(i + 1) ) { |
|||
primes.add(i); |
|||
}; |
|||
return primes; |
|||
} |
|||
} |
|||
</syntaxhighlight> |
|||
{{ out }} |
|||
<pre> |
|||
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149] |
|||
Number of primes up to 1,000,000 is 78498. |
|||
Number of primes up to 100,000,000 is 5761455. Obtained in a time of 0.648 seconds. |
|||
</pre> |
|||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Line 190: | Line 865: | ||
""" Pritchard sieve of primes up to limit. Uses type of `limit` arg for type of primes """ |
""" 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} |
function pritchard(limit::T, verbose=false) where {T<:Integer} |
||
members = falses(limit |
members = falses(limit) |
||
members[1] = true |
members[1] = true |
||
steplength = 1 # wheel size |
steplength = 1 # wheel size |
||
Line 199: | Line 874: | ||
rc = 1 # removed count, since 1 will be removed at the end |
rc = 1 # removed count, since 1 will be removed at the end |
||
rtlim = T(isqrt(limit)) # this allows the main loop to go |
rtlim = T(isqrt(limit)) # this allows the main loop to go |
||
while prime < |
while prime < rtlim # one extra time, eliminating the follow-up for |
||
# the last partial wheel (if present) |
# the last partial wheel (if present) |
||
if steplength < limit |
if steplength < limit |
||
for w in 1: |
for w in 1:steplength |
||
if members[w] |
if members[w] |
||
n = w + steplength |
n = w + steplength |
||
Line 216: | Line 891: | ||
np = 5 |
np = 5 |
||
mcopy = copy(members) |
mcopy = copy(members) |
||
for w in 1: |
for w in 1:nlimit |
||
if mcopy[w] |
if mcopy[w] |
||
np == 5 && w > prime && (np = w) |
np == 5 && w > prime && (np = w) |
||
Line 236: | Line 911: | ||
length(primes) + length(newprimes), |
length(primes) + length(newprimes), |
||
) |
) |
||
return |
return append!(primes, newprimes) |
||
end |
end |
||
Line 244: | Line 919: | ||
up to 1000000, added 186825, removed 108494, prime count 78498 |
up to 1000000, added 186825, removed 108494, prime count 78498 |
||
</pre> |
</pre> |
||
=={{header|Nim}}== |
|||
{{trans|Python}} |
|||
<syntaxhighlight lang="Nim">import std/[algorithm, math, sugar] |
|||
proc pritchard(limit: Natural): seq[int] = |
|||
## Pritchard sieve of primes up to "limit". |
|||
var members = newSeq[bool](limit + 1) |
|||
members[1] = true |
|||
var |
|||
stepLength = 1 |
|||
prime = 2 |
|||
rtlim = sqrt(limit.toFloat).int |
|||
nlimit = 2 |
|||
primes: seq[int] |
|||
while prime <= rtlim: |
|||
if stepLength < limit: |
|||
for w in 1..members.high: |
|||
if members[w]: |
|||
var n = w + stepLength |
|||
while n <= nlimit: |
|||
members[n] = true |
|||
inc n, stepLength |
|||
stepLength = nlimit |
|||
var np = 5 |
|||
var mcpy = members |
|||
for w in 1..members.high: |
|||
if mcpy[w]: |
|||
if np == 5 and w > prime: |
|||
np = w |
|||
let n = prime * w |
|||
if n > limit: |
|||
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.add prime |
|||
prime = if prime == 2: 3 else: np |
|||
nlimit = min(stepLength * prime, limit) # Advance wheel limit. |
|||
let newPrimes = collect: |
|||
for i in 2..members.high: |
|||
if members[i]: i |
|||
result = sorted(primes & newPrimes) |
|||
echo pritchard(150) |
|||
echo "Number of primes up to 1_000_000: ", pritchard(1_000_000).len |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre>@[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149] |
|||
Number of primes up to 1_000_000: 78498 |
|||
</pre> |
|||
=={{header|Pascal}}== |
|||
{{works with|Free Pascal}} |
|||
A console program written in Delphi 7. It follows the algorithm in the Wikipedia article "Sieve of Pritchard". |
|||
<syntaxhighlight lang="pascal"> |
|||
program Pritchard_console; |
|||
{$APPTYPE CONSOLE} |
|||
uses |
|||
Math, SysUtils, Types; |
|||
// Function to return an array of all primes <= N |
|||
function PritchardSieve( const N : integer) : Types.TIntegerDynArray; |
|||
var |
|||
j, j_max, k, len, nrPrimes, p : integer; |
|||
marked : Types.TBooleanDynArray; |
|||
smallPrimes : Types.TIntegerDynArray; // i.e. primes <= sqrt( N) |
|||
spi : integer; // index into array smallPrimes |
|||
const |
|||
SP_STEP = 16; // step when extending dynamic array smallPrimes |
|||
begin |
|||
// Deal with trivial input |
|||
result := nil; |
|||
if (N <= 1) then exit; |
|||
// Initialize |
|||
SetLength( marked, N + 1); // 0..N for convenience; marked[0] is not used |
|||
marked[1] := true; // no other initialization of "marked" is needed |
|||
len := 1; |
|||
p := 2; |
|||
SetLength( smallPrimes, SP_STEP); |
|||
spi := 0; |
|||
while p*p <= N do begin |
|||
// Roll the wheel |
|||
if len < N then begin |
|||
j_max := Math.Min( p*len, N); |
|||
for j := len + 1 to j_max do marked[j] := marked[j - len]; |
|||
len := j_max; |
|||
end; |
|||
// Unmark multiples of p |
|||
for k := len div p downto 1 do |
|||
if marked[k] then marked[p*k] := false; |
|||
// Store the prime p, extending the array if necessary |
|||
if spi = Length( smallPrimes) then |
|||
SetLength( smallPrimes, spi + SP_STEP); |
|||
smallPrimes[spi] := p; |
|||
inc(spi); |
|||
// Find the next prime p |
|||
if p = 2 then p := 3 |
|||
else repeat inc(p) until (p > N) or marked[p]; |
|||
// Condition p > N is a safety net; should always hit a marked value |
|||
Assert(p <= N); |
|||
end; // while |
|||
// Final roll, if needed. It is not needed if N >= 49. This is because |
|||
// 2 < 3^2, 2*3 < 5^2, 2*3*5 < 7^2, but thereafter 2*3*5*7 > 11^2, etc. |
|||
if len < N then |
|||
for j := len + 1 to N do marked[j] := marked[j - len]; |
|||
// Remove 1 and put the small primes back |
|||
marked[1] := false; |
|||
for k := 0 to spi - 1 do marked[smallPrimes[k]] := true; |
|||
// Use the boolean array to return an array of prime integers |
|||
nrPrimes := 0; |
|||
for j := 2 to N do |
|||
if marked[j] then inc( nrPrimes); |
|||
SetLength( result, nrPrimes); |
|||
k := 0; |
|||
for j := 2 to N do |
|||
if marked[j] then begin result[k] := j; inc(k); end; |
|||
end; |
|||
// Main routine. User types the program name, |
|||
// optionally followed by the limit N (defaults to 150) |
|||
var |
|||
N, j : integer; |
|||
primes : Types.TIntegerDynArray; |
|||
begin |
|||
if ParamCount = 0 then N := 150 |
|||
else N := SysUtils.StrToInt( ParamStr(1)); |
|||
primes := PritchardSieve(N); |
|||
WriteLn( 'Number of primes = ', Length(primes)); |
|||
for j := 0 to Length(primes) - 1 do begin |
|||
Write( ' ', primes[j]:4); |
|||
if j mod 10 = 9 then WriteLn; |
|||
end; |
|||
end. |
|||
</syntaxhighlight> |
|||
{{out}} |
|||
<pre> |
|||
Number of primes = 35 |
|||
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 |
|||
</pre> |
|||
==={{header|Free Pascal}}=== |
|||
see [[Sieve_of_Eratosthenes#alternative_using_wheel]] |
|||
=={{header|Perl}}== |
|||
{{trans|Raku}} |
|||
<syntaxhighlight lang="perl" line>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";</syntaxhighlight> |
|||
{{out}} |
|||
<pre>2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127 131 137 139 149</pre> |
|||
=={{header|Phix}}== |
=={{header|Phix}}== |
||
Line 459: | Line 1,348: | ||
{{libheader|Wren-sort}} |
{{libheader|Wren-sort}} |
||
{{libheader|Wren-fmt}} |
{{libheader|Wren-fmt}} |
||
<syntaxhighlight lang=" |
<syntaxhighlight lang="wren">import "./sort" for SortedList |
||
import "./fmt" for Fmt |
import "./fmt" for Fmt |
||
Latest revision as of 14:56, 25 February 2024
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 considers 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.
- 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 wheel : {2}
property extension : missing value
end script
set {x, newCircumference, currentPrime, mv} to {0, 2, 1, missing value}
repeat until (currentPrime * currentPrime > limit)
-- Get the next confirmed prime from the wheel.
set x to x + 1
set currentPrime to o's wheel's item x
-- Get an extension list nominally expanding the wheel to the lesser of
-- its current circumference * currentPrime and the limit parameter.
-- It'll be far longer than needed, but hey.
set oldCircumference to newCircumference
set newCircumference to oldCircumference * currentPrime
if (newCircumference > limit) then set newCircumference to limit
set o's extension to makeList(newCircumference - oldCircumference, mv)
-- Insert numbers that are oldCircumference added to 1 and to each number currently in the
-- unpartitioned part of the wheel, except where the results are multiples of currentPrime.
set k to 0
set listLen to (count o's wheel)
repeat with augend from oldCircumference to (newCircumference - 1) by oldCircumference
set n to augend + 1
if (n mod currentPrime > 0) then
set k to k + 1
set o's extension's item k to n
end if
repeat with i from x to listLen
set n to augend + (o's wheel's item i)
if (n > newCircumference) then exit repeat
if (n mod currentPrime > 0) then
set k to k + 1
set o's extension's item k to n
end if
end repeat
end repeat
-- Find and delete multiples of the current prime which occur in the old part of the wheel.
set maxMultiple to oldCircumference div currentPrime
set i to x
repeat while ((o's wheel's item i) < maxMultiple)
set i to i + 1
end repeat
repeat with i from i to x by -1
set j to binarySearch((o's wheel's item i) * currentPrime, o's wheel, i, listLen)
if (j > 0) then
set o's wheel's item j to mv
set listLen to j - 1
end if
end repeat
-- Keep the undeleted numbers and any in the extension list.
set o's wheel to o's wheel's numbers
if (k > 0) then set o's wheel to o's wheel & o's extension's items 1 thru k
end repeat
return o's wheel
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
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 (o's lst's item m < n) then
set l to m + 1
else
set r to m
end if
end repeat
if (o's lst's item l 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}
BASIC
BASIC256
arraybase 1
call sieveOfPritchard(150, true)
call sieveOfPritchard(1e6, false)
end
function min(a, b)
if a < b then return a else return b
end function
subroutine sieveOfPritchard(limit, imprime)
dim members[limit + 1] fill false
members[1] = true
ub = members[?]
stepLength = 1
prime = 2
rtlim = sqr(limit)
nlimit = 2
dim primes[1]
cont = 0
while prime <= rtlim
if stepLength < limit then
for w = 1 to ub
if members[w] then
dim n = w + stepLength
while n <= nlimit
members[n] = true
n += stepLength
end while
end if
next
stepLength = nlimit
end if
np = 5
dim mcpy[ub]
for i = 1 to ub
mcpy[i] = members[i]
next
for i = 1 to ub
if mcpy[i] then
if np = 5 and i > prime then np = i
dim n = prime * i
if n > limit then exit for
members[n] = false
end if
next
if np < prime then exit while
cont += 1
redim primes(cont)
primes[cont] = prime
if prime = 2 then prime = 3 else prime = np
nlimit = min(stepLength * prime, limit)
end while
dim newPrimes(ub)
for i = 2 to ub
if members[i] then newPrimes[i] = i
next
cont = 0
for i = 1 to primes[?]
if imprime then print " "; primes[i];
cont += 1
next
for i = 1 to ub
if newPrimes[i] then
cont += 1
if imprime then print " "; i;
end if
next
if not imprime then print : print "Number of primes up to "; limit; ": "; cont
end subroutine
- Output:
Similar to FreeBASIC entry.
FreeBASIC
#define min(a, b) iif((a) < (b), (a), (b))
Sub sieveOfPritchard(limit As Uinteger, imprime As Boolean)
Dim As Boolean members(1 To limit + 1)
members(1) = True
Dim As Uinteger ub = Ubound(members)
Dim As Uinteger stepLength = 1
Dim As Uinteger prime = 2
Dim As Uinteger rtlim = Sqr(limit)
Dim As Uinteger nlimit = 2
Dim As Integer primes()
Dim As Integer i, cont = 0
While prime <= rtlim
If stepLength < limit Then
For w As Integer = 1 To ub
If members(w) Then
Dim As Integer n = w + stepLength
While n <= nlimit
members(n) = True
n += stepLength
Wend
End If
Next
stepLength = nlimit
End If
Dim As Uinteger np = 5
Dim As Boolean mcpy(ub)
For i = 1 To ub
mcpy(i) = members(i)
Next
For i = 1 To ub
If mcpy(i) Then
If np = 5 And i > prime Then np = i
Dim As Uinteger n = prime * i
If n > limit Then Exit For there.
members(n) = False
End If
Next
If np < prime Then Exit While
cont += 1
Redim Preserve primes(cont)
primes(cont) = prime
prime = Iif(prime = 2, 3, np)
nlimit = min(stepLength * prime, limit)
Wend
Dim As Integer newPrimes(ub)
For i = 2 To ub
If members(i) Then newPrimes(i) = i
Next
cont = 0
For i = 1 To Ubound(primes)
If imprime Then Print primes(i);
cont += 1
Next
For i = 1 To ub
If newPrimes(i) Then
cont += 1
If imprime Then Print i;
End If
Next
If Not imprime Then Print !"\nNumber of primes up to "; limit; ":"; cont
End Sub
sieveOfPritchard(150, True)
sieveOfPritchard(1e6, False)
Sleep
- 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 1000000: 78498
True BASIC
FUNCTION min(a, b)
IF a < b THEN LET min = a ELSE LET min = b
END FUNCTION
SUB sieveofpritchard (limit)
DIM members(0)
MAT REDIM members(limit+1)
LET members(1) = 1
LET ub = UBOUND(members)
LET steplength = 1
LET prime = 2
LET rtlim = SQR(limit)
LET nlimit = 2
DIM primes(10)
LET cnt = 0
DIM mcpy(0)
MAT REDIM mcpy(ub)
DO WHILE prime <= rtlim
IF steplength < limit THEN
FOR w = 1 TO ub
IF members(w)<>0 THEN
LET n = w+steplength
DO WHILE n <= nlimit
LET members(n) = 1
LET n = n+steplength
LOOP
END IF
NEXT w
LET steplength = nlimit
END IF
LET np = 5
FOR i = 1 TO ub
LET mcpy(i) = members(i)
NEXT i
FOR i = 1 TO ub
IF mcpy(i)<>0 THEN
IF np = 5 AND i > prime THEN LET np = i
LET n = prime*i
IF n > limit THEN EXIT FOR
LET members(n) = 0
END IF
NEXT i
IF np < prime THEN EXIT DO
LET cnt = cnt+1
MAT REDIM primes(cnt)
LET primes(cnt) = prime
IF prime = 2 THEN LET prime = 3 ELSE LET prime = np
LET nlimit = min(steplength*prime, limit)
LOOP
DIM newprimes(0)
MAT REDIM newprimes(ub)
FOR i = 2 TO ub
IF members(i)<>0 THEN LET newprimes(i) = i
NEXT i
LET cnt = 0
FOR i = 1 TO UBOUND(primes)
PRINT primes(i);
LET cnt = cnt+1
NEXT i
FOR i = 1 TO ub
IF newprimes(i)<>0 THEN
PRINT i;
LET cnt = cnt+1
END IF
NEXT i
END SUB
CLEAR
CALL sieveofpritchard (150)
END
- Output:
Similar to FreeBASIC entry.
Yabasic
sieveOfPritchard(150, true)
sieveOfPritchard(1e6, false)
end
sub sieveOfPritchard(limit, imprime)
dim members(limit + 1)
members(1) = true
ub = arraysize(members(),1)
stepLength = 1
prime = 2
rtlim = sqrt(limit)
nlimit = 2
dim primes(1)
cont = 0
while prime <= rtlim
if stepLength < limit then
for w = 1 to ub
if members(w) then
n = w + stepLength
while n <= nlimit
members(n) = true
n = n + stepLength
wend
fi
next
stepLength = nlimit
fi
np = 5
dim mcpy(ub)
for i = 1 to ub
mcpy(i) = members(i)
next
for i = 1 to ub
if mcpy(i) then
if np = 5 and i > prime np = i
n = prime * i
if n > limit break
members(n) = false
fi
next
if np < prime break
cont = cont + 1
redim primes(cont)
primes(cont) = prime
if prime = 2 then prime = 3 else prime = np : fi
nlimit = min(stepLength * prime, limit)
wend
dim newPrimes(ub)
for i = 2 to ub
if members(i) newPrimes(i) = i
next
cont = 0
for i = 1 to arraysize(primes(),1)
if imprime print " ", primes(i);
cont = cont + 1
next
for i = 1 to ub
if newPrimes(i) then
cont = cont + 1
if imprime print " ", i;
fi
next
if not imprime then print : print "Number of primes up to ", limit, ": ", cont : fi
end sub
- Output:
Similar to FreeBASIC entry.
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++
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]
EasyLang
proc pritchard limit . primes[] .
len members[] limit
members[1] = 1
steplength = 1
prime = 2
rtlimit = sqrt limit
nlimit = 2
primes[] = [ ]
while prime <= rtlimit
if steplength < limit
for w to len members[]
if members[w] = 1
n = w + steplength
while n <= nlimit
members[n] = 1
n += steplength
.
.
.
steplength = nlimit
.
np = 5
mcpy[] = members[]
for w to nlimit
if mcpy[w] = 1
if np = 5 and w > prime
np = w
.
n = prime * w
if n > nlimit
break 1
.
members[n] = 0
.
.
if np < prime
break 1
.
primes[] &= prime
if prime = 2
prime = 3
else
prime = np
.
nlimit = lower (steplength * prime) limit
.
for i = 2 to len members[]
if members[i] = 1
primes[] &= i
.
.
.
pritchard 150 p[]
print p[]
- 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 ]
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.
Java
import java.util.ArrayList;
import java.util.BitSet;
import java.util.List;
public final class SieveOfPritchard {
public static void main(String[] args) {
System.out.println(sieveOfPritchard(150) + System.lineSeparator());
System.out.println("Number of primes up to 1,000,000 is " + sieveOfPritchard(1_000_000).size() + ".");
System.out.println();
final long start = System.currentTimeMillis();
System.out.print("Number of primes up to 100,000,000 is " + sieveOfPritchard(100_000_000).size());
final long finish = System.currentTimeMillis();
System.out.println(". Obtained in a time of " + ( (double) finish - start ) / 1_000 + " seconds.");
}
private static List<Integer> sieveOfPritchard(int limit) {
List<Integer> primes = new ArrayList<Integer>();
BitSet members = new BitSet(limit + 1);
members.set(1);
List<Integer> deletions = new ArrayList<Integer>();
final int rootLimit = (int) Math.sqrt(limit);
int nLimit = 2;
int stepLength = 1;
int prime = 2;
while ( prime <= rootLimit ) {
if ( stepLength < limit ) {
for ( int w = 1; w >= 0; w = members.nextSetBit(w + 1) ) {
int n = w + stepLength;
while ( n <= nLimit ) {
members.set(n);
n += stepLength;
}
}
stepLength = nLimit;
}
deletions.clear();
int nextPrime = 5;
for ( int w = 1; w < nLimit; w = members.nextSetBit(w + 1) ) {
if ( nextPrime == 5 && w > prime ) {
nextPrime = w;
}
final int n = prime * w;
if ( n > nLimit ) {
break;
}
deletions.add(n);
}
for ( int deletion : deletions ) {
members.clear(deletion);
}
if ( nextPrime < prime ) {
break;
}
primes.add(prime);
prime = ( prime == 2 ) ? 3 : nextPrime;
nLimit = (int) Math.min((long) stepLength * prime, limit);
}
if ( stepLength < limit ) {
for ( int w = 1; w >= 0; w = members.nextSetBit(w + 1) ) {
int n = w + stepLength;
while ( n <= limit ) {
members.set(n);
n += stepLength;
}
}
}
members.clear(1);
for ( int i = members.nextSetBit(0); i >= 0; i = members.nextSetBit(i + 1) ) {
primes.add(i);
};
return 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] Number of primes up to 1,000,000 is 78498. Number of primes up to 100,000,000 is 5761455. Obtained in a time of 0.648 seconds.
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
Nim
import std/[algorithm, math, sugar]
proc pritchard(limit: Natural): seq[int] =
## Pritchard sieve of primes up to "limit".
var members = newSeq[bool](limit + 1)
members[1] = true
var
stepLength = 1
prime = 2
rtlim = sqrt(limit.toFloat).int
nlimit = 2
primes: seq[int]
while prime <= rtlim:
if stepLength < limit:
for w in 1..members.high:
if members[w]:
var n = w + stepLength
while n <= nlimit:
members[n] = true
inc n, stepLength
stepLength = nlimit
var np = 5
var mcpy = members
for w in 1..members.high:
if mcpy[w]:
if np == 5 and w > prime:
np = w
let n = prime * w
if n > limit:
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.add prime
prime = if prime == 2: 3 else: np
nlimit = min(stepLength * prime, limit) # Advance wheel limit.
let newPrimes = collect:
for i in 2..members.high:
if members[i]: i
result = sorted(primes & newPrimes)
echo pritchard(150)
echo "Number of primes up to 1_000_000: ", pritchard(1_000_000).len
- 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
Pascal
A console program written in Delphi 7. It follows the algorithm in the Wikipedia article "Sieve of Pritchard".
program Pritchard_console;
{$APPTYPE CONSOLE}
uses
Math, SysUtils, Types;
// Function to return an array of all primes <= N
function PritchardSieve( const N : integer) : Types.TIntegerDynArray;
var
j, j_max, k, len, nrPrimes, p : integer;
marked : Types.TBooleanDynArray;
smallPrimes : Types.TIntegerDynArray; // i.e. primes <= sqrt( N)
spi : integer; // index into array smallPrimes
const
SP_STEP = 16; // step when extending dynamic array smallPrimes
begin
// Deal with trivial input
result := nil;
if (N <= 1) then exit;
// Initialize
SetLength( marked, N + 1); // 0..N for convenience; marked[0] is not used
marked[1] := true; // no other initialization of "marked" is needed
len := 1;
p := 2;
SetLength( smallPrimes, SP_STEP);
spi := 0;
while p*p <= N do begin
// Roll the wheel
if len < N then begin
j_max := Math.Min( p*len, N);
for j := len + 1 to j_max do marked[j] := marked[j - len];
len := j_max;
end;
// Unmark multiples of p
for k := len div p downto 1 do
if marked[k] then marked[p*k] := false;
// Store the prime p, extending the array if necessary
if spi = Length( smallPrimes) then
SetLength( smallPrimes, spi + SP_STEP);
smallPrimes[spi] := p;
inc(spi);
// Find the next prime p
if p = 2 then p := 3
else repeat inc(p) until (p > N) or marked[p];
// Condition p > N is a safety net; should always hit a marked value
Assert(p <= N);
end; // while
// Final roll, if needed. It is not needed if N >= 49. This is because
// 2 < 3^2, 2*3 < 5^2, 2*3*5 < 7^2, but thereafter 2*3*5*7 > 11^2, etc.
if len < N then
for j := len + 1 to N do marked[j] := marked[j - len];
// Remove 1 and put the small primes back
marked[1] := false;
for k := 0 to spi - 1 do marked[smallPrimes[k]] := true;
// Use the boolean array to return an array of prime integers
nrPrimes := 0;
for j := 2 to N do
if marked[j] then inc( nrPrimes);
SetLength( result, nrPrimes);
k := 0;
for j := 2 to N do
if marked[j] then begin result[k] := j; inc(k); end;
end;
// Main routine. User types the program name,
// optionally followed by the limit N (defaults to 150)
var
N, j : integer;
primes : Types.TIntegerDynArray;
begin
if ParamCount = 0 then N := 150
else N := SysUtils.StrToInt( ParamStr(1));
primes := PritchardSieve(N);
WriteLn( 'Number of primes = ', Length(primes));
for j := 0 to Length(primes) - 1 do begin
Write( ' ', primes[j]:4);
if j mod 10 = 9 then WriteLn;
end;
end.
- Output:
Number of primes = 35 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
Free Pascal
see Sieve_of_Eratosthenes#alternative_using_wheel
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