1implement Sieve; 2 3include "sys.m"; 4 sys: Sys; 5include "draw.m"; 6include "arg.m"; 7 arg: Arg; 8 9M: con 16*1024*1024; 10N: con 8*M; 11T: con 2*1024*1024; 12 13limit := array[5] of { M, N, 2*N, 3*N, 15*(N/4) }; 14 15Sieve: module 16{ 17 init: fn(nil: ref Draw->Context, argv: list of string); 18}; 19 20init(nil: ref Draw->Context, argv: list of string) 21{ 22 sys = load Sys Sys->PATH; 23 arg = load Arg Arg->PATH; 24 25 np := 0; 26 alg := 3; 27 arg->init(argv); 28 while((c := arg->opt()) != 0){ 29 case (c){ 30 'a' => 31 alg = int arg->arg(); 32 } 33 } 34 if(alg < 0 || alg > 4) 35 alg = 3; 36 lim := limit[alg]; 37 argv = arg->argv(); 38 if(argv != nil) 39 lim = int hd argv; 40 if(lim < 0 || lim > limit[alg]) 41 lim = limit[alg]; 42 if(lim < 6){ 43 if(lim > 2){ 44 sys->print("2\n"); 45 np++; 46 } 47 if(lim > 3){ 48 sys->print("3\n"); 49 np++; 50 } 51 } 52 else{ 53 case (alg){ 54 0 => np = init0(lim); 55 1 => np = init1(lim); 56 2 => np = init2(lim); 57 3 => np = init3(lim); 58 4 => np = init4(lim); 59 } 60 } 61 sys->print("%d primes < %d\n", np, lim); 62} 63 64init0(lim: int): int 65{ 66 p := array[lim] of byte; 67 for(i := 0; i < lim; i++) 68 p[i] = byte 1; 69 p[0] = p[1] = byte 0; 70 np := 0; 71 for(i = 0; i < lim; i++){ 72 if(p[i] == byte 1){ 73 np++; 74 sys->print("%d\n", i); 75 for(j := i+i; j < lim; j += i) 76 p[j] = byte 0; 77 } 78 } 79 return np; 80} 81 82init1(lim: int): int 83{ 84 n := (lim+31)/32; 85 p := array[n] of int; 86 for(i := 0; i < n; i++) 87 p[i] = int 16rffffffff; 88 p[0] = int 16rfffffffc; 89 np := 0; 90 for(i = 0; i < lim; i++){ 91 if(p[i>>5] & (1<<(i&31))){ 92 np++; 93 sys->print("%d\n", i); 94 for(j := i+i; j < lim; j += i) 95 p[j>>5] &= ~(1<<(j&31)); 96 } 97 } 98 return np; 99} 100 101init2(lim: int): int 102{ 103 n := ((lim+1)/2+31)/32; 104 p := array[n] of int; 105 for(i := 0; i < n; i++) 106 p[i] = int 16rffffffff; 107 p[0] = int 16rfffffffe; 108 np := 1; 109 sys->print("%d\n", 2); 110 for(i = 1; i < lim; i += 2){ 111 k := (i-1)>>1; 112 if(p[k>>5] & (1<<(k&31))){ 113 np++; 114 sys->print("%d\n", i); 115 inc := i+i; 116 for(j := i+i+i; j < lim; j += inc){ 117 k = (j-1)>>1; 118 p[k>>5] &= ~(1<<(k&31)); 119 } 120 } 121 } 122 return np; 123} 124 125init3(lim: int): int 126{ 127 n := ((lim+2)/3+31)/32; 128 p := array[n] of int; 129 for(i := 0; i < n; i++) 130 p[i] = int 16rffffffff; 131 p[0] = int 16rfffffffe; 132 np := 2; 133 sys->print("%d\n", 2); 134 sys->print("%d\n", 3); 135 d := 2; 136 for(i = 1; i < lim; i += d){ 137 k := (i-1)/3; 138 if(p[k>>5] & (1<<(k&31))){ 139 np++; 140 sys->print("%d\n", i); 141 inc := 6*i; 142 for(j := 5*i; j > 0 && j < lim; j += inc){ 143 k = (j-1)/3; 144 p[k>>5] &= ~(1<<(k&31)); 145 } 146 for(j = 7*i; j > 0 && j < lim; j += inc){ 147 k = (j-1)/3; 148 p[k>>5] &= ~(1<<(k&31)); 149 } 150 } 151 d = 6-d; 152 } 153 return np; 154} 155 156init4(lim: int): int 157{ 158 n := (4*((lim+14)/15)+31)/32; 159 p := array[n] of int; 160 for(i := 0; i < n; i++) 161 p[i] = int 16rffffffff; 162 p[0] = int 16rfffffffe; 163 np := 3; 164 sys->print("%d\n", 2); 165 sys->print("%d\n", 3); 166 sys->print("%d\n", 5); 167 m := -1; 168 d := array[8] of { 6, 4, 2, 4, 2, 4, 6, 2 }; 169 for(i = 1; i < lim; i += d[m]){ 170 k := (17*(i%30-1))/60+8*(i/30); 171 if(p[k>>5] & (1<<(k&31))){ 172 np++; 173 sys->print("%d\n", i); 174 inc := 30*i; 175 for(j := 7*i; j > 0 && j < lim; j += inc){ 176 k = (17*(j%30-1))/60+8*(j/30); 177 p[k>>5] &= ~(1<<(k&31)); 178 } 179 for(j = 11*i; j > 0 && j < lim; j += inc){ 180 k = (17*(j%30-1))/60+8*(j/30); 181 p[k>>5] &= ~(1<<(k&31)); 182 } 183 for(j = 13*i; j > 0 && j < lim; j += inc){ 184 k = (17*(j%30-1))/60+8*(j/30); 185 p[k>>5] &= ~(1<<(k&31)); 186 } 187 for(j = 17*i; j > 0 && j < lim; j += inc){ 188 k = (17*(j%30-1))/60+8*(j/30); 189 p[k>>5] &= ~(1<<(k&31)); 190 } 191 for(j = 19*i; j > 0 && j < lim; j += inc){ 192 k = (17*(j%30-1))/60+8*(j/30); 193 p[k>>5] &= ~(1<<(k&31)); 194 } 195 for(j = 23*i; j > 0 && j < lim; j += inc){ 196 k = (17*(j%30-1))/60+8*(j/30); 197 p[k>>5] &= ~(1<<(k&31)); 198 } 199 for(j = 29*i; j > 0 && j < lim; j += inc){ 200 k = (17*(j%30-1))/60+8*(j/30); 201 p[k>>5] &= ~(1<<(k&31)); 202 } 203 for(j = 31*i; j > 0 && j < lim; j += inc){ 204 k = (17*(j%30-1))/60+8*(j/30); 205 p[k>>5] &= ~(1<<(k&31)); 206 } 207 } 208 m++; 209 if(m == 8) 210 m = 0; 211 } 212 return np; 213} 214 215init5(lim: int): int 216{ 217 # you must be joking 218 lim = 0; 219 return 0; 220} 221