code: purgatorio

ref: eb0e026c241bd3afc6a10a6e2c7417de1b91f7e8
dir: /appl/math/primes.b/

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";
}