Sieve of Pritchard: Difference between revisions
Added Easylang
(Added Easylang) |
|||
| (18 intermediate revisions by 9 users not shown) | |||
Line 2:
[[Category:Simple]]
The [[wp:Sieve_of_Pritchard|Sieve of Pritchard]] is
Conceptually, it works by
For example, the second-order wheel has
[https://www.youtube.com/watch?v=
;Task:
Line 31:
if (limit < 2) then return {}
script o
property
property
end script
set {x, newCircumference, currentPrime, mv} to {0, 2, 1, missing value}
repeat until (currentPrime * currentPrime > limit)
set
set
-- Get an extension list nominally expanding the wheel to the lesser of
if (newCircumference > limit) then set newCircumference to limit
-- 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 i to x
repeat while
set i to i + 1
end repeat
set j to binarySearch((o's wheel's item i) * currentPrime, o's wheel, i, listLen)
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
end repeat
return
end sieveOfPritchard
Line 115 ⟶ 102:
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)</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}}==
Line 170 ⟶ 491:
=={{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;
Line 191 ⟶ 513:
<syntaxhighlight lang="cpp">/* Sieve of Pritchard in C++, as described at https://en.wikipedia.org/wiki/Sieve_of_Pritchard */
/*
/* 2 <= N <= 1000000000 */
/* (like the standard Sieve of Eratosthenes, this algorithm is not suitable for very large N due to memory requirements) */
Line 202 ⟶ 524:
#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
uint32_t i, x;
i = 0;
Line 233 ⟶ 550:
}
void Compress(uint32_t w[]
/* Removes deleted values in w[0..
uint32_t i, j;
j = 0;
for (i=1; i <=
if (!d[w[i]]) {
w[++j] = w[i];
}
}
if (
w_end = j;
} else {
for (uint32_t k=j+1; k <=
}
}
Line 258 ⟶ 575:
/* 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;
Line 274 ⟶ 591:
if (printPrimes) printf(" %d", p);
if (length < N) {
/* Extend W
Extend (w, w_end, length, std::min(p*length,N)
}
Delete(w, length, p, d, imaxf);
Compress(w
/* p = next(W, 1): */
/* k++ */
}
if (length < N) {
/* Extend full wheel W,length to N: */
Extend (w, w_end, length
}
/* gather remaining primes: */
for (uint32_t i=1; i <=
if (w[i] == 0 || d[w[i]]) continue;
if (printPrimes) printf(" %d", w[i]);
Line 301 ⟶ 615:
int main (int argc, char *argw[]) {
bool error = false;
bool printPrimes = false;
uint32_t N, nrPrimes, vBound;
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,
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>
{{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
35 primes up to 150 found in 0.
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>
Line 363 ⟶ 762:
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}}==
Line 424 ⟶ 919:
up to 1000000, added 186825, removed 108494, prime count 78498
</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}}==
Line 691 ⟶ 1,348:
{{libheader|Wren-sort}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
import "./fmt" for Fmt
| |||