ref: eb0e026c241bd3afc6a10a6e2c7417de1b91f7e8
dir: /appl/math/sieve.b/
implement Sieve;
include "sys.m";
sys: Sys;
include "draw.m";
include "arg.m";
arg: Arg;
M: con 16*1024*1024;
N: con 8*M;
T: con 2*1024*1024;
limit := array[5] of { M, N, 2*N, 3*N, 15*(N/4) };
Sieve: module
{
init: fn(nil: ref Draw->Context, argv: list of string);
};
init(nil: ref Draw->Context, argv: list of string)
{
sys = load Sys Sys->PATH;
arg = load Arg Arg->PATH;
np := 0;
alg := 3;
arg->init(argv);
while((c := arg->opt()) != 0){
case (c){
'a' =>
alg = int arg->arg();
}
}
if(alg < 0 || alg > 4)
alg = 3;
lim := limit[alg];
argv = arg->argv();
if(argv != nil)
lim = int hd argv;
if(lim < 0 || lim > limit[alg])
lim = limit[alg];
if(lim < 6){
if(lim > 2){
sys->print("2\n");
np++;
}
if(lim > 3){
sys->print("3\n");
np++;
}
}
else{
case (alg){
0 => np = init0(lim);
1 => np = init1(lim);
2 => np = init2(lim);
3 => np = init3(lim);
4 => np = init4(lim);
}
}
sys->print("%d primes < %d\n", np, lim);
}
init0(lim: int): int
{
p := array[lim] of byte;
for(i := 0; i < lim; i++)
p[i] = byte 1;
p[0] = p[1] = byte 0;
np := 0;
for(i = 0; i < lim; i++){
if(p[i] == byte 1){
np++;
sys->print("%d\n", i);
for(j := i+i; j < lim; j += i)
p[j] = byte 0;
}
}
return np;
}
init1(lim: int): int
{
n := (lim+31)/32;
p := array[n] of int;
for(i := 0; i < n; i++)
p[i] = int 16rffffffff;
p[0] = int 16rfffffffc;
np := 0;
for(i = 0; i < lim; i++){
if(p[i>>5] & (1<<(i&31))){
np++;
sys->print("%d\n", i);
for(j := i+i; j < lim; j += i)
p[j>>5] &= ~(1<<(j&31));
}
}
return np;
}
init2(lim: int): int
{
n := ((lim+1)/2+31)/32;
p := array[n] of int;
for(i := 0; i < n; i++)
p[i] = int 16rffffffff;
p[0] = int 16rfffffffe;
np := 1;
sys->print("%d\n", 2);
for(i = 1; i < lim; i += 2){
k := (i-1)>>1;
if(p[k>>5] & (1<<(k&31))){
np++;
sys->print("%d\n", i);
inc := i+i;
for(j := i+i+i; j < lim; j += inc){
k = (j-1)>>1;
p[k>>5] &= ~(1<<(k&31));
}
}
}
return np;
}
init3(lim: int): int
{
n := ((lim+2)/3+31)/32;
p := array[n] of int;
for(i := 0; i < n; i++)
p[i] = int 16rffffffff;
p[0] = int 16rfffffffe;
np := 2;
sys->print("%d\n", 2);
sys->print("%d\n", 3);
d := 2;
for(i = 1; i < lim; i += d){
k := (i-1)/3;
if(p[k>>5] & (1<<(k&31))){
np++;
sys->print("%d\n", i);
inc := 6*i;
for(j := 5*i; j > 0 && j < lim; j += inc){
k = (j-1)/3;
p[k>>5] &= ~(1<<(k&31));
}
for(j = 7*i; j > 0 && j < lim; j += inc){
k = (j-1)/3;
p[k>>5] &= ~(1<<(k&31));
}
}
d = 6-d;
}
return np;
}
init4(lim: int): int
{
n := (4*((lim+14)/15)+31)/32;
p := array[n] of int;
for(i := 0; i < n; i++)
p[i] = int 16rffffffff;
p[0] = int 16rfffffffe;
np := 3;
sys->print("%d\n", 2);
sys->print("%d\n", 3);
sys->print("%d\n", 5);
m := -1;
d := array[8] of { 6, 4, 2, 4, 2, 4, 6, 2 };
for(i = 1; i < lim; i += d[m]){
k := (17*(i%30-1))/60+8*(i/30);
if(p[k>>5] & (1<<(k&31))){
np++;
sys->print("%d\n", i);
inc := 30*i;
for(j := 7*i; j > 0 && j < lim; j += inc){
k = (17*(j%30-1))/60+8*(j/30);
p[k>>5] &= ~(1<<(k&31));
}
for(j = 11*i; j > 0 && j < lim; j += inc){
k = (17*(j%30-1))/60+8*(j/30);
p[k>>5] &= ~(1<<(k&31));
}
for(j = 13*i; j > 0 && j < lim; j += inc){
k = (17*(j%30-1))/60+8*(j/30);
p[k>>5] &= ~(1<<(k&31));
}
for(j = 17*i; j > 0 && j < lim; j += inc){
k = (17*(j%30-1))/60+8*(j/30);
p[k>>5] &= ~(1<<(k&31));
}
for(j = 19*i; j > 0 && j < lim; j += inc){
k = (17*(j%30-1))/60+8*(j/30);
p[k>>5] &= ~(1<<(k&31));
}
for(j = 23*i; j > 0 && j < lim; j += inc){
k = (17*(j%30-1))/60+8*(j/30);
p[k>>5] &= ~(1<<(k&31));
}
for(j = 29*i; j > 0 && j < lim; j += inc){
k = (17*(j%30-1))/60+8*(j/30);
p[k>>5] &= ~(1<<(k&31));
}
for(j = 31*i; j > 0 && j < lim; j += inc){
k = (17*(j%30-1))/60+8*(j/30);
p[k>>5] &= ~(1<<(k&31));
}
}
m++;
if(m == 8)
m = 0;
}
return np;
}
init5(lim: int): int
{
# you must be joking
lim = 0;
return 0;
}