xref: /inferno-os/appl/math/powers.b (revision 37da2899f40661e3e9631e497da8dc59b971cbd0)

implement Powers;

include "sys.m";
	sys: Sys;
include "draw.m";
include "arg.m";
include "lock.m";
	lockm: Lock;
	Semaphore: import lockm;

Powers: module
{
	init: fn(nil: ref Draw->Context, nil: list of string);
};

MAXNODES: con (1<<20)/4;

verbose: int;

# Doing
# 	powers -p 3
# gives
# 	[2] 1729 = 1**3 + 12**3 = 9**3 + 10**3
# 	[2] 4104 = 2**3 + 16**3 = 9**3 + 15**3

# ie 1729 can be written in two ways as the sum of 2 cubes as can 4104.

# The options are

# -p	the power to use - default 2
# -n	the number of powers summed - default 2
# -f	the minimum number of ways found before reporting it - default 2
# -l	the least number to consider - default 0
# -m	the greatest number to consider - default 8192

# Thus
# 	pow -p 4 -n 3 -f 3 -l 0 -m 1000000
# gives
# 	[3] 811538 = 12**4 + 17**4 + 29**4 = 7**4 + 21**4 + 28**4 = 4**4 + 23**4 + 27**4

# ie fourth powers, 3 in each sum, minimum of 3 representations, numbers from 0-1000000.

# [2] 25
# [3] 325
# [4] 1105
# [5] 4225
# [6] 5525
# [7] 203125
# [8] 27625
# [9] 71825
# [10] 138125
# [11] 2640625
# [12] 160225
# [13] 17850625
# [14] 1221025
# [15] 1795625
# [16] 801125
# [18] 2082925
# [20] 4005625
# [23] 30525625
# [24] 5928325
# [32] 29641625

# [24] 5928325 = 63**2 + 2434**2 = 94**2 + 2433**2 = 207**2 + 2426**2 = 294**2 + 2417**2 = 310**2 + 2415**2 = 465**2 + 2390**2 = 490**2 + 2385**2 = 591**2 + 2362**2 = 690**2 + 2335**2 = 742**2 + 2319**2 = 849**2 + 2282**2 = 878**2 + 2271**2 = 959**2 + 2238**2 = 1039**2 + 2202**2 = 1062**2 + 2191**2 = 1201**2 + 2118**2 = 1215**2 + 2110**2 = 1290**2 + 2065**2 = 1410**2 + 1985**2 = 1454**2 + 1953**2 = 1535**2 + 1890**2 = 1614**2 + 1823**2 = 1633**2 + 1806**2 = 1697**2 + 1746**2

# [32] 29641625 = 67**2 + 5444**2 = 124**2 + 5443**2 = 284**2 + 5437**2 = 320**2 + 5435**2 = 515**2 + 5420**2 = 584**2 + 5413**2 = 835**2 + 5380**2 = 955**2 + 5360**2 = 1180**2 + 5315**2 = 1405**2 + 5260**2 = 1460**2 + 5245**2 = 1648**2 + 5189**2 = 1795**2 + 5140**2 = 1829**2 + 5128**2 = 1979**2 + 5072**2 = 2012**2 + 5059**2 = 2032**2 + 5051**2 = 2245**2 + 4960**2 = 2308**2 + 4931**2 = 2452**2 + 4861**2 = 2560**2 + 4805**2 = 2621**2 + 4772**2 = 2840**2 + 4645**2 = 3005**2 + 4540**2 = 3035**2 + 4520**2 = 3320**2 + 4315**2 = 3365**2 + 4280**2 = 3517**2 + 4156**2 = 3544**2 + 4133**2 = 3664**2 + 4027**2 = 3715**2 + 3980**2 = 3803**2 + 3896**2

# [2] 1729 = 1**3 + 12**3 = 9**3 + 10**3
# [2] 4104 = 2**3 + 16**3 = 9**3 + 15**3
# [3] 87539319 = 167**3 + 436**3 = 228**3 + 423**3 = 255**3 + 414**3

# [2] 635318657 = 59**4 + 158**4 = 133**4 + 134**4
# [2] 3262811042 = 7**4 + 239**4 = 157**4 + 227**4
# [2] 8657437697 = 193**4 + 292**4 = 256**4 + 257**4
# [2] 68899596497 = 271**4 + 502**4 = 298**4 + 497**4
# [2] 86409838577 = 103**4 + 542**4 = 359**4 + 514**4
# [2] 160961094577 = 222**4 + 631**4 = 503**4 + 558**4
# [2] 2094447251857 = 76**4 + 1203**4 = 653**4 + 1176**4
# [2] 4231525221377 = 878**4 + 1381**4 = 997**4 + 1342**4
# [2] 26033514998417 = 1324**4 + 2189**4 = 1784**4 + 1997**4
# [2] 37860330087137 = 1042**4 + 2461**4 = 2026**4 + 2141**4
# [2] 61206381799697 = 248**4 + 2797**4 = 2131**4 + 2524**4
# [2] 76773963505537 = 1034**4 + 2949**4 = 1797**4 + 2854**4
# [2] 109737827061041 = 1577**4 + 3190**4 = 2345**4 + 2986**4
# [2] 155974778565937 = 1623**4 + 3494**4 = 2338**4 + 3351**4
# [2] 156700232476402 = 661**4 + 3537**4 = 2767**4 + 3147**4
# [2] 621194785437217 = 2694**4 + 4883**4 = 3966**4 + 4397**4
# [2] 652057426144337 = 604**4 + 5053**4 = 1283**4 + 5048**4
# [2] 680914892583617 = 3364**4 + 4849**4 = 4288**4 + 4303**4
# [2] 1438141494155441 = 2027**4 + 6140**4 = 4840**4 + 5461**4
# [2] 1919423464573697 = 274**4 + 6619**4 = 5093**4 + 5942**4
# [2] 2089568089060657 = 498**4 + 6761**4 = 5222**4 + 6057**4
# [2] 2105144161376801 = 2707**4 + 6730**4 = 3070**4 + 6701**4
# [2] 3263864585622562 = 1259**4 + 7557**4 = 4661**4 + 7269**4
# [2] 4063780581008977 = 5181**4 + 7604**4 = 6336**4 + 7037**4
# [2] 6315669699408737 = 1657**4 + 8912**4 = 7432**4 + 7559**4
# [2] 6884827518602786 = 635**4 + 9109**4 = 3391**4 + 9065**4
# [2] 7191538859126257 = 4903**4 + 9018**4 = 6842**4 + 8409**4
# [2] 7331928977565937 = 1104**4 + 9253**4 = 5403**4 + 8972**4
# [2] 7362748995747617 = 5098**4 + 9043**4 = 6742**4 + 8531**4
# [2] 7446891977980337 = 1142**4 + 9289**4 = 4946**4 + 9097**4
# [2] 7532132844821777 = 173**4 + 9316**4 = 4408**4 + 9197**4
# [2] 7985644522300177 = 6262**4 + 8961**4 = 7234**4 + 8511**4

# 5, 6, 7, 8, 9, 10, 11 none

Btree: adt{
	sum: big;
	left: cyclic ref Btree;
	right: cyclic ref Btree;
};

Dtree: adt{
	sum: big;
	freq: int;
	lst: list of array of int;
	left: cyclic ref Dtree;
	right: cyclic ref Dtree;
};

nCr(n: int, r: int): int
{
	if(r > n-r)
		r = n-r;

	# f := g := 1;
	# for(i := 0; i < r; i++){
	# 	f *= n-i;
	# 	g *= i+1;
	# }
	# return f/g;

	num := array[r] of int;
	den := array[r] of int;
	for(i := 0; i < r; i++){
		num[i] = n-i;
		den[i] = i+1;
	}
	for(i = 0; i < r; i++){
		for(j := 0; den[i] != 1; j++){
			if(num[j] == 1)
				continue;
			k := hcf(num[j], den[i]);
			if(k != 1){
				num[j] /= k;
				den[i] /= k;
			}
		}
	}
	f := 1;
	for(i = 0; i < r; i++)
		f *= num[i];
	return f;
}

nHr(n: int, r: int): int
{
	if(n == 0)
		return 0;
	return nCr(n+r-1, r);
}

nSr(n: int, i: int, j: int): int
{
	return nHr(j, n)-nHr(i, n);
	# s := 0;
	# for(k := i; k < j; k++)
	# 	s += nHr(k+1, n-1);
	# return s;
}

nSrmax(n: int, i: int, m: int): int
{
	s := 0;
	for(k := i; ; k++){
		s += nHr(k+1, n-1);
		if(s > m)
			break;
	}
	if(k == i)
		return i+1;
	return k;
}

kth(c: array of int, n: int, i: int, j: int, k: int)
{
	l, u: int;

	m := nSr(n, i, j);
	if(k < 0)
		k = 0;
	if(k >= m)
		k = m-1;
	p := 0;
	for(q := 0; q < n; q++){
		if(q == 0){
			l = i;
			u = j-1;
		}
		else{
			l = 0;
			u = c[q-1];
		}
		for(x := l; x <= u; x++){
			m = nHr(x+1, n-q-1);
			p += m;
			if(p > k){
				p -= m;
				break;
			}
		}
		c[q] = x;
	}
}

pos(c: array of int, n: int): int
{
	p := 0;
	for(q := 0; q < n; q++)
		p += nSr(n-q, 0, c[q]);
	return p;
}

min(c: array of int, n: int, p: int): big
{
	s := big(0);
	for(i := 0; i < n; i++)
		s += big(c[i])**p;
	m := s;
	for(i = n-1; i > 0; i--){
		s -= big(c[i])**p;
		s -= big(c[i-1])**p;
		c[i]--;
		c[i-1]++;
		s += big(c[i-1])**p;
		if(s < m)
			m = s;
	}
	c[0]--;
	c[n-1]++;
	# m--;
	return m;
}

hcf(a, b: int): int
{
	if(b == 0)
		return a;
	for(;;){
		if(a == 0)
			break;
		if(a < b)
			(a, b) = (b, a);
		a %= b;
		# a -= (a/b)*b;
	}
	return b;
}

gcd(l: list of array of int): int
{
	g := (hd l)[0];
	for(; l != nil; l = tl l){
		d := hd l;
		n := len d;
		for(i := 0; i < n; i++)
			g = hcf(d[i], g);
	}
	return g;
}

adddup(s: big, root: ref Dtree): int
{
	n, p, lp: ref Dtree;

	p = root;
	while(p != nil){
		if(s == p.sum)
			return ++p.freq;
		lp = p;
		if(s < p.sum)
			p = p.left;
		else
			p = p.right;
	}
	n = ref Dtree(s, 2, nil, nil, nil);
	if(s < lp.sum)
		lp.left = n;
	else
		lp.right = n;
	return n.freq;
}

cp(c: array of int): array of int
{
	n := len c;
	m := 0;
	for(i := 0; i < n; i++)
		if(c[i] != 0)
			m++;
	nc := array[m] of int;
	nc[0: ] = c[0: m];
	return nc;
}

finddup(s: big, c: array of int, root: ref Dtree, f: int)
{
	p: ref Dtree;

	p = root;
	while(p != nil){
		if(s == p.sum){
			if(p.freq >= f)
				p.lst = cp(c) :: p.lst;
			return;
		}
		if(s < p.sum)
			p = p.left;
		else
			p = p.right;
	}
}

printdup(p: ref Dtree, pow: int, ix: int)
{
	if(p == nil)
		return;
	printdup(p.left, pow, ix);
	if((l := p.lst) != nil){
		if(gcd(l) == 1){
			min1 := min2 := 16r7fffffff;
			for(; l != nil; l = tl l){
				n := len hd l;
				if(n < min1){
					min2 = min1;
					min1 = n;
				}
				else if(n < min2)
					min2 = n;
			}
			i := min1+min2-pow;
			if(i <= ix){
				sys->print("[%d, %d] %bd", i, p.freq, p.sum);
				for(l = p.lst; l != nil; l = tl l){
					d := hd l;
					n := len d;
					sys->print(" = ");
					for(j := n-1; j >= 0; j--){
						sys->print("%d**%d", d[j], pow);
						if(j > 0)
							sys->print(" + ");
					}
				}
				sys->print("\n");
				if(i < 0){
					sys->print("****************\n");
					exit;
				}
			}
		}
	}
	printdup(p.right, pow, ix);
}

addsum(s: big, root: ref Btree, root1: ref Dtree): int
{
	n, p, lp: ref Btree;

	p = root;
	while(p != nil){
		if(s == p.sum)
			return adddup(s, root1);
		lp = p;
		if(s < p.sum)
			p = p.left;
		else
			p = p.right;
	}
	n = ref Btree(s, nil, nil);
	if(s < lp.sum)
		lp.left = n;
	else
		lp.right = n;
	return 1;
}

oiroot(x: big, p: int): int
{
	for(i := 0; ; i++){
		n := big(i)**p;
		if(n > x)
			break;
	}
	return i-1;
}

iroot(x: big, p: int): int
{
	m: big;

	if(x == big(0) || x == big(1))
		return int x;
	v := x;
	n := 0;
	for(i := 32; i > 0; i >>= 1){
		m = ((big(1)<<i)-big(1))<<i;
		if((v&m) != big(0)){
			n += i;
			v >>= i;
		}
	}
	a := big(1) << (n/p);
	b := a<<1;
	while(a < b){
		m = (a+b+big(1))/big(2);
		y := m**p;
		if(y > x)
			b = m-big(1);
		else if(y < x)
			a = m;
		else
			a = b = m;
	}
	if(a**p <= x && (a+big(1))**p > x)
		;
	else{
		sys->print("fatal: %bd %d -> %bd\n", x, p, a);
		exit;
	}
	return int a;
}

initval(c: array of int, n: int, p: int, v: int): big
{
	for(i := 0; i < n; i++)
		c[i] = 0;
	c[0] = v;
	return big(v)**p;
}

nxtval(c: array of int, n: int, p: int, s: big): big
{
	for(k := n-1; k >= 0; k--){
		s -= big(c[k])**p;
		c[k]++;
		if(k == 0){
			s += big(c[k])**p;
			break;
		}
		else{
			if(c[k] <= c[k-1]){
				s += big(c[k])**p;
				break;
			}
			c[k] = 0;
		}
	}
	return s;
}

powers(p: int, n: int, f: int, ix: int, lim0: big, lim: big, ch: chan of int, lock: ref Semaphore)
{
	root := ref Btree(big(-1), nil, nil);
	root1 := ref Dtree(big(-1), 0, nil, nil, nil);

	min := max := lim0;

	c := array[n] of int;

	for(;;){
		imin := iroot((min+big(n-1))/big(n), p);
		imax := nSrmax(n, imin, MAXNODES);
		max = big(imax)**p - big(1);
		while(max <= min){	# could do better
			imax++;
			max = big(imax)**p - big(1);
		}
		if(max > lim){
			max = lim;
			imax = iroot(max, p)+1;
		}

		if(verbose)
			sys->print("searching in %d-%d(%bd-%bd)\n", imin, imax, min, max);

		m := mm := 0;
		maxf := 0;
		s := initval(c, n, p, imin);
		for(;;){
			mm++;
			if(s >= min && s < max){
				fr := addsum(s, root, root1);
				if(fr > maxf)
					maxf = fr;
				m++;
			}
			s = nxtval(c, n, p, s);
			if(c[0] == imax)
				break;
		}

		root.left = root.right = nil;

		if(maxf >= f){
			if(verbose)
				sys->print("finding duplicates\n");

			s = initval(c, n, p, imin);
			for(;;){
				if(s >= min && s < max)
					finddup(s, c, root1, f);
				s = nxtval(c, n, p, s);
				if(c[0] == imax)
					break;
			}

			if(lock != nil)
				lock.obtain();
			printdup(root1, p, ix);
			if(lock != nil)
				lock.release();

			root1.left = root1.right = nil;
		}

		if(verbose)
			sys->print("%d(%d) nodes searched\n", m, mm);

		if(mm != nSr(n, imin, imax)){
			sys->print("**fatal**\n");
			exit;
		}

		min = max;
		if(min >= lim)
			break;
	}
	if(ch != nil)
		ch <-= 0;
}

usage()
{
	sys->print("usage: powers -p power -n number -f frequency -i index -l minimum -m maximum -s procs -v\n");
	exit;
}

partition(p: int, n: int, l: big, m: big, s: int): array of big
{
	a := array[s+1] of big;
	a[0] = big(iroot(l, p))**n;
	a[s] = (big(iroot(m, p))+big(1))**n;
	nn := a[s]-a[0];
	q := nn/big(s);
	r := nn-q*big(s);
	t := big(0);
	lb := a[0];
	for(i := 0; i < s; i++){
		ub := lb+q;
		t += r;
		if(t >= big(s)){
			ub++;
			t -= big(s);
		}
		a[i+1] = ub;
		lb = ub;
	}
	if(a[s] != a[0]+nn){
		sys->print("fatal: a[s]\n");
		exit;
	}
	for(i = 0; i < s; i++){
		# sys->print("%bd %bd\n", a[i], a[i]**p);
		a[i] = big(iroot(a[i], n))**p;
	}
	a[0] = l;
	a[s] = m;
	while(a[0] >= a[1]){
		a[1] = a[0];
		a = a[1: ];
		--s;
	}
	while(a[s] <= a[s-1]){
		a[s-1] = a[s];
		a = a[0: s];
		--s;
	}
	return a;
}

init(nil: ref Draw->Context, args: list of string)
{
	sys = load Sys Sys->PATH;
	arg := load Arg Arg->PATH;
	lockm = load Lock Lock->PATH;

	lockm->init();
	lock := Semaphore.new();

	p := n := f := 2;
	ix := 1<<30;
	l := m := big(0);
	s := 1;

	arg->init(args);
	while((c := arg->opt()) != 0){
		case c {
			'p' =>
				p = int arg->arg();
			'n' =>
				n = int arg->arg();
			'f' =>
				f = int arg->arg();
			'i' =>
				ix = int arg->arg();
			'l' =>
				l = big(arg->arg());
			'm' =>
				m = big(arg->arg())+big(1);
			's' =>
				s = int arg->arg();
			'v' =>
				verbose = 1;
			* =>
				usage();
		}
	}
	if(arg->argv() != nil)
		usage();

	if(p < 2){
		p = 2;
		sys->print("setting p = %d\n", p);
	}
	if(n < 2){
		n = 2;
		sys->print("setting n = %d\n", n);
	}
	if(f < 2){
		f = 2;
		sys->print("setting f = %d\n", f);
	}
	if(l < big(0)){
		l = big(0);
		sys->print("setting l = %bd\n", l);
	}
	if(m <= big(0)){
		m = big((1<<13)+1);
		sys->print("setting m = %bd\n", m-big(1));
	}
	if(l >= m)
		exit;

	if(s <= 1)
		powers(p, n, f, ix, l, m, nil, nil);
	else{
		nproc := 0;
		ch := chan of int;
		a := partition(p, n, l, m, s);
		lb := a[0];
		for(i := 0; i < s; i++){
			ub := a[i+1];
			if(lb < ub){
				nproc++;
				spawn powers(p, n, f, ix, lb, ub, ch, lock);
			}
			lb = ub;
		}
		for( ; nproc != 0; nproc--)
			<- ch;
	}
}