diff options
Diffstat (limited to 'appl/math/sieve.b')
| -rw-r--r-- | appl/math/sieve.b | 220 |
1 files changed, 220 insertions, 0 deletions
diff --git a/appl/math/sieve.b b/appl/math/sieve.b new file mode 100644 index 00000000..25d4b0f1 --- /dev/null +++ b/appl/math/sieve.b @@ -0,0 +1,220 @@ +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; +} |