1 files changed, 228 insertions, 0 deletions
diff --git a/appl/math/primes.b b/appl/math/primes.b
new file mode 100644
index 00000000..70bbc651
--- /dev/null
+++ b/appl/math/primes.b
@@ -0,0 +1,228 @@
+implement Primes;
+
+#
+# primes starting [ending]
+#
+# Subject to the Lucent Public License 1.02
+#
+
+include "draw.m";
+
+Primes: module
+{
+	init: fn(nil: ref Draw->Context, argl: list of string);
+};
+
+include "sys.m";
+	sys: Sys;
+include "math.m";
+	maths: Math;
+
+bigx: con 9.007199254740992e15;
+pt := array[] of {
+	2,
+	3,
+	5,
+	7,
+	11,
+	13,
+	17,
+	19,
+	23,
+	29,
+	31,
+	37,
+	41,
+	43,
+	47,
+	53,
+	59,
+	61,
+	67,
+	71,
+	73,
+	79,
+	83,
+	89,
+	97,
+	101,
+	103,
+	107,
+	109,
+	113,
+	127,
+	131,
+	137,
+	139,
+	149,
+	151,
+	157,
+	163,
+	167,
+	173,
+	179,
+	181,
+	191,
+	193,
+	197,
+	199,
+	211,
+	223,
+	227,
+	229,
+};
+wheel := array[] of {
+	10.0,
+	2.0,
+	4.0,
+	2.0,
+	4.0,
+	6.0,
+	2.0,
+	6.0,
+	4.0,
+	2.0,
+	4.0,
+	6.0,
+	6.0,
+	2.0,
+	6.0,
+	4.0,
+	2.0,
+	6.0,
+	4.0,
+	6.0,
+	8.0,
+	4.0,
+	2.0,
+	4.0,
+	2.0,
+	4.0,
+	8.0,
+	6.0,
+	4.0,
+	6.0,
+	2.0,
+	4.0,
+	6.0,
+	2.0,
+	6.0,
+	6.0,
+	4.0,
+	2.0,
+	4.0,
+	6.0,
+	2.0,
+	6.0,
+	4.0,
+	2.0,
+	4.0,
+	2.0,
+	10.0,
+	2.0,
+};
+BITS: con 8;
+TABLEN: con 1000;
+table := array[TABLEN] of byte;
+bittab := array[8] of {
+	byte 1,
+	byte 2,
+	byte 4,
+	byte 8,
+	byte 16,
+	byte 32,
+	byte 64,
+	byte 128,
+};
+
+init(nil: ref Draw->Context, args: list of string)
+{
+	sys = load Sys Sys->PATH;
+	maths = load Math Math->PATH;
+
+	if(len args <= 1){
+		sys->fprint(sys->fildes(2), "usage: primes starting [ending]\n");
+		raise "fail:usage";
+	}
+	args = tl args;
+	nn := real hd args;
+	limit := bigx;
+	if(tl args != nil){
+		limit = real hd tl args;
+		if(limit < nn)
+			exit;
+		if(limit > bigx)
+			ouch();
+	}
+	if(nn < 0.0 || nn > bigx)
+		ouch();
+	if(nn == 0.0)
+		nn = 1.0;
+	if(nn < 230.0){
+		for(i := 0; i < len pt; i++){
+			r := real pt[i];
+			if(r < nn)
+				continue;
+			if(r > limit)
+				exit;
+			sys->print("%d\n", pt[i]);
+			if(limit >= bigx)
+				exit;
+		}
+		nn = 230.0;
+	}
+	(t, nil) := maths->modf(nn/2.0);
+	nn = 2.0*real t+1.0;
+	for(;;){
+		# 
+		# clear the sieve table.
+		#  
+		for(i := 0; i < len table; i++)
+			table[i] = byte 0;
+		# 
+		# run the sieve
+		#  
+		v := maths->sqrt(nn+real (TABLEN*BITS));
+		mark(nn, 3);
+		mark(nn, 5);
+		mark(nn, 7);
+		i = 0;
+		for(k := 11.0; k <= v; k += wheel[i]){
+			mark(nn, int k);
+			i++;
+			if(i >= len wheel)
+				i = 0;
+		}
+		# 
+		# now get the primes from the table and print them
+		#  
+		for(i = 0; i < TABLEN*BITS; i += 2){
+			if(int table[i>>3]&int bittab[i&8r7])
+				continue;
+			temp := nn+real i;
+			if(temp > limit)
+				exit;
+			sys->print("%d\n", int temp);
+			if(limit >= bigx)
+				exit;
+		}
+		nn += real (TABLEN*BITS);
+	}
+}
+
+mark(nn: real, k: int)
+{
+	(it1, nil) := maths->modf(nn/real k);
+	j := int (real k*real it1-nn);
+	if(j < 0)
+		j += k;
+	for(; j < len table*BITS; j += k)
+		table[j>>3] |= bittab[j&8r7];
+}
+
+ouch()
+{
+	sys->fprint(sys->fildes(2), "primes: limits exceeded\n");
+	raise "fail:ouch";
+}
+