Sieve of Pritchard: Difference between revisions
Added Easylang
m (→{{header|AppleScript}}: Bug fix.) |
(Added Easylang) |
||
| (13 intermediate revisions by 8 users not shown) | |||
Line 2:
[[Category:Simple]]
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
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.)
Line 8:
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=GxgGMwLfTjE
;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 (
set
set
-- Get an extension list nominally expanding the wheel to the lesser of
-- its current circumference * currentPrime and the limit parameter.
set newCircumference
if (newCircumference > limit)
set o's extension to makeList(newCircumference - oldCircumference, mv)
-- Insert numbers that are oldCircumference added to
-- 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 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
end sieveOfPritchard
Line 115 ⟶ 101:
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)</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 326 ⟶ 646:
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>
Line 379 ⟶ 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 440 ⟶ 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 707 ⟶ 1,348:
{{libheader|Wren-sort}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
import "./fmt" for Fmt
| |||