| ' Sieve of Eratosthenes, Sifting thru 1 Million Prime Numbers
' BYTE, October 1985, Special IBM Edition, page 244
' Turbo Pascal 1.1, Macintosh Version
program Sieve_of_Eratosthenes;
uses
memtypes, quickdraw, osintf, toolintf; { Macintosh Specific }
label
branchA, branchB, branchC, branchD;
type
arraytype1 = array[1..10000] of integer;
arraytype2 = array[1..550] of integer;
var
timezero, totaltime: integer;
c, i, j, k, n, p: integer;
x, y, z: longint;
dummy: eventrecord;
array1: arraytype1;
array2: arraytype2;
array3: arraytype2;
begin
writeln('Sifting Thru 1 Million Prime Numbers');
timezero := tickcount; { 60 Ticks = 1 Second }
p := 2;
{ writeln(p:10); }
c := 2;
for i := 1 to 10000 do
array1[i] := 1;
for i := 1 to 10000 do
begin
if array1[i] = 0 then
begin
array1[i] := 1;
goto branchA
end;
p := i + i + 1;
{ writeln(p:10); }
if p > 4000 then
goto branchB;
k := i + p;
while k <= 10000 do
begin
array1[k] := 0;
k := k + p
end;
array3[c] := p;
array2[c] := k;
branchB:
c := c + 1;
branchA:
end;
x := 1;
z := c;
for n := 1 to 799 do
begin
for j := 2 to 550 do
begin
p := array3[j];
k := array2[j] - 10000;
while k <= 10000 do
begin
array1[k] := 0;
k := k + p
end;
array2[j] := k
end;
x := x + 20000;
for i := 1 to 10000 do
begin
if array1[i] = 0 then
begin
array1[i] := 1;
goto branchC
end;
y := x + i + i;
{ writeln(y:10); }
if z = 1000000 then { try 2^31-1 or 2,147,483,647 }
goto branchD;
z := z + 1;
branchC:
end;
end;
branchD:
totaltime := (tickcount - timezero);
writeln(y, ' is the ', z, 'th prime number.');
writeln('Elapse time: ', totaltime/60, ' seconds.');
while not oseventavail(mdownmask + keydownmask, dummy) do
{ Macintosh Specific }
end.
|
| Re .0, .1
To find the primes to some number like 2^32, you only have to apply to the
sieve the multiples of primes up to the squareroot of your limit, which
for 2^32 would be 2^16.
If I understand the program in .1, although the sieve represents 10,000
odd numbers, the first loop appears to only save the primes less than
4000 or so,
[p is i+i+1, and if p > 4000 it appears to quit saving primes. The
second loop appears to assume that there are only 550 primes whose
multiples need to be applied to the sieve, so it had better be right.]
It then begins to use the primes saved to find the rest of
the primes. Which would reliably only allow it to find the primes less
than 4000*4000 or 16 million or so. And I think the 1 millionth prime is
in the 14-15 million range so the program appears to work for the 1st
million primes, but won't work if you go past 16 million.
Anyway, to sieve up to greater than 2^32, you have to calculate all the
primes up to the square root of your limit and save them. .1, uses the
descriptive names of array2 and array3. array2 as the sieve position
to begin marking multiples on the next pass, and array3 the distance to
step by while marking multiples. You may find that the distance you
step by becomes larger than the sieve size, and you'll have to worry about
handling the case that every prime won't have a multiple in every pass.
One way would be to keep four arrays, two used for number of passes and step
distance of each prime, and two to report how far to go to make the wipe
out the next multiple. Or it'd be easier if the sieve size is large
enough so you don't have to do this.
There's probably not more than 10,000 primes between 1 and 2^16, but if
you want to go higher than 2^32, you'll have to figure out how many to
save.
Any particular prime you find is the number of passes over the sieve *
the sieve size + the index into the sieve it is found at.
So if your sieve has 8K entries (2^13) then after 2^40/2^13 or 2^27 passes
you would have counted the number primes less than 2^40. (It'd only take
a few months). Notice this won't require any use of variables greater than
32 bits, you just think of them as having two parts. The number of passes
it was found at, and the index into the sieve. So, in .1 change x to number
of passes, and print the primes out as something like x*8K + index.
Simon
|