xref: /inferno-os/appl/math/primes.b (revision 37da2899f40661e3e9631e497da8dc59b971cbd0) 
1implement Primes;
2
3#
4# primes starting [ending]
5#
6# Subject to the Lucent Public License 1.02
7#
8
9include "draw.m";
10
11Primes: module
12{
13	init: fn(nil: ref Draw->Context, argl: list of string);
14};
15
16include "sys.m";
17	sys: Sys;
18include "math.m";
19	maths: Math;
20
21bigx: con 9.007199254740992e15;
22pt := array[] of {
23	2,
24	3,
25	5,
26	7,
27	11,
28	13,
29	17,
30	19,
31	23,
32	29,
33	31,
34	37,
35	41,
36	43,
37	47,
38	53,
39	59,
40	61,
41	67,
42	71,
43	73,
44	79,
45	83,
46	89,
47	97,
48	101,
49	103,
50	107,
51	109,
52	113,
53	127,
54	131,
55	137,
56	139,
57	149,
58	151,
59	157,
60	163,
61	167,
62	173,
63	179,
64	181,
65	191,
66	193,
67	197,
68	199,
69	211,
70	223,
71	227,
72	229,
73};
74wheel := array[] of {
75	10.0,
76	2.0,
77	4.0,
78	2.0,
79	4.0,
80	6.0,
81	2.0,
82	6.0,
83	4.0,
84	2.0,
85	4.0,
86	6.0,
87	6.0,
88	2.0,
89	6.0,
90	4.0,
91	2.0,
92	6.0,
93	4.0,
94	6.0,
95	8.0,
96	4.0,
97	2.0,
98	4.0,
99	2.0,
100	4.0,
101	8.0,
102	6.0,
103	4.0,
104	6.0,
105	2.0,
106	4.0,
107	6.0,
108	2.0,
109	6.0,
110	6.0,
111	4.0,
112	2.0,
113	4.0,
114	6.0,
115	2.0,
116	6.0,
117	4.0,
118	2.0,
119	4.0,
120	2.0,
121	10.0,
122	2.0,
123};
124BITS: con 8;
125TABLEN: con 1000;
126table := array[TABLEN] of byte;
127bittab := array[8] of {
128	byte 1,
129	byte 2,
130	byte 4,
131	byte 8,
132	byte 16,
133	byte 32,
134	byte 64,
135	byte 128,
136};
137
138init(nil: ref Draw->Context, args: list of string)
139{
140	sys = load Sys Sys->PATH;
141	maths = load Math Math->PATH;
142
143	if(len args <= 1){
144		sys->fprint(sys->fildes(2), "usage: primes starting [ending]\n");
145		raise "fail:usage";
146	}
147	args = tl args;
148	nn := real hd args;
149	limit := bigx;
150	if(tl args != nil){
151		limit = real hd tl args;
152		if(limit < nn)
153			exit;
154		if(limit > bigx)
155			ouch();
156	}
157	if(nn < 0.0 || nn > bigx)
158		ouch();
159	if(nn == 0.0)
160		nn = 1.0;
161	if(nn < 230.0){
162		for(i := 0; i < len pt; i++){
163			r := real pt[i];
164			if(r < nn)
165				continue;
166			if(r > limit)
167				exit;
168			sys->print("%d\n", pt[i]);
169			if(limit >= bigx)
170				exit;
171		}
172		nn = 230.0;
173	}
174	(t, nil) := maths->modf(nn/2.0);
175	nn = 2.0*real t+1.0;
176	for(;;){
177		#
178		# clear the sieve table.
179		#
180		for(i := 0; i < len table; i++)
181			table[i] = byte 0;
182		#
183		# run the sieve
184		#
185		v := maths->sqrt(nn+real (TABLEN*BITS));
186		mark(nn, 3);
187		mark(nn, 5);
188		mark(nn, 7);
189		i = 0;
190		for(k := 11.0; k <= v; k += wheel[i]){
191			mark(nn, int k);
192			i++;
193			if(i >= len wheel)
194				i = 0;
195		}
196		#
197		# now get the primes from the table and print them
198		#
199		for(i = 0; i < TABLEN*BITS; i += 2){
200			if(int table[i>>3]&int bittab[i&8r7])
201				continue;
202			temp := nn+real i;
203			if(temp > limit)
204				exit;
205			sys->print("%d\n", int temp);
206			if(limit >= bigx)
207				exit;
208		}
209		nn += real (TABLEN*BITS);
210	}
211}
212
213mark(nn: real, k: int)
214{
215	(it1, nil) := maths->modf(nn/real k);
216	j := int (real k*real it1-nn);
217	if(j < 0)
218		j += k;
219	for(; j < len table*BITS; j += k)
220		table[j>>3] |= bittab[j&8r7];
221}
222
223ouch()
224{
225	sys->fprint(sys->fildes(2), "primes: limits exceeded\n");
226	raise "fail:ouch";
227}
228
229