summaryrefslogtreecommitdiff
path: root/appl/math
diff options
context:
space:
mode:
Diffstat (limited to 'appl/math')
-rw-r--r--appl/math/ack.b43
-rw-r--r--appl/math/crackerbarrel.b133
-rw-r--r--appl/math/doc.txt48
-rw-r--r--appl/math/factor.b180
-rw-r--r--appl/math/ffts.b639
-rw-r--r--appl/math/fibonacci.b59
-rw-r--r--appl/math/fit.b264
-rw-r--r--appl/math/genprimes.b64
-rw-r--r--appl/math/geodesy.b849
-rw-r--r--appl/math/gr.b557
-rw-r--r--appl/math/graph0.b84
-rw-r--r--appl/math/hist0.b99
-rw-r--r--appl/math/linalg.b197
-rw-r--r--appl/math/linbench.b197
-rw-r--r--appl/math/mersenne.b99
-rw-r--r--appl/math/mkfile43
-rw-r--r--appl/math/parts.b125
-rw-r--r--appl/math/perms.b131
-rw-r--r--appl/math/pi.b405
-rw-r--r--appl/math/polyfill.b440
-rw-r--r--appl/math/polyhedra.b195
-rw-r--r--appl/math/powers.b672
-rw-r--r--appl/math/primes.b228
-rw-r--r--appl/math/sieve.b220
24 files changed, 5971 insertions, 0 deletions
diff --git a/appl/math/ack.b b/appl/math/ack.b
new file mode 100644
index 00000000..cbed6fdd
--- /dev/null
+++ b/appl/math/ack.b
@@ -0,0 +1,43 @@
+implement Ackermann;
+
+include "sys.m";
+ sys: Sys;
+include "draw.m";
+ draw: Draw;
+
+Ackermann: module
+{
+ init: fn(ctxt: ref Draw->Context, argv: list of string);
+};
+
+init(nil: ref Draw->Context, argv: list of string)
+{
+ sys = load Sys Sys->PATH;
+ draw = load Draw Draw->PATH;
+ argv = tl argv; # remove program name
+ m := n := 0;
+ if(argv != nil){
+ m = int hd argv;
+ argv = tl argv;
+ }
+ if(m < 0)
+ m = 0;
+ if(argv != nil)
+ n = int hd argv;
+ if(n < 0)
+ n = 0;
+ t0 := sys->millisec();
+ a := ack(m, n);
+ t1 := sys->millisec();
+ sys->print("A(%d, %d) = %d (t = %d ms)\n", m, n, a, t1-t0);
+}
+
+ack(m, n: int) : int
+{
+ if(m == 0)
+ return n+1;
+ else if(n == 0)
+ return ack(m-1, 1);
+ else
+ return ack(m-1, ack(m, n-1));
+}
diff --git a/appl/math/crackerbarrel.b b/appl/math/crackerbarrel.b
new file mode 100644
index 00000000..e5e1357b
--- /dev/null
+++ b/appl/math/crackerbarrel.b
@@ -0,0 +1,133 @@
+implement CBPuzzle;
+
+# Cracker Barrel Puzzle
+#
+# Holes are drilled in a triangular arrangement into which all but one
+# are seated pegs. A 6th order puzzle appears in the diagram below.
+# Note, the hole in the lower left corner of the triangle is empty.
+#
+# V
+# V V
+# V V V
+# V V V V
+# V V V V V
+# O V V V V V
+#
+# Pegs are moved by jumping over a neighboring peg thereby removing the
+# jumped peg. A peg can only be moved if a neighboring hole contains a
+# peg and the hole on the other side of the neighbor is empty. The last
+# peg cannot be removed.
+#
+# The object is to remove as many pegs as possible.
+
+include "sys.m";
+ sys: Sys;
+include "draw.m";
+
+CBPuzzle: module {
+ init: fn(nil: ref Draw->Context, args: list of string);
+};
+
+ORDER: con 6;
+
+Move: adt {
+ x, y: int;
+};
+
+valid:= array[] of {Move (1,0), (0,1), (-1,1), (-1,0), (0,-1), (1,-1)};
+
+board:= array[ORDER*ORDER] of int;
+pegs, minpegs: int;
+
+puzzle(): int
+{
+ if (pegs < minpegs)
+ minpegs = pegs;
+
+ if (pegs == 1)
+ return 1;
+
+ # Check each row of puzzle
+ for (r := 0; r < ORDER; r += 1)
+ # Check each column
+ for (c := 0; c < ORDER-r; c += 1) {
+ fromx := r*ORDER + c;
+ # Is a peg in this hole?
+ if (board[fromx])
+ # Check valid moves from this hole
+ for (m := 0; m < len valid; m += 1) {
+ tor := r + 2*valid[m].y;
+ toc := c + 2*valid[m].x;
+
+ # Is new location still on the board?
+ if (tor + toc < ORDER && tor >= 0 && toc >= 0) {
+ jumpr := r + valid[m].y;
+ jumpc := c + valid[m].x;
+ jumpx := jumpr*ORDER + jumpc;
+
+ # Is neighboring hole occupied?
+ if (board[jumpx]) {
+ # Is new location empty?
+ tox := tor*ORDER + toc;
+
+ if (! board[tox]) {
+ # Jump neighboring hole
+ board[fromx] = 0;
+ board[jumpx] = 0;
+ board[tox] = 1;
+ pegs -= 1;
+
+ # Try solving puzzle from here
+ if (puzzle()) {
+ #sys->print("(%d,%d) - (%d,%d)\n", r, c, tor, toc);
+ return 1;
+ }
+ # Dead end, put pegs back and try another move
+ board[fromx] = 1;
+ board[jumpx] = 1;
+ board[tox] = 0;
+ pegs += 1;
+ } # empty location
+ } # occupied neighbor
+ } # still on board
+ } # valid moves
+ }
+ return 0;
+}
+
+solve(): int
+{
+ minpegs = pegs = (ORDER+1)*ORDER/2 - 1;
+
+ # Put pegs on board
+ for (r := 0; r < ORDER; r += 1)
+ for (c := 0; c < ORDER - r; c += 1)
+ board[r*ORDER + c] = 1;
+
+ # Remove one peg
+ board[0] = 0;
+
+ return puzzle();
+}
+
+init(nil: ref Draw->Context, args: list of string)
+{
+ sys = load Sys Sys->PATH;
+
+ TRIALS: int;
+ if (len args < 2)
+ TRIALS = 1;
+ else
+ TRIALS = int hd tl args;
+
+ start := sys->millisec();
+ for (trials := 0; trials < TRIALS; trials += 1)
+ solved := solve();
+ end := sys->millisec();
+
+ sys->print("%d ms\n", end - start);
+
+ if (! solved)
+ sys->print("No solution\n");
+ sys->print("Minimum pegs: %d\n", minpegs);
+}
diff --git a/appl/math/doc.txt b/appl/math/doc.txt
new file mode 100644
index 00000000..821e3a22
--- /dev/null
+++ b/appl/math/doc.txt
@@ -0,0 +1,48 @@
+==========graph.b, gr.b================
+
+I believe scientific authors and readers are best served by a
+minimalist approach that makes all plots look boringly alike, except
+for the data content. Here is a library that does the task of
+accumulating data to find a bounding box and then draws the curves,
+points, and text surrounded by readable axes.
+
+The command appl/math/graph.dis is meant to be launched from wm and
+then provided a file containing (x,y) pairs. There are no options.
+
+The library version, appl/math/gr.b, is a little more flexible.
+
+Include gr.m and call p := GR->open(). To draw a curve, call
+ p.graph(x,y)
+where x and y are of real arrays. To add the string s at a point (u,v), call
+ p.text(j,s,u,v)
+where j is LJUST, CENTER, or RJUST to indicate whether the left, middle,
+or right of the string should be at (u,v) plus one of HIGH, MED, BASE, or LOW
+to indicate where the baseline of the text should be relative to (u,v).
+To get text running in the y direction, add UP. To enforce the
+same scaling in x and y, call p.equalxy().
+
+To change from the default solid line, call
+ p.pen(j) where j is one
+of the symbols: DASHED, DOTTED, REFERENCE, SOLID,
+CIRCLE, CROSS, or INVIS. A CIRCLE and CROSS
+"pen" just puts markers at the points, and doesn't connect with line
+segments. A REFERENCE line is lighter than other lines. DASHED
+curve follows the curve even within one dash, and preserves arclength of
+dashes and spaces. An INVIS line is sometimes handy as "strut" for
+maintaining a consistent scale across a series of plots.
+
+Finally, call
+ p.out(p,"xlabel","xunit","ylabel","yunit")
+to put out the curves and text in PostScript on standard output.
+Axes are produced with "Guggenheim slash notation," in which the user
+variable is divided by scaled units to get dimensionless numbers of
+reasonable magnitude.
+
+The function
+ name := p.getfilename();
+opens a dialog box to get a filename from the user and update the titlebar;
+ p.bye();
+pauses until the user clicks the "X" exit button in the titlebar.
+
+
+<ehg@bell-labs.com> 16 May 1996
diff --git a/appl/math/factor.b b/appl/math/factor.b
new file mode 100644
index 00000000..5f2aa775
--- /dev/null
+++ b/appl/math/factor.b
@@ -0,0 +1,180 @@
+#
+# initially generated by c2l
+#
+
+implement Factor;
+
+include "draw.m";
+
+Factor: module
+{
+ init: fn(nil: ref Draw->Context, argl: list of string);
+};
+
+include "sys.m";
+ sys: Sys;
+include "bufio.m";
+ bufio: Bufio;
+ Iobuf: import bufio;
+include "math.m";
+ maths: Math;
+ modf: import maths;
+
+init(nil: ref Draw->Context, argl: list of string)
+{
+ sys = load Sys Sys->PATH;
+ bufio = load Bufio Bufio->PATH;
+ maths = load Math Math->PATH;
+ main(len argl, argl);
+}
+
+WHLEN: con 48;
+wheel := array[WHLEN] of {
+ real 2,
+ real 10,
+ real 2,
+ real 4,
+ real 2,
+ real 4,
+ real 6,
+ real 2,
+ real 6,
+ real 4,
+ real 2,
+ real 4,
+ real 6,
+ real 6,
+ real 2,
+ real 6,
+ real 4,
+ real 2,
+ real 6,
+ real 4,
+ real 6,
+ real 8,
+ real 4,
+ real 2,
+ real 4,
+ real 2,
+ real 4,
+ real 8,
+ real 6,
+ real 4,
+ real 6,
+ real 2,
+ real 4,
+ real 6,
+ real 2,
+ real 6,
+ real 6,
+ real 4,
+ real 2,
+ real 4,
+ real 6,
+ real 2,
+ real 6,
+ real 4,
+ real 2,
+ real 4,
+ real 2,
+ real 10,
+};
+bin: ref Iobuf;
+
+main(argc: int, argv: list of string)
+{
+ n: real;
+ i: int;
+ l: string;
+
+ if(argc > 1){
+ argv = tl argv;
+ for(i = 1; i < argc; i++){
+ n = real hd argv;
+ factor(n);
+ argv = tl argv;
+ }
+ exit;
+ }
+ bin = bufio->fopen(sys->fildes(0), Sys->OREAD);
+ for(;;){
+ l = bin.gets('\n');
+ if(l == nil)
+ break;
+ n = real l;
+ if(n <= real 0)
+ break;
+ factor(n);
+ }
+ exit;
+}
+
+factor(n: real)
+{
+ quot, d, s: real;
+ i: int;
+
+ sys->print("%d\n", int n);
+ if(n == real 0)
+ return;
+ s = maths->sqrt(n)+real 1;
+ for(;;){
+ (iquot, frac) := modf(n/real 2);
+ if(frac != real 0)
+ break;
+ quot = real iquot;
+ sys->print(" 2\n");
+ n = quot;
+ s = maths->sqrt(n)+real 1;
+ }
+ for(;;){
+ (iquot, frac) := modf(n/real 3);
+ if(frac != real 0)
+ break;
+ quot = real iquot;
+ sys->print(" 3\n");
+ n = quot;
+ s = maths->sqrt(n)+real 1;
+ }
+ for(;;){
+ (iquot, frac) := modf(n/real 5);
+ if(frac != real 0)
+ break;
+ quot = real iquot;
+ sys->print(" 5\n");
+ n = quot;
+ s = maths->sqrt(n)+real 1;
+ }
+ for(;;){
+ (iquot, frac) := modf(n/real 7);
+ if(frac != real 0)
+ break;
+ quot = real iquot;
+ sys->print(" 7\n");
+ n = quot;
+ s = maths->sqrt(n)+real 1;
+ }
+ d = real 1;
+ for(i = 1;;){
+ d += wheel[i];
+ for(;;){
+ (iquot, frac) := modf(n/d);
+ if(frac != real 0)
+ break;
+ quot = real iquot;
+ sys->print(" %d\n", int d);
+ n = quot;
+ s = maths->sqrt(n)+real 1;
+ }
+ i++;
+ if(i >= WHLEN){
+ i = 0;
+ if(d > s)
+ break;
+ }
+ }
+ if(n > real 1)
+ sys->print(" %d\n", int n);
+ sys->print("\n");
+}
+
diff --git a/appl/math/ffts.b b/appl/math/ffts.b
new file mode 100644
index 00000000..e6132289
--- /dev/null
+++ b/appl/math/ffts.b
@@ -0,0 +1,639 @@
+implement FFTs;
+include "sys.m";
+ sys: Sys;
+ print: import sys;
+include "math.m";
+ math: Math;
+ cos, sin, Degree, Pi: import math;
+include "ffts.m";
+
+# by r. c. singleton, stanford research institute, sept. 1968
+# translated to limbo by eric grosse, jan 1997
+# arrays at(maxf), ck(maxf), bt(maxf), sk(maxf), and np(maxp)
+# are used for temporary storage. if the available storage
+# is insufficient, the program exits.
+# maxf must be >= the maximum prime factor of n.
+# maxp must be > the number of prime factors of n.
+# in addition, if the square-free portion k of n has two or
+# more prime factors, then maxp must be >= k-1.
+# array storage in nfac for a maximum of 15 prime factors of n.
+# if n has more than one square-free factor, the product of the
+# square-free factors must be <= 210
+
+ffts(a,b:array of real, ntot,n,nspan,isn:int){
+ maxp: con 209;
+ i,ii,inc,j,jc,jf,jj,k,k1,k2,k3,k4,kk:int;
+ ks,kspan,kspnn,kt,m,maxf,nn,nt:int;
+ aa,aj,ajm,ajp,ak,akm,akp,bb,bj,bjm,bjp,bk,bkm,bkp:real;
+ c1,c2,c3,c72,cd,rad,radf,s1,s2,s3,s72,s120,sd:real;
+ maxf = 23;
+ if(math == nil){
+ sys = load Sys Sys->PATH;
+ math = load Math Math->PATH;
+ }
+ nfac := array[12] of int;
+ np := array[maxp] of int;
+ at := array[23] of real;
+ ck := array[23] of real;
+ bt := array[23] of real;
+ sk := array[23] of real;
+
+ if(n<2) return;
+ inc = isn;
+ c72 = cos(72.*Degree);
+ s72 = sin(72.*Degree);
+ s120 = sin(120.*Degree);
+ rad = 2.*Pi;
+ if(isn<0){
+ s72 = -s72;
+ s120 = -s120;
+ rad = -rad;
+ inc = -inc;
+ }
+ nt = inc*ntot;
+ ks = inc*nspan;
+ kspan = ks;
+ nn = nt-inc;
+ jc = ks/n;
+ radf = rad*real(jc)*0.5;
+ i = 0;
+ jf = 0;
+
+ # determine the factors of n
+ m = 0;
+ k = n;
+ while(k==k/16*16){
+ m = m+1;
+ nfac[m] = 4;
+ k = k/16;
+ }
+ j = 3;
+ jj = 9;
+ for(;;)
+ if(k%jj==0){
+ m = m+1;
+ nfac[m] = j;
+ k = k/jj;
+ }else{
+ j = j+2;
+ jj = j*j;
+ if(jj>k)
+ break;
+ }
+ if(k<=4){
+ kt = m;
+ nfac[m+1] = k;
+ if(k!=1)
+ m = m+1;
+ }else{
+ if(k==k/4*4){
+ m = m+1;
+ nfac[m] = 2;
+ k = k/4;
+ }
+ kt = m;
+ j = 2;
+ do{
+ if(k%j==0){
+ m = m+1;
+ nfac[m] = j;
+ k = k/j;
+ }
+ j = ((j+1)/2)*2+1;
+ }while(j<=k);
+ }
+ if(kt!=0){
+ j = kt;
+ do{
+ m = m+1;
+ nfac[m] = nfac[j];
+ j = j-1;
+ }while(j!=0);
+ }
+
+ for(;;){ # compute fourier transform
+ sd = radf/real(kspan);
+ cd = sin(sd);
+ cd = 2.0*cd*cd;
+ sd = sin(sd+sd);
+ kk = 1;
+ i = i+1;
+ if(nfac[i]==2){ # transform for factor of 2 (including rotation factor)
+ kspan = kspan/2;
+ k1 = kspan+2;
+ for(;;){
+ k2 = kk+kspan;
+ ak = a[k2-1];
+ bk = b[k2-1];
+ a[k2-1] = a[kk-1]-ak;
+ b[k2-1] = b[kk-1]-bk;
+ a[kk-1] = a[kk-1]+ak;
+ b[kk-1] = b[kk-1]+bk;
+ kk = k2+kspan;
+ if(kk>nn){
+ kk = kk-nn;
+ if(kk>jc)
+ break;
+ }
+ }
+ if(kk>kspan)
+ break;
+ do{
+ c1 = 1.0-cd;
+ s1 = sd;
+ for(;;){
+ k2 = kk+kspan;
+ ak = a[kk-1]-a[k2-1];
+ bk = b[kk-1]-b[k2-1];
+ a[kk-1] = a[kk-1]+a[k2-1];
+ b[kk-1] = b[kk-1]+b[k2-1];
+ a[k2-1] = c1*ak-s1*bk;
+ b[k2-1] = s1*ak+c1*bk;
+ kk = k2+kspan;
+ if(kk>=nt){
+ k2 = kk-nt;
+ c1 = -c1;
+ kk = k1-k2;
+ if(kk<=k2){
+ ak = c1-(cd*c1+sd*s1);
+ s1 = (sd*c1-cd*s1)+s1;
+ c1 = 2.0-(ak*ak+s1*s1);
+ s1 = c1*s1;
+ c1 = c1*ak;
+ kk = kk+jc;
+ if(kk>=k2)
+ break;
+ }
+ }
+ }
+ k1 = k1+inc+inc;
+ kk = (k1-kspan)/2+jc;
+ }while(kk<=jc+jc);
+ }else{ # transform for factor of 4
+ if(nfac[i]!=4){
+ # transform for odd factors
+ k = nfac[i];
+ kspnn = kspan;
+ kspan = kspan/k;
+ if(k==3)
+ for(;;){
+ # transform for factor of 3 (optional code)
+ k1 = kk+kspan;
+ k2 = k1+kspan;
+ ak = a[kk-1];
+ bk = b[kk-1];
+ aj = a[k1-1]+a[k2-1];
+ bj = b[k1-1]+b[k2-1];
+ a[kk-1] = ak+aj;
+ b[kk-1] = bk+bj;
+ ak = -0.5*aj+ak;
+ bk = -0.5*bj+bk;
+ aj = (a[k1-1]-a[k2-1])*s120;
+ bj = (b[k1-1]-b[k2-1])*s120;
+ a[k1-1] = ak-bj;
+ b[k1-1] = bk+aj;
+ a[k2-1] = ak+bj;
+ b[k2-1] = bk-aj;
+ kk = k2+kspan;
+ if(kk>=nn){
+ kk = kk-nn;
+ if(kk>kspan)
+ break;
+ }
+ }
+ else if(k==5){
+ # transform for factor of 5 (optional code)
+ c2 = c72*c72-s72*s72;
+ s2 = 2.0*c72*s72;
+ for(;;){
+ k1 = kk+kspan;
+ k2 = k1+kspan;
+ k3 = k2+kspan;
+ k4 = k3+kspan;
+ akp = a[k1-1]+a[k4-1];
+ akm = a[k1-1]-a[k4-1];
+ bkp = b[k1-1]+b[k4-1];
+ bkm = b[k1-1]-b[k4-1];
+ ajp = a[k2-1]+a[k3-1];
+ ajm = a[k2-1]-a[k3-1];
+ bjp = b[k2-1]+b[k3-1];
+ bjm = b[k2-1]-b[k3-1];
+ aa = a[kk-1];
+ bb = b[kk-1];
+ a[kk-1] = aa+akp+ajp;
+ b[kk-1] = bb+bkp+bjp;
+ ak = akp*c72+ajp*c2+aa;
+ bk = bkp*c72+bjp*c2+bb;
+ aj = akm*s72+ajm*s2;
+ bj = bkm*s72+bjm*s2;
+ a[k1-1] = ak-bj;
+ a[k4-1] = ak+bj;
+ b[k1-1] = bk+aj;
+ b[k4-1] = bk-aj;
+ ak = akp*c2+ajp*c72+aa;
+ bk = bkp*c2+bjp*c72+bb;
+ aj = akm*s2-ajm*s72;
+ bj = bkm*s2-bjm*s72;
+ a[k2-1] = ak-bj;
+ a[k3-1] = ak+bj;
+ b[k2-1] = bk+aj;
+ b[k3-1] = bk-aj;
+ kk = k4+kspan;
+ if(kk>=nn){
+ kk = kk-nn;
+ if(kk>kspan)
+ break;
+ }
+ }
+ }else{
+ if(k!=jf){
+ jf = k;
+ s1 = rad/real(k);
+ c1 = cos(s1);
+ s1 = sin(s1);
+ if(jf>maxf){
+ sys->fprint(sys->fildes(2),"too many primes for fft");
+ exit;
+ }
+ ck[jf-1] = 1.0;
+ sk[jf-1] = 0.0;
+ j = 1;
+ do{
+ ck[j-1] = ck[k-1]*c1+sk[k-1]*s1;
+ sk[j-1] = ck[k-1]*s1-sk[k-1]*c1;
+ k = k-1;
+ ck[k-1] = ck[j-1];
+ sk[k-1] = -sk[j-1];
+ j = j+1;
+ }while(j<k);
+ }
+ for(;;){
+ k1 = kk;
+ k2 = kk+kspnn;
+ aa = a[kk-1];
+ bb = b[kk-1];
+ ak = aa;
+ bk = bb;
+ j = 1;
+ k1 = k1+kspan;
+ do{
+ k2 = k2-kspan;
+ j = j+1;
+ at[j-1] = a[k1-1]+a[k2-1];
+ ak = at[j-1]+ak;
+ bt[j-1] = b[k1-1]+b[k2-1];
+ bk = bt[j-1]+bk;
+ j = j+1;
+ at[j-1] = a[k1-1]-a[k2-1];
+ bt[j-1] = b[k1-1]-b[k2-1];
+ k1 = k1+kspan;
+ }while(k1<k2);
+ a[kk-1] = ak;
+ b[kk-1] = bk;
+ k1 = kk;
+ k2 = kk+kspnn;
+ j = 1;
+ do{
+ k1 = k1+kspan;
+ k2 = k2-kspan;
+ jj = j;
+ ak = aa;
+ bk = bb;
+ aj = 0.0;
+ bj = 0.0;
+ k = 1;
+ do{
+ k = k+1;
+ ak = at[k-1]*ck[jj-1]+ak;
+ bk = bt[k-1]*ck[jj-1]+bk;
+ k = k+1;
+ aj = at[k-1]*sk[jj-1]+aj;
+ bj = bt[k-1]*sk[jj-1]+bj;
+ jj = jj+j;
+ if(jj>jf)
+ jj = jj-jf;
+ }while(k<jf);
+ k = jf-j;
+ a[k1-1] = ak-bj;
+ b[k1-1] = bk+aj;
+ a[k2-1] = ak+bj;
+ b[k2-1] = bk-aj;
+ j = j+1;
+ }while(j<k);
+ kk = kk+kspnn;
+ if(kk>nn){
+ kk = kk-nn;
+ if(kk>kspan)
+ break;
+ }
+ }
+ }
+ # multiply by rotation factor (except for factors of 2 and 4)
+ if(i==m)
+ break;
+ kk = jc+1;
+ do{
+ c2 = 1.0-cd;
+ s1 = sd;
+ do{
+ c1 = c2;
+ s2 = s1;
+ kk = kk+kspan;
+ for(;;){
+ ak = a[kk-1];
+ a[kk-1] = c2*ak-s2*b[kk-1];
+ b[kk-1] = s2*ak+c2*b[kk-1];
+ kk = kk+kspnn;
+ if(kk>nt){
+ ak = s1*s2;
+ s2 = s1*c2+c1*s2;
+ c2 = c1*c2-ak;
+ kk = kk-nt+kspan;
+ if(kk>kspnn)
+ break;
+ }
+ }
+ c2 = c1-(cd*c1+sd*s1);
+ s1 = s1+(sd*c1-cd*s1);
+ c1 = 2.0-(c2*c2+s1*s1);
+ s1 = c1*s1;
+ c2 = c1*c2;
+ kk = kk-kspnn+jc;
+ }while(kk<=kspan);
+ kk = kk-kspan+jc+inc;
+ }while(kk<=jc+jc);
+ }else{
+ kspnn = kspan;
+ kspan = kspan/4;
+ do{
+ c1 = 1.;
+ s1 = 0.;
+ for(;;){
+ k1 = kk+kspan;
+ k2 = k1+kspan;
+ k3 = k2+kspan;
+ akp = a[kk-1]+a[k2-1];
+ akm = a[kk-1]-a[k2-1];
+ ajp = a[k1-1]+a[k3-1];
+ ajm = a[k1-1]-a[k3-1];
+ a[kk-1] = akp+ajp;
+ ajp = akp-ajp;
+ bkp = b[kk-1]+b[k2-1];
+ bkm = b[kk-1]-b[k2-1];
+ bjp = b[k1-1]+b[k3-1];
+ bjm = b[k1-1]-b[k3-1];
+ b[kk-1] = bkp+bjp;
+ bjp = bkp-bjp;
+ do10 := 0;
+ if(isn<0){
+ akp = akm+bjm;
+ akm = akm-bjm;
+ bkp = bkm-ajm;
+ bkm = bkm+ajm;
+ if(s1!=0.) do10 = 1;
+ }else{
+ akp = akm-bjm;
+ akm = akm+bjm;
+ bkp = bkm+ajm;
+ bkm = bkm-ajm;
+ if(s1!=0.) do10 = 1;
+ }
+ if(do10){
+ a[k1-1] = akp*c1-bkp*s1;
+ b[k1-1] = akp*s1+bkp*c1;
+ a[k2-1] = ajp*c2-bjp*s2;
+ b[k2-1] = ajp*s2+bjp*c2;
+ a[k3-1] = akm*c3-bkm*s3;
+ b[k3-1] = akm*s3+bkm*c3;
+ kk = k3+kspan;
+ if(kk<=nt)
+ continue;
+ }else{
+ a[k1-1] = akp;
+ b[k1-1] = bkp;
+ a[k2-1] = ajp;
+ b[k2-1] = bjp;
+ a[k3-1] = akm;
+ b[k3-1] = bkm;
+ kk = k3+kspan;
+ if(kk<=nt)
+ continue;
+ }
+ c2 = c1-(cd*c1+sd*s1);
+ s1 = (sd*c1-cd*s1)+s1;
+ c1 = 2.0-(c2*c2+s1*s1);
+ s1 = c1*s1;
+ c1 = c1*c2;
+ c2 = c1*c1-s1*s1;
+ s2 = 2.0*c1*s1;
+ c3 = c2*c1-s2*s1;
+ s3 = c2*s1+s2*c1;
+ kk = kk-nt+jc;
+ if(kk>kspan)
+ break;
+ }
+ kk = kk-kspan+inc;
+ }while(kk<=jc);
+ if(kspan==jc)
+ break;
+ }
+ }
+ } # end "compute fourier transform"
+
+ # permute the results to normal order---done in two stages
+ # permutation for square factors of n
+ np[0] = ks;
+ if(kt!=0){
+ k = kt+kt+1;
+ if(m<k)
+ k = k-1;
+ j = 1;
+ np[k] = jc;
+ do{
+ np[j] = np[j-1]/nfac[j];
+ np[k-1] = np[k]*nfac[j];
+ j = j+1;
+ k = k-1;
+ }while(j<k);
+ k3 = np[k];
+ kspan = np[1];
+ kk = jc+1;
+ k2 = kspan+1;
+ j = 1;
+ if(n!=ntot){
+ for(;;){
+ # permutation for multivariate transform
+ k = kk+jc;
+ do{
+ ak = a[kk-1];
+ a[kk-1] = a[k2-1];
+ a[k2-1] = ak;
+ bk = b[kk-1];
+ b[kk-1] = b[k2-1];
+ b[k2-1] = bk;
+ kk = kk+inc;
+ k2 = k2+inc;
+ }while(kk<k);
+ kk = kk+ks-jc;
+ k2 = k2+ks-jc;
+ if(kk>=nt){
+ k2 = k2-nt+kspan;
+ kk = kk-nt+jc;
+ if(k2>=ks)
+ permm: for(;;){
+ k2 = k2-np[j-1];
+ j = j+1;
+ k2 = np[j]+k2;
+ if(k2<=np[j-1]){
+ j = 1;
+ do{
+ if(kk<k2)
+ break permm;
+ kk = kk+jc;
+ k2 = kspan+k2;
+ }while(k2<ks);
+ if(kk>=ks)
+ break permm;
+ }
+ }
+ }
+ }
+ jc = k3;
+ }else{
+ for(;;){
+ # permutation for single-variate transform (optional code)
+ ak = a[kk-1];
+ a[kk-1] = a[k2-1];
+ a[k2-1] = ak;
+ bk = b[kk-1];
+ b[kk-1] = b[k2-1];
+ b[k2-1] = bk;
+ kk = kk+inc;
+ k2 = kspan+k2;
+ if(k2>=ks)
+ perms: for(;;){
+ k2 = k2-np[j-1];
+ j = j+1;
+ k2 = np[j]+k2;
+ if(k2<=np[j-1]){
+ j = 1;
+ do{
+ if(kk<k2)
+ break perms;
+ kk = kk+inc;
+ k2 = kspan+k2;
+ }while(k2<ks);
+ if(kk>=ks)
+ break perms;
+ }
+ }
+ }
+ jc = k3;
+ }
+ }
+ if(2*kt+1>=m)
+ return;
+ kspnn = np[kt];
+ # permutation for square-free factors of n
+ j = m-kt;
+ nfac[j+1] = 1;
+ do{
+ nfac[j] = nfac[j]*nfac[j+1];
+ j = j-1;
+ }while(j!=kt);
+ kt = kt+1;
+ nn = nfac[kt]-1;
+ if(nn<=maxp){
+ jj = 0;
+ j = 0;
+ for(;;){
+ k2 = nfac[kt];
+ k = kt+1;
+ kk = nfac[k];
+ j = j+1;
+ if(j>nn)
+ break;
+ for(;;){
+ jj = kk+jj;
+ if(jj<k2)
+ break;
+ jj = jj-k2;
+ k2 = kk;
+ k = k+1;
+ kk = nfac[k];
+ }
+ np[j-1] = jj;
+ }
+ # determine the permutation cycles of length greater than 1
+ j = 0;
+ for(;;){
+ j = j+1;
+ kk = np[j-1];
+ if(kk>=0)
+ if(kk==j){
+ np[j-1] = -j;
+ if(j==nn)
+ break;
+ }else{
+ do{
+ k = kk;
+ kk = np[k-1];
+ np[k-1] = -kk;
+ }while(kk!=j);
+ k3 = kk;
+ }
+ }
+ maxf = inc*maxf;
+ for(;;){
+ j = k3+1;
+ nt = nt-kspnn;
+ ii = nt-inc+1;
+ if(nt<0)
+ break;
+ for(;;){
+ j = j-1;
+ if(np[j-1]>=0){
+ jj = jc;
+ do{
+ kspan = jj;
+ if(jj>maxf)
+ kspan = maxf;
+ jj = jj-kspan;
+ k = np[j-1];
+ kk = jc*k+ii+jj;
+ k1 = kk+kspan;
+ k2 = 0;
+ do{
+ k2 = k2+1;
+ at[k2-1] = a[k1-1];
+ bt[k2-1] = b[k1-1];
+ k1 = k1-inc;
+ }while(k1!=kk);
+ do{
+ k1 = kk+kspan;
+ k2 = k1-jc*(k+np[k-1]);
+ k = -np[k-1];
+ do{
+ a[k1-1] = a[k2-1];
+ b[k1-1] = b[k2-1];
+ k1 = k1-inc;
+ k2 = k2-inc;
+ }while(k1!=kk);
+ kk = k2;
+ }while(k!=j);
+ k1 = kk+kspan;
+ k2 = 0;
+ do{
+ k2 = k2+1;
+ a[k1-1] = at[k2-1];
+ b[k1-1] = bt[k2-1];
+ k1 = k1-inc;
+ }while(k1!=kk);
+ }while(jj!=0);
+ if(j==1)
+ break;
+ }
+ }
+ }
+ }
+}
diff --git a/appl/math/fibonacci.b b/appl/math/fibonacci.b
new file mode 100644
index 00000000..af4e769c
--- /dev/null
+++ b/appl/math/fibonacci.b
@@ -0,0 +1,59 @@
+implement Fibonacci;
+
+include "sys.m";
+include "draw.m";
+
+Fibonacci: module
+{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+init(nil: ref Draw->Context, nil: list of string)
+{
+ sys := load Sys Sys->PATH;
+ for(i := 0; ; i++){
+ f := fibonacci(i);
+ if(f < 0)
+ break;
+ sys->print("F(%d) = %d\n", i, f);
+ }
+}
+
+FIB: exception(int, int);
+HELP: con "help";
+
+NOVAL: con -1000000000;
+
+fibonacci(n: int): int
+{
+ {
+ fib(1, n, 1, 1);
+ }
+ exception e{
+ FIB =>
+ (x, nil) := e;
+ return x;
+ * =>
+ return NOVAL;
+ }
+ return NOVAL;
+}
+
+fib(n: int, m: int, x: int, y: int) raises (FIB)
+{
+ if(n >= m)
+ raise FIB(x, y);
+
+ {
+ fib(n+1, m, x, y);
+ }
+ exception e{
+ FIB =>
+ (x, y) = e;
+ x = x+y;
+ y = x-y;
+ raise FIB(x, y);
+ * =>
+ raise HELP;
+ }
+}
diff --git a/appl/math/fit.b b/appl/math/fit.b
new file mode 100644
index 00000000..b512c13b
--- /dev/null
+++ b/appl/math/fit.b
@@ -0,0 +1,264 @@
+# fit a polynomial to a set of points
+# fit -dn [-v]
+# where n is the degree of the polynomial
+
+implement Fit;
+
+include "sys.m";
+ sys: Sys;
+include "draw.m";
+include "math.m";
+ maths: Math;
+include "bufio.m";
+ bufio: Bufio;
+include "arg.m";
+
+Fit: module
+{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+MAXPTS: con 512;
+MAXDEG: con 16;
+EPS: con 0.0000005;
+
+init(nil: ref Draw->Context, argv: list of string)
+{
+ sys = load Sys Sys->PATH;
+ maths = load Math Math->PATH;
+ if(maths == nil)
+ fatal(sys->sprint("cannot load maths library"));
+ bufio = load Bufio Bufio->PATH;
+ if(bufio == nil)
+ fatal(sys->sprint("cannot load bufio"));
+ main(argv);
+}
+
+isn(r: real, n: int): int
+{
+ s := r - real n;
+ if(s < 0.0)
+ s = -s;
+ return s < EPS;
+}
+
+fact(n: int): real
+{
+ f := 1.0;
+ for(i := 1; i <= n; i++)
+ f *= real i;
+ return f;
+}
+
+comb(n: int, r: int): real
+{
+ f := 1.0;
+ for(i := 0; i < r; i++)
+ f *= real (n-i);
+ return f/fact(r);
+}
+
+matalloc(n: int): array of array of real
+{
+ mat := array[n] of array of real;
+ for(i := 0; i < n; i++)
+ mat[i] = array[n] of real;
+ return mat;
+}
+
+matsalloc(n: int): array of array of array of real
+{
+ mats := array[n+1] of array of array of real;
+ for(i := 0; i <= n; i++)
+ mats[i] = matalloc(i);
+ return mats;
+}
+
+det(mat: array of array of real, n: int, mats: array of array of array of real): real
+{
+ # easy cases first
+ if(n == 0)
+ return 1.0;
+ if(n == 1)
+ return mat[0][0];
+ if(n == 2)
+ return mat[0][0]*mat[1][1]-mat[0][1]*mat[1][0];
+ d := 0.0;
+ s := 1;
+ m := mats[n-1];
+ for(k := 0; k < n; k++){
+ for(i := 0; i < n-1; i++){
+ for(j := 0; j < n-1; j++){
+ if(j < k)
+ m[i][j] = mat[i+1][j];
+ else
+ m[i][j] = mat[i+1][j+1];
+ }
+ }
+ d += (real s)*mat[0][k]*det(m, n-1, mats);
+ s = -s;
+ }
+ return d;
+}
+
+main(argv: list of string)
+{
+ i, j: int;
+ x, y, z: real;
+ fb: ref Bufio->Iobuf;
+
+ n := 0;
+ p := 1;
+ arg := load Arg Arg->PATH;
+ if(arg == nil)
+ fatal(sys->sprint("cannot load %s: %r", Arg->PATH));
+ arg->init(argv);
+ verbose := 0;
+ while((o := arg->opt()) != 0)
+ case o{
+ 'd' =>
+ p = int arg->arg();
+ 'v' =>
+ verbose = 1;
+ * =>
+ fatal(sys->sprint("bad option %c", o));
+ }
+ args := arg->argv();
+ arg = nil;
+ if(args != nil){
+ s := hd args;
+ fb = bufio->open(s, bufio->OREAD);
+ if(fb == nil)
+ fatal(sys->sprint("cannot open %s", s));
+ }
+ else{
+ fb = bufio->open("/dev/cons", bufio->OREAD);
+ if(fb == nil)
+ fatal(sys->sprint("missing data file name"));
+ }
+ a := array[p+1] of real;
+ b := array[p+1] of real;
+ sx := array[2*p+1] of real;
+ sxy := array[p+1] of real;
+ xd := array[MAXPTS] of real;
+ yd := array[MAXPTS] of real;
+ while(1){
+ xs := ss(bufio->fb.gett(" \t\r\n"));
+ if(xs == nil)
+ break;
+ ys := ss(bufio->fb.gett(" \t\r\n"));
+ if(ys == nil)
+ fatal(sys->sprint("missing value"));
+ if(n >= MAXPTS)
+ fatal(sys->sprint("too many points"));
+ xd[n] = real xs;
+ yd[n] = real ys;
+ n++;
+ }
+ if(p < 0)
+ fatal(sys->sprint("negative power"));
+ if(p > MAXDEG)
+ fatal(sys->sprint("power too large"));
+ if(n < p+1)
+ fatal(sys->sprint("not enough points"));
+ # use x-xbar, y-ybar to avoid overflow
+ for(i = 0; i <= p; i++)
+ sxy[i] = 0.0;
+ for(i = 0; i <= 2*p; i++)
+ sx[i] = 0.0;
+ xbar := ybar := 0.0;
+ for(i = 0; i < n; i++){
+ xbar += xd[i];
+ ybar += yd[i];
+ }
+ xbar = xbar/(real n);
+ ybar = ybar/(real n);
+ for(i = 0; i < n; i++){
+ x = xd[i]-xbar;
+ y = yd[i]-ybar;
+ for(j = 0; j <= p; j++)
+ sxy[j] += y*x**j;
+ for(j = 0; j <= 2*p; j++)
+ sx[j] += x**j;
+ }
+ mats := matsalloc(p+1);
+ mat := mats[p+1];
+ for(i = 0; i <= p; i++)
+ for(j = 0; j <= p; j++)
+ mat[i][j] = sx[i+j];
+ d := det(mat, p+1, mats);
+ if(isn(d, 0))
+ fatal(sys->sprint("points not independent"));
+ for(j = 0; j <= p; j++){
+ for(i = 0; i <= p; i++)
+ mat[i][j] = sxy[i];
+ a[j] = det(mat, p+1, mats)/d;
+ for(i = 0; i <= p; i++)
+ mat[i][j] = sx[i+j];
+ }
+ if(verbose)
+ sys->print("\npt actual x actual y predicted y\n");
+ e := 0.0;
+ for(i = 0; i < n; i++){
+ x = xd[i]-xbar;
+ y = yd[i]-ybar;
+ z = 0.0;
+ for(j = 0; j <= p; j++)
+ z += a[j]*x**j;
+ z += ybar;
+ e += (z-yd[i])*(z-yd[i]);
+ if(verbose)
+ sys->print("%d. %f %f %f\n", i+1, xd[i], yd[i], z);
+ }
+ if(verbose)
+ sys->print("root mean squared error = %f\n", maths->sqrt(e/(real n)));
+ for(i = 0; i <= p; i++)
+ b[i] = 0.0;
+ b[0] += ybar;
+ for(i = 0; i <= p; i++)
+ for(j = 0; j <= i; j++)
+ b[j] += a[i]*comb(i, j)*(-xbar)**(i-j);
+ pr := 0;
+ sys->print("y = ");
+ for(i = p; i >= 0; i--){
+ if(!isn(b[i], 0) || (i == 0 && pr == 0)){
+ if(b[i] < 0.0){
+ sys->print("-");
+ b[i] = -b[i];
+ }
+ else if(pr)
+ sys->print("+");
+ pr = 1;
+ if(i == 0)
+ sys->print("%f", b[i]);
+ else{
+ if(!isn(b[i], 1))
+ sys->print("%f*", b[i]);
+ sys->print("x");
+ if(i > 1)
+ sys->print("^%d", i);
+ }
+ }
+ }
+ sys->print("\n");
+}
+
+ss(s: string): string
+{
+ l := len s;
+ while(l > 0 && (s[0] == ' ' || s[0] == '\t' || s[0] == '\r' || s[0] == '\n')){
+ s = s[1: ];
+ l--;
+ }
+ while(l > 0 && (s[l-1] == ' ' || s[l-1] == '\t' || s[l-1] == '\r' || s[l-1] == '\n')){
+ s = s[0: l-1];
+ l--;
+ }
+ return s;
+}
+
+fatal(s: string)
+{
+ sys->fprint(sys->fildes(2), "fit: %s\n", s);
+ exit;
+}
diff --git a/appl/math/genprimes.b b/appl/math/genprimes.b
new file mode 100644
index 00000000..19c32938
--- /dev/null
+++ b/appl/math/genprimes.b
@@ -0,0 +1,64 @@
+implement Primes;
+
+include "draw.m";
+
+Primes: module
+{
+ init: fn(nil: ref Draw->Context, argl: list of string);
+};
+
+include "sys.m";
+ sys: Sys;
+include "arg.m";
+ arg: Arg;
+
+LIM: con 1729;
+MAX: con 1000000;
+BUFSZ: con 256;
+
+init(nil: ref Draw->Context, argl: list of string)
+{
+ sys = load Sys Sys->PATH;
+ arg = load Arg Arg->PATH;
+ arg->init(argl);
+ quiet := 0;
+ lim := LIM;
+ while((ch := arg->opt()) != 0){
+ case (ch){
+ 'q' =>
+ quiet = 1;
+ * =>
+ ;
+ }
+ }
+ argv := arg->argv();
+ if(argv != nil)
+ lim = int hd argv;
+ if(lim < 2)
+ lim = 2;
+ if(lim > MAX)
+ lim = MAX;
+ c := chan[BUFSZ] of int;
+ spawn prime(c, !quiet);
+ for(n := 2; n <= lim; n++)
+ c <-= n;
+ c <-= 1;
+}
+
+prime(c: chan of int, pr: int)
+{
+ p := <-c;
+ if(p == 1)
+ exit;
+ if(pr)
+ sys->print("%d\n", p);
+ nc := chan[BUFSZ] of int;
+ spawn prime(nc, pr);
+ for(;;){
+ n := <-c;
+ if(n%p)
+ nc <-= n;
+ if(n == 1)
+ exit;
+ }
+}
diff --git a/appl/math/geodesy.b b/appl/math/geodesy.b
new file mode 100644
index 00000000..7c4cd152
--- /dev/null
+++ b/appl/math/geodesy.b
@@ -0,0 +1,849 @@
+implement Geodesy;
+
+include "sys.m";
+ sys: Sys;
+include "math.m";
+ maths: Math;
+ Pi: import Math;
+ sin, cos, tan, asin, acos, atan, atan2, sqrt, fabs: import maths;
+include "math/geodesy.m";
+
+Approx: con 0;
+
+Epsilon: con 0.000001;
+Mperft: con 0.3048;
+Earthrad: con 10800.0/Pi*6076.115*Mperft; # in feet (about 4000 miles) : now metres
+Δt: con 16.0; # now-1989
+
+# lalo0: con "53:57:45N 01:04:55W";
+# os0: con "SE6022552235";
+
+# ellipsoids
+Airy1830, Airy1830m, Int1924, GRS80: con iota;
+
+Ngrid: con 100000; # in metres
+
+Vector: adt{
+ x, y, z: real;
+};
+
+Latlong: adt{
+ la: real; # -Pi to Pi
+ lo: real; # -Pi to Pi
+ x: real;
+ y: real;
+};
+
+Ellipsoid: adt{
+ name: string;
+ a: real;
+ b: real;
+};
+
+Datum: adt{
+ name: string;
+ e: int;
+ # X, Y, Z axes etc
+};
+
+Mercator: adt{
+ name: string;
+ F0: real;
+ φ0λ0: string;
+ E0: real;
+ N0: real;
+ e: int;
+};
+
+Helmert: adt{
+ tx, ty, tz: real; # metres
+ s: real; # ppm
+ rx, ry, rz: real; # secs
+};
+
+Format: adt{
+ dat: int; # datum
+ cdat: int; # converting datum
+ prj: int; # projection
+ tmp: ref Mercator; # actual projection
+ orig: Lalo; # origin of above projection
+ zone: int; # UTM zone
+};
+
+# ellipsoids
+ells := array[] of {
+ Airy1830 => Ellipsoid("Airy1830", 6377563.396, 6356256.910),
+ Airy1830m => Ellipsoid("Airy1830 modified", 6377340.189, 6356034.447),
+ Int1924 => Ellipsoid("International 1924", 6378388.000, 6356911.946),
+ GRS80 => Ellipsoid("GRS80", 6378137.000, 6356752.3141),
+ };
+
+# datums
+dats := array[] of {
+ OSGB36 => Datum("OSGB36", Airy1830),
+ Ireland65 => Datum("Ireland65", Airy1830m),
+ ED50 => Datum("ED50", Int1924),
+ WGS84 => Datum("WGS84", GRS80),
+ ITRS2000 => Datum("ITRS2000", GRS80),
+ ETRS89 => Datum("ETRS89", GRS80),
+ };
+
+# transverse Mercator projections
+tmps := array[] of {
+ Natgrid => Mercator("National Grid", 0.9996012717, "49:00:00N 02:00:00W", real(4*Ngrid), real(-Ngrid), Airy1830),
+ IrishNatgrid => Mercator("Irish National Grid", 1.000035, "53:30:00N 08:00:00W", real(2*Ngrid), real(5*Ngrid/2), Airy1830m),
+ UTMEur => Mercator("UTM Europe", 0.9996, nil, real(5*Ngrid), real(0), Int1924),
+ UTM => Mercator("UTM", 0.9996, nil, real(5*Ngrid), real(0), GRS80),
+ };
+
+# Helmert tranformations
+HT_WGS84_OSGB36: con Helmert(-446.448, 125.157, -542.060, 20.4894, -0.1502, -0.2470, -0.8421);
+HT_ITRS2000_ETRS89: con Helmert(0.054, 0.051, -0.048, 0.0, 0.000081*Δt, 0.00049*Δt, -0.000792*Δt);
+
+# Helmert matrices
+HM_WGS84_OSGB36, HM_OSGB36_WGS84, HM_ITRS2000_ETRS89, HM_ETRS89_ITRS2000, HM_ETRS89_OSGB36, HM_OSGB36_ETRS89, HM_IDENTITY: array of array of real;
+
+fmt: ref Format;
+
+# latlong: ref Latlong;
+
+init(d: int, t: int, z: int)
+{
+ sys = load Sys Sys->PATH;
+ maths = load Math Math->PATH;
+
+ helmertinit();
+ format(d, t, z);
+ # (nil, (la, lo)) := str2lalo(lalo0);
+ # (nil, (E, N)) := os2en(os0);
+ # latlong = ref Latlong(la, lo, real E, real N);
+}
+
+format(d: int, t: int, z: int)
+{
+ if(fmt == nil)
+ fmt = ref Format(WGS84, 0, Natgrid, nil, (0.0, 0.0), 30);
+ if(d >= 0 && d <= ETRS89)
+ fmt.dat = d;
+ if(t >= 0 && t <= UTM)
+ fmt.prj = t;
+ if(z >= 1 && z <= 60)
+ fmt.zone = z;
+ fmt.cdat = fmt.dat;
+ fmt.tmp = ref Mercator(tmps[fmt.prj]);
+ if(fmt.tmp.φ0λ0 == nil)
+ fmt.orig = utmlaloz(fmt.zone);
+ else
+ (nil, fmt.orig) = str2lalo(fmt.tmp.φ0λ0);
+ e := fmt.tmp.e;
+ if(e != dats[fmt.dat].e){
+ for(i := 0; i <= ETRS89; i++)
+ if(e == dats[i].e){
+ fmt.cdat = i;
+ break;
+ }
+ }
+}
+
+str2en(s: string): (int, Eano)
+{
+ s = trim(s, " \t\n\r");
+ if(s == nil)
+ return (0, (0.0, 0.0));
+ os := s[0] >= 'A' && s[0] <= 'Z' || strchrs(s, "NSEW:") < 0;
+ en: Eano;
+ if(os){
+ (ok, p) := os2en(s);
+ if(!ok)
+ return (0, (0.0, 0.0));
+ en = p;
+ }
+ else{
+ (ok, lalo) := str2lalo(s);
+ if(!ok)
+ return (0, (0.0, 0.0));
+ en = lalo2en(lalo);
+ }
+ return (1, en);
+}
+
+str2ll(s: string, pos: int, neg: int): (int, real)
+{
+ (n, ls) := sys->tokenize(s, ": \t");
+ if(n < 1 || n > 3)
+ return (0, 0.0);
+ t := hd ls; ls = tl ls;
+ v := real t;
+ if(ls != nil){
+ t = hd ls; ls = tl ls;
+ v += (real t)/60.0;
+ }
+ if(ls != nil){
+ t = hd ls; ls = tl ls;
+ v += (real t)/3600.0;
+ }
+ c := t[len t-1];
+ if(c == pos)
+ ;
+ else if(c == neg)
+ v = -v;
+ else
+ return (0, 0.0);
+ return (1, norm(deg2rad(v)));
+}
+
+str2lalo(s: string): (int, Lalo)
+{
+ s = trim(s, " \t\n\r");
+ p := strchr(s, 'N');
+ if(p < 0)
+ p = strchr(s, 'S');
+ if(p < 0)
+ return (0, (0.0, 0.0));
+ (ok1, la) := str2ll(s[0: p+1], 'N', 'S');
+ (ok2, lo) := str2ll(s[p+1: ], 'E', 'W');
+ if(!ok1 || !ok2 || la < -Pi/2.0 || la > Pi/2.0)
+ return (0, (0.0, 0.0));
+ return (1, (la, lo));
+}
+
+ll2str(ll: int, dir: string): string
+{
+ d := ll/360000;
+ ll -= 360000*d;
+ m := ll/6000;
+ ll -= 6000*m;
+ s := ll/100;
+ ll -= 100*s;
+ return d2(d) + ":" + d2(m) + ":" + d2(s) + "." + d2(ll) + dir;
+}
+
+lalo2str(lalo: Lalo): string
+{
+ la := int(360000.0*rad2deg(lalo.la));
+ lo := int(360000.0*rad2deg(lalo.lo));
+ lad := "N";
+ lod := "E";
+ if(la < 0){
+ lad = "S";
+ la = -la;
+ }
+ if(lo < 0){
+ lod = "W";
+ lo = -lo;
+ }
+ return ll2str(la, lad) + " " + ll2str(lo, lod);
+}
+
+en2os(p: Eano): string
+{
+ E := trunc(p.e);
+ N := trunc(p.n);
+ es := E/Ngrid;
+ ns := N/Ngrid;
+ e := E-Ngrid*es;
+ n := N-Ngrid*ns;
+ d1 := 5*(4-ns/5)+es/5+'A'-3;
+ d2 := 5*(4-ns%5)+es%5+'A';
+ # now account for 'I' missing
+ if(d1 >= 'I')
+ d1++;
+ if(d2 >= 'I')
+ d2++;
+ return sys->sprint("%c%c%5.5d%5.5d", d1, d2, e, n);
+}
+
+os2en(s: string): (int, Eano)
+{
+ s = trim(s, " \t\n\r");
+ if((m := len s) != 4 && m != 6 && m != 8 && m != 10 && m != 12)
+ return (0, (0.0, 0.0));
+ m = m/2-1;
+ u := Ngrid/10**m;
+ d1 := s[0];
+ d2 := s[1];
+ if(d1 < 'A' || d2 < 'A' || d1 > 'Z' || d2 > 'Z'){
+ # error(sys->sprint("bad os reference %s", s));
+ e := u*int s[0: 1+m];
+ n := u*int s[1+m: 2+2*m];
+ return (1, (real e, real n));
+ }
+ e := u*int s[2: 2+m];
+ n := u*int s[2+m: 2+2*m];
+ if(d1 >= 'I')
+ d1--;
+ if(d2 >= 'I')
+ d2--;
+ d1 -= 'A'-3;
+ d2 -= 'A';
+ es := 5*(d1%5)+d2%5;
+ ns := 5*(4-d1/5)+4-d2/5;
+ return (1, (real(Ngrid*es+e), real(Ngrid*ns+n)));
+}
+
+utmlalo(lalo: Lalo): Lalo
+{
+ (nil, zn) := utmzone(lalo);
+ return utmlaloz(zn);
+}
+
+utmlaloz(zn: int): Lalo
+{
+ return (0.0, deg2rad(real(6*zn-183)));
+}
+
+utmzone(lalo: Lalo): (int, int)
+{
+ (la, lo) := lalo;
+ la = rad2deg(la);
+ lo = rad2deg(lo);
+ zlo := trunc(lo+180.0)/6+1;
+ if(la < -80.0)
+ zla := 'B';
+ else if(la >= 84.0)
+ zla = 'Y';
+ else if(la >= 72.0)
+ zla = 'X';
+ else{
+ zla = trunc(la+80.0)/8+'C';
+ if(zla >= 'I')
+ zla++;
+ if(zla >= 'O')
+ zla++;
+ }
+ return (zla, zlo);
+}
+
+helmertinit()
+{
+ (HM_WGS84_OSGB36, HM_OSGB36_WGS84) = helminit(HT_WGS84_OSGB36);
+ (HM_ITRS2000_ETRS89, HM_ETRS89_ITRS2000) = helminit(HT_ITRS2000_ETRS89);
+ HM_ETRS89_OSGB36 = mulmm(HM_WGS84_OSGB36, HM_ETRS89_ITRS2000);
+ HM_OSGB36_ETRS89 = mulmm(HM_ITRS2000_ETRS89, HM_OSGB36_WGS84);
+ HM_IDENTITY = m := matrix(3, 4);
+ m[0][0] = m[1][1] = m[2][2] = 1.0;
+ # mprint(HM_WGS84_OSGB36);
+ # mprint(HM_OSGB36_WGS84);
+}
+
+helminit(h: Helmert): (array of array of real, array of array of real)
+{
+ m := matrix(3, 4);
+
+ s := 1.0+h.s/1000000.0;
+ rx := sec2rad(h.rx);
+ ry := sec2rad(h.ry);
+ rz := sec2rad(h.rz);
+
+ m[0][0] = s;
+ m[0][1] = -rz;
+ m[0][2] = ry;
+ m[0][3] = h.tx;
+ m[1][0] = rz;
+ m[1][1] = s;
+ m[1][2] = -rx;
+ m[1][3] = h.ty;
+ m[2][0] = -ry;
+ m[2][1] = rx;
+ m[2][2] = s;
+ m[2][3] = h.tz;
+
+ return (m, inv(m));
+}
+
+trans(f: int, t: int): array of array of real
+{
+ case(f){
+ WGS84 =>
+ case(t){
+ WGS84 =>
+ return HM_IDENTITY;
+ OSGB36 =>
+ return HM_WGS84_OSGB36;
+ ITRS2000 =>
+ return HM_IDENTITY;
+ ETRS89 =>
+ return HM_ITRS2000_ETRS89;
+ }
+ OSGB36 =>
+ case(t){
+ WGS84 =>
+ return HM_OSGB36_WGS84;
+ OSGB36 =>
+ return HM_IDENTITY;
+ ITRS2000 =>
+ return HM_OSGB36_WGS84;
+ ETRS89 =>
+ return HM_OSGB36_ETRS89;
+ }
+ ITRS2000 =>
+ case(t){
+ WGS84 =>
+ return HM_IDENTITY;
+ OSGB36 =>
+ return HM_WGS84_OSGB36;
+ ITRS2000 =>
+ return HM_IDENTITY;
+ ETRS89 =>
+ return HM_ITRS2000_ETRS89;
+ }
+ ETRS89 =>
+ case(t){
+ WGS84 =>
+ return HM_ETRS89_ITRS2000;
+ OSGB36 =>
+ return HM_ETRS89_OSGB36;
+ ITRS2000 =>
+ return HM_ETRS89_ITRS2000;
+ ETRS89 =>
+ return HM_IDENTITY;
+ }
+ }
+ return HM_IDENTITY; # Ireland65, ED50 not done
+}
+
+datum2datum(lalo: Lalo, f: int, t: int): Lalo
+{
+ if(f == t)
+ return lalo;
+ (la, lo) := lalo;
+ v := laloh2xyz(la, lo, 0.0, dats[f].e);
+ v = mulmv(trans(f, t), v);
+ (la, lo, nil) = xyz2laloh(v, dats[t].e);
+ return (la, lo);
+}
+
+laloh2xyz(φ: real, λ: real, H: real, e: int): Vector
+{
+ a := ells[e].a;
+ b := ells[e].b;
+ e2 := 1.0-(b/a)**2;
+
+ s := sin(φ);
+ c := cos(φ);
+
+ ν := a/sqrt(1.0-e2*s*s);
+ x := (ν+H)*c*cos(λ);
+ y := (ν+H)*c*sin(λ);
+ z := ((1.0-e2)*ν+H)*s;
+
+ return (x, y, z);
+}
+
+xyz2laloh(v: Vector, e: int): (real, real, real)
+{
+ x := v.x;
+ y := v.y;
+ z := v.z;
+
+ a := ells[e].a;
+ b := ells[e].b;
+ e2 := 1.0-(b/a)**2;
+
+ λ := atan2(y, x);
+
+ p := sqrt(x*x+y*y);
+ φ := φ1 := atan(z/(p*(1.0-e2)));
+ ν := 0.0;
+ do{
+ φ = φ1;
+ s := sin(φ);
+ ν = a/sqrt(1.0-e2*s*s);
+ φ1 = atan((z+e2*ν*s)/p);
+ }while(!small(fabs(φ-φ1)));
+
+ φ = φ1;
+ H := p/cos(φ)-ν;
+
+ return (φ, λ, H);
+}
+
+lalo2en(lalo: Lalo): Eano
+{
+ (φ, λ) := lalo;
+ if(fmt.cdat != fmt.dat)
+ (φ, λ) = datum2datum(lalo, fmt.dat, fmt.cdat);
+
+ s := sin(φ);
+ c := cos(φ);
+ t2 := tan(φ)**2;
+
+ (nil, F0, φ0λ0, E0, N0, e) := *fmt.tmp;
+ a := ells[e].a;
+ b := ells[e].b;
+ e2 := 1.0-(b/a)**2;
+
+ if(φ0λ0 == nil) # UTM
+ (φ0, λ0) := utmlalo((φ, λ)); # don't use fmt.zone here
+ else
+ (φ0, λ0) = fmt.orig;
+
+ n := (a-b)/(a+b);
+ ν := a*F0/sqrt(1.0-e2*s*s);
+ ρ := ν*(1.0-e2)/(1.0-e2*s*s);
+ η2 := ν/ρ-1.0;
+
+ φ1 := φ-φ0;
+ φ2 := φ+φ0;
+ M := b*F0*((1.0+n*(1.0+1.25*n*(1.0+n)))*φ1 - (3.0*n*(1.0+n*(1.0+0.875*n)))*sin(φ1)*cos(φ2) + 1.875*n*n*(1.0+n)*sin(2.0*φ1)*cos(2.0*φ2) - 35.0/24.0*n**3*sin(3.0*φ1)*cos(3.0*φ2));
+
+ I := M+N0;
+ II := ν*s*c/2.0;
+ III := ν*s*c**3*(5.0-t2+9.0*η2)/24.0;
+ IIIA := ν*s*c**5*(61.0+t2*(t2-58.0))/720.0;
+ IV := ν*c;
+ V := ν*c**3*(ν/ρ-t2)/6.0;
+ VI := ν*c**5*(5.0+14.0*η2+t2*(t2-18.0-58.0*η2))/120.0;
+
+ λ -= λ0;
+ λ2 := λ*λ;
+ N := I+λ2*(II+λ2*(III+IIIA*λ2));
+ E := E0+λ*(IV+λ2*(V+VI*λ2));
+
+ # if(E < 0.0 || E >= real(7*Ngrid))
+ # E = 0.0;
+ # if(N < 0.0 || N >= real(13*Ngrid))
+ # N = 0.0;
+ return (E, N);
+}
+
+en2lalo(en: Eano): Lalo
+{
+ E := en.e;
+ N := en.n;
+
+ (nil, F0, nil, E0, N0, e) := *fmt.tmp;
+ a := ells[e].a;
+ b := ells[e].b;
+ e2 := 1.0-(b/a)**2;
+
+ (φ0, λ0) := fmt.orig;
+
+ n := (a-b)/(a+b);
+
+ M0 := 1.0+n*(1.0+1.25*n*(1.0+n));
+ M1 := 3.0*n*(1.0+n*(1.0+0.875*n));
+ M2 := 1.875*n*n*(1.0+n);
+ M3 := 35.0/24.0*n**3;
+
+ N -= N0;
+ M := 0.0;
+ φ := φold := φ0;
+ do{
+ φ = (N-M)/(a*F0)+φold;
+ φ1 := φ-φ0;
+ φ2 := φ+φ0;
+ M = b*F0*(M0*φ1 - M1*sin(φ1)*cos(φ2) + M2*sin(2.0*φ1)*cos(2.0*φ2) - M3*sin(3.0*φ1)*cos(3.0*φ2));
+ φold = φ;
+ }while(fabs(N-M) >= 0.01);
+
+ s := sin(φ);
+ c := cos(φ);
+ t := tan(φ);
+ t2 := t*t;
+
+ ν := a*F0/sqrt(1.0-e2*s*s);
+ ρ := ν*(1.0-e2)/(1.0-e2*s*s);
+ η2 := ν/ρ-1.0;
+
+ VII := t/(2.0*ρ*ν);
+ VIII := VII*(5.0+η2+3.0*t2*(1.0-3.0*η2))/(12.0*ν*ν);
+ IX := VII*(61.0+45.0*t2*(2.0+t2))/(360.0*ν**4);
+ X := 1.0/(ν*c);
+ XI := X*(ν/ρ+2.0*t2)/(6.0*ν*ν);
+ XII := X*(5.0+4.0*t2*(7.0+6.0*t2))/(120.0*ν**4);
+ XIIA := X*(61.0+2.0*t2*(331.0+60.0*t2*(11.0+6.0*t2)))/(5040.0*ν**6);
+
+ E -= E0;
+ E2 := E*E;
+ φ = φ-E2*(VII-E2*(VIII-E2*IX));
+ λ := λ0+E*(X-E2*(XI-E2*(XII-E2*XIIA)));
+
+ if(fmt.cdat != fmt.dat)
+ (φ, λ) = datum2datum((φ, λ), fmt.cdat, fmt.dat);
+ return (φ, λ);
+}
+
+mulmm(m1: array of array of real, m2: array of array of real): array of array of real
+{
+ m := matrix(3, 4);
+ mul3x3(m, m1, m2);
+ for(i := 0; i < 3; i++){
+ sum := 0.0;
+ for(k := 0; k < 3; k++)
+ sum += m1[i][k]*m2[k][3];
+ m[i][3] = sum+m1[i][3];
+ }
+ return m;
+}
+
+mulmv(m: array of array of real, v: Vector): Vector
+{
+ x := v.x;
+ y := v.y;
+ z := v.z;
+ v.x = m[0][0]*x + m[0][1]*y + m[0][2]*z + m[0][3];
+ v.y = m[1][0]*x + m[1][1]*y + m[1][2]*z + m[1][3];
+ v.z = m[2][0]*x + m[2][1]*y + m[2][2]*z + m[2][3];
+ return v;
+}
+
+inv(m: array of array of real): array of array of real
+{
+ n := matrix(3, 4);
+ inv3x3(m, n);
+ (n[0][3], n[1][3], n[2][3]) = mulmv(n, (-m[0][3], -m[1][3], -m[2][3]));
+ return n;
+}
+
+mul3x3(m: array of array of real, m1: array of array of real, m2: array of array of real)
+{
+ for(i := 0; i < 3; i++){
+ for(j := 0; j < 3; j++){
+ sum := 0.0;
+ for(k := 0; k < 3; k++)
+ sum += m1[i][k]*m2[k][j];
+ m[i][j] = sum;
+ }
+ }
+}
+
+inv3x3(m: array of array of real, n: array of array of real)
+{
+ t00 := m[0][0];
+ t01 := m[0][1];
+ t02 := m[0][2];
+ t10 := m[1][0];
+ t11 := m[1][1];
+ t12 := m[1][2];
+ t20 := m[2][0];
+ t21 := m[2][1];
+ t22 := m[2][2];
+
+ n[0][0] = t11*t22-t12*t21;
+ n[1][0] = t12*t20-t10*t22;
+ n[2][0] = t10*t21-t11*t20;
+ n[0][1] = t02*t21-t01*t22;
+ n[1][1] = t00*t22-t02*t20;
+ n[2][1] = t01*t20-t00*t21;
+ n[0][2] = t01*t12-t02*t11;
+ n[1][2] = t02*t10-t00*t12;
+ n[2][2] = t00*t11-t01*t10;
+
+ d := t00*n[0][0]+t01*n[1][0]+t02*n[2][0];
+ for(i := 0; i < 3; i++)
+ for(j := 0; j < 3; j++)
+ n[i][j] /= d;
+}
+
+matrix(rows: int, cols: int): array of array of real
+{
+ m := array[rows] of array of real;
+ for(i := 0; i < rows; i++)
+ m[i] = array[cols] of { * => 0.0 };
+ return m;
+}
+
+vprint(v: Vector)
+{
+ sys->print(" %f %f %f\n", v.x, v.y, v.z);
+}
+
+mprint(m: array of array of real)
+{
+ for(i := 0; i < len m; i++){
+ for(j := 0; j < len m[i]; j++)
+ sys->print(" %f", m[i][j]);
+ sys->print("\n");
+ }
+}
+
+# lalo2xy(la: real, lo: real, lalo: ref Latlong): Eano
+# {
+# x, y: real;
+#
+# la0 := lalo.la;
+# lo0 := lalo.lo;
+# if(Approx){
+# x = Earthrad*cos(la0)*(lo-lo0)+lalo.x;
+# y = Earthrad*(la-la0)+lalo.y;
+# }
+# else{
+# x = Earthrad*cos(la)*sin(lo-lo0)+lalo.x;
+# y = Earthrad*(sin(la)*cos(la0)-sin(la0)*cos(la)*cos(lo-lo0))+lalo.y;
+# }
+# return (x, y);
+# }
+
+# lalo2xyz(la: real, lo: real, lalo: ref Latlong): (int, int, int)
+# {
+# z: real;
+#
+# la0 := lalo.la;
+# lo0 := lalo.lo;
+# (x, y) := lalo2xy(la, lo, lalo);
+# if(Approx)
+# z = Earthrad;
+# else
+# z = Earthrad*(sin(la)*sin(la0)+cos(la)*cos(la0)*cos(lo-lo0));
+# return (x, y, int z);
+# }
+
+# xy2lalo(p: Eano, lalo: ref Latlong): (real, real)
+# {
+# la, lo: real;
+#
+# x := p.e;
+# y := p.n;
+# la0 := lalo.la;
+# lo0 := lalo.lo;
+# if(Approx){
+# la = la0 + (y-lalo.y)/Earthrad;
+# lo = lo0 + (x-lalo.x)/(Earthrad*cos(la0));
+# }
+# else{
+# a, b, c, d, bestd, r, r1, r2, lat, lon, tmp: real;
+# i, n: int;
+#
+# bestd = -1.0;
+# la = lo = 0.0;
+# a = (x-lalo.x)/Earthrad;
+# b = (y-lalo.y)/Earthrad;
+# (n, r1, r2) = quad(1.0, -2.0*b*cos(la0), (a*a-1.0)*sin(la0)*sin(la0)+b*b);
+# if(n == 0)
+# return (la, lo);
+# while(--n >= 0){
+# if(n == 1)
+# r = r2;
+# else
+# r = r1;
+# if(fabs(r) <= 1.0){
+# lat = asin(r);
+# c = cos(lat);
+# if(small(c))
+# tmp = 0.0; # lat = +90, -90, lon = lo0
+# else
+# tmp = a/c;
+# if(fabs(tmp) <= 1.0){
+# for(i = 0; i < 2; i++){
+# if(i == 0)
+# lon = norm(asin(tmp)+lo0);
+# else
+# lon = norm(Pi-asin(tmp)+lo0);
+# (X, Y, Z) := lalo2xyz(lat, lon, lalo);
+# # eliminate non-roots by d, root on other side of earth by Z
+# d = (real X-x)**2+(real Y-y)**2;
+# if(Z >= 0 && (bestd < 0.0 || d < bestd)){
+# bestd = d;
+# la = lat;
+# lo = lon;
+# }
+# }
+# }
+# }
+# }
+# }
+# return (la, lo);
+# }
+
+# quad(a: real, b: real, c: real): (int, real, real)
+# {
+# r1, r2: real;
+#
+# D := b*b-4.0*a*c;
+# if(small(a)){
+# if(small(b))
+# return (0, r1, r2);
+# r1 = r2 = -c/b;
+# return (1, r1, r2);
+# }
+# if(D < 0.0)
+# return (0, r1, r2);
+# D = sqrt(D);
+# r1 = (-b+D)/(2.0*a);
+# r2 = (-b-D)/(2.0*a);
+# if(small(D))
+# return (1, r1, r2);
+# else
+# return (2, r1, r2);
+# }
+
+d2(v: int): string
+{
+ s := string v;
+ if(v < 10)
+ s = "0" + s;
+ return s;
+}
+
+trim(s: string, t: string): string
+{
+ while(s != nil && strchr(t, s[0]) >= 0)
+ s = s[1: ];
+ while(s != nil && strchr(t, s[len s-1]) >= 0)
+ s = s[0: len s-1];
+ return s;
+}
+
+strchrs(s: string, t: string): int
+{
+ for(i := 0; i < len t; i++){
+ p := strchr(s, t[i]);
+ if(p >= 0)
+ return p;
+ }
+ return -1;
+}
+
+strchr(s: string, c: int): int
+{
+ for(i := 0; i < len s; i++)
+ if(s[i] == c)
+ return i;
+ return -1;
+}
+
+deg2rad(d: real): real
+{
+ return d*Pi/180.0;
+}
+
+rad2deg(r: real): real
+{
+ return r*180.0/Pi;
+}
+
+sec2rad(s: real): real
+{
+ return deg2rad(s/3600.0);
+}
+
+norm(r: real): real
+{
+ while(r > Pi)
+ r -= 2.0*Pi;
+ while(r < -Pi)
+ r += 2.0*Pi;
+ return r;
+}
+
+small(r: real): int
+{
+ return r > -Epsilon && r < Epsilon;
+}
+
+trunc(r: real): int
+{
+ # down : assumes r >= 0
+ i := int r;
+ if(real i > r)
+ i--;
+ return i;
+}
+
+abs(x: int): int
+{
+ if(x < 0)
+ return -x;
+ return x;
+}
diff --git a/appl/math/gr.b b/appl/math/gr.b
new file mode 100644
index 00000000..46899353
--- /dev/null
+++ b/appl/math/gr.b
@@ -0,0 +1,557 @@
+implement GR;
+
+include "sys.m";
+ sys: Sys;
+ print, sprint: import sys;
+include "math.m";
+ math: Math;
+ ceil, fabs, floor, Infinity, log10, pow10, sqrt: import math;
+include "draw.m";
+ screen: ref Draw->Screen;
+include "tk.m";
+ tk: Tk;
+ Toplevel: import tk;
+include "tkclient.m";
+ tkclient: Tkclient;
+include "gr.m";
+
+gr_cfg := array[] of {
+ "frame .fc",
+ "frame .fc.b",
+ "label .fc.b.xy -text {0 0} -anchor e",
+ "pack .fc.b.xy -fill x",
+ "pack .fc.b -fill both -expand 1",
+ "canvas .fc.c -relief sunken -bd 2 -width 600 -height 480 -bg white"+
+ " -font /fonts/lucidasans/unicode.8.font",
+ "pack .fc.c -fill both -expand 1",
+ "pack .Wm_t -fill x",
+ "pack .fc -fill both -expand 1",
+ "pack propagate . 0",
+ "bind .fc.c <ButtonPress-1> {send grcmd down1,%x,%y}",
+};
+
+TkCmd(t: ref Toplevel, arg: string): string
+{
+ rv := tk->cmd(t,arg);
+ if(rv!=nil && rv[0]=='!')
+ print("tk->cmd(%s): %s\n",arg,rv);
+ return rv;
+}
+
+
+open(ctxt: ref Draw->Context, title: string): ref Plot
+{
+ if(sys==nil){
+ sys = load Sys Sys->PATH;
+ math = load Math Math->PATH;
+ tk = load Tk Tk->PATH;
+ tkclient = load Tkclient Tkclient->PATH;
+ tkclient->init();
+ }
+ textsize := 8.; # textsize is in points, if no user transform
+ (t, tb) := tkclient->toplevel(ctxt, "", title, Tkclient->Appl);
+ cc := chan of string;
+ tk->namechan(t, cc, "grcmd");
+ p := ref Plot(nil, Infinity,-Infinity,Infinity,-Infinity, textsize, t, tb, cc);
+ for (i:=0; i<len gr_cfg; i++)
+ tk->cmd(p.t,gr_cfg[i]);
+ tkclient->onscreen(p.t, nil);
+ tkclient->startinput(p.t, "kbd"::"ptr"::nil);
+ return p;
+}
+
+Plot.bye(p: self ref Plot)
+{
+ cmdloop: for(;;) alt {
+ s := <-p.t.ctxt.kbd =>
+ tk->keyboard(p.t, s);
+ s := <-p.t.ctxt.ptr =>
+ tk->pointer(p.t, *s);
+ s := <-p.t.ctxt.ctl or
+ s = <-p.t.wreq or
+ s = <-p.titlechan =>
+ if(s == "exit")
+ break cmdloop;
+ tkclient->wmctl(p.t, s);
+ case s{
+ "size" =>
+ canvw := int TkCmd(p.t, ".fc.c cget -width");
+ canvh := int TkCmd(p.t, ".fc.c cget -height");
+ TkCmd(p.t,".fc.b.xy configure -text {"+sprint("%d %d",canvw,canvh)+"}");
+ }
+ press := <-p.canvaschan =>
+ (nil,cmds) := sys->tokenize(press,",");
+ if(cmds==nil) continue;
+ case hd cmds {
+ "down1" =>
+ xpos := real(hd tl cmds);
+ ypos := real(hd tl tl cmds);
+ x := (xpos-bx)/ax;
+ y := -(ypos-tky+by)/ay;
+ TkCmd(p.t,".fc.b.xy configure -text {"+sprint("%.3g %.3g",x,y)+"}");
+ }
+ }
+ TkCmd(p.t,"destroy .;update");
+ p.t = nil;
+}
+
+Plot.equalxy(p: self ref Plot)
+{
+ r := 0.;
+ if( r < p.xmax - p.xmin ) r = p.xmax - p.xmin;
+ if( r < p.ymax - p.ymin ) r = p.ymax - p.ymin;
+ m := (p.xmax + p.xmin)/2.;
+ p.xmax = m + r/2.;
+ p.xmin = m - r/2.;
+ m = (p.ymax + p.ymin)/2.;
+ p.ymax = m + r/2.;
+ p.ymin = m - r/2.;
+}
+
+Plot.graph(p: self ref Plot, x, y: array of real)
+{
+ n := len x;
+ op := OP(GR->GRAPH, n, array[n] of real, array[n] of real, nil);
+ while(n--){
+ t := x[n];
+ op.x[n] = t;
+ if(t < p.xmin)
+ p.xmin = t;
+ if(t > p.xmax)
+ p.xmax = t;
+ t = y[n];
+ op.y[n] = t;
+ if(t < p.ymin)
+ p.ymin = t;
+ if(t > p.ymax)
+ p.ymax = t;
+ }
+ p.op = op :: p.op;
+}
+
+Plot.text(p: self ref Plot, justify: int, s: string, x, y: real)
+{
+ op := OP(GR->TEXT, justify, array[1] of real, array[1] of real, s);
+ op.x[0] = x;
+ op.y[0] = y;
+ p.op = op :: p.op;
+}
+
+Plot.pen(p: self ref Plot, nib: int)
+{
+ p.op = OP(GR->PEN, nib, nil, nil, nil) :: p.op;
+}
+
+
+#---------------------------------------------------------
+# The rest of this file is concerned with sending the "display list"
+# to Tk. The only interesting parts of the problem are picking axes
+# and drawing dashed lines properly.
+
+ax, bx, ay, by: real; # transform user to pixels
+tky: con 630.; # Tk_y = tky - y
+nseg: int; # how many segments in current stroke path
+pendown: int; # is pen currently drawing?
+xoff := array[] of{"w","","e"}; # LJUST, CENTER, RJUST
+yoff := array[] of{"n","","s","s"}; # HIGH, MED, BASE, LOW
+linewidth: real;
+toplevel: ref Toplevel; # p.t
+tkcmd: string;
+
+mv(x, y: real)
+{
+ tkcmd = sprint(".fc.c create line %.1f %.1f", ax*x+bx, tky-(ay*y+by));
+}
+
+stroke()
+{
+ if(pendown){
+ tkcmd += " -width 3"; # -capstyle round -joinstyle round
+ TkCmd(toplevel,tkcmd);
+ tkcmd = nil;
+ pendown = 0;
+ nseg = 0;
+ }
+}
+
+vec(x, y: real)
+{
+ tkcmd += sprint(" %.1f %.1f", ax*x+bx, tky-(ay*y+by));
+ pendown = 1;
+ nseg++;
+ if(nseg>1000){
+ stroke();
+ mv(x,y);
+ }
+}
+
+circle(u, v, radius: real)
+{
+ x := ax*u+bx;
+ y := tky-(ay*v+by);
+ r := radius*(ax+ay)/2.;
+ tkcmd = sprint(".fc.c create oval %.1f %.1f %.1f %.1f -width 3",
+ x-r, y-r, x+r, y+r);
+ TkCmd(toplevel,tkcmd);
+ tkcmd = nil;
+}
+
+text(s: string, x, y: real, xoff, yoff: string)
+{
+ # rot = rotation in degrees. 90 is used for y-axis
+ # x,y are in PostScript coordinate system, not user
+ anchor := yoff + xoff;
+ if(anchor!="")
+ anchor = "-anchor " + anchor + " ";
+ tkcmd = sprint(".fc.c create text %.1f %.1f %s-text '%s",
+ ax*x+bx,
+ tky-(ay*y+by), anchor, s);
+ TkCmd(toplevel,tkcmd);
+ tkcmd = nil;
+}
+
+datarange(xmin, xmax, margin: real): (real,real)
+{
+ r := 1.e-30;
+ if( r < 0.001*fabs(xmin) )
+ r = 0.001*fabs(xmin);
+ if( r < 0.001*fabs(xmax) )
+ r = 0.001*fabs(xmax);
+ if( r < xmax-xmin )
+ r = xmax-xmin;
+ r *= 1.+2.*margin;
+ x0 :=(xmin+xmax)/2. - r/2.;
+ return ( x0, x0 + r);
+}
+
+dashed(ndash: int, x, y: array of real)
+{
+ cx, cy: real; # current position
+ d: real; # length undone in p[i],p[i+1]
+ t: real; # length undone in current dash
+ n := len x;
+ if(n!=len y || n<=0)
+ return;
+
+ # choose precise dashlen
+ s := 0.;
+ for(i := 0; i < n - 1; i += 1){
+ u := x[i+1] - x[i];
+ v := y[i+1] - y[i];
+ s += sqrt(u*u + v*v);
+ }
+ i = int floor(real ndash * s);
+ if(i < 2)
+ i = 2;
+ dashlen := s / real(2 * i - 1);
+
+ t = dashlen;
+ ink := 1;
+ mv(x[0], y[0]);
+ cx = x[0];
+ cy = y[0];
+ for(i = 0; i < n - 1; i += 1){
+ u := x[i+1] - x[i];
+ v := y[i+1] - y[i];
+ d = sqrt(u * u + v * v);
+ if(d > 0.){
+ u /= d;
+ v /= d;
+ while(t <= d){
+ cx += t * u;
+ cy += t * v;
+ if(ink){
+ vec(cx, cy);
+ stroke();
+ }else{
+ mv(cx, cy);
+ }
+ d -= t;
+ t = dashlen;
+ ink = 1 - ink;
+ }
+ cx = x[i+1];
+ cy = y[i+1];
+ if(ink){
+ vec(cx, cy);
+ }else{
+ mv(cx, cy);
+ }
+ t -= d;
+ }
+ }
+ stroke();
+}
+
+labfmt(x:real): string
+{
+ lab := sprint("%.6g",x);
+ if(len lab>2){
+ if(lab[0]=='0' && lab[1]=='.')
+ lab = lab[1:];
+ else if(lab[0]=='-' && len lab>3 && lab[1]=='0' && lab[2]=='.')
+ lab = "-"+lab[2:];
+ }
+ return lab;
+}
+
+Plot.paint(p: self ref Plot, xlabel, xunit, ylabel, yunit: string)
+{
+ oplist: list of OP;
+
+ # tunable parameters for dimensions of graph (fraction of box side)
+ margin: con 0.075; # separation of data from box boundary
+ ticksize := 0.02;
+ sep := ticksize; # separation of text from box boundary
+
+ # derived coordinates of various feature points...
+ x0, x1, y0, y1: real; # box corners, in original coord
+ # radius := 0.2*p.textsize; # radius for circle marker
+ radius := 0.8*p.textsize; # radius for circle marker
+
+ Pen := SOLID;
+ width := SOLID;
+ linewidth = 2.;
+ nseg = 0;
+ pendown = 0;
+
+ if(xunit=="") xunit = nil;
+ if(yunit=="") yunit = nil;
+
+ (x0,x1) = datarange(p.xmin,p.xmax,margin);
+ ax = (400.-2.*p.textsize)/((x1-x0)*(1.+2.*sep));
+ bx = 506.-ax*x1;
+ (y0,y1) = datarange(p.ymin,p.ymax,margin);
+ ay = (400.-2.*p.textsize)/((y1-y0)*(1.+2.*sep));
+ by = 596.-ay*y1;
+ # PostScript version
+ # magic numbers here come from BoundingBox: 106 196 506 596
+ # (x0,x1) = datarange(p.xmin,p.xmax,margin);
+ # ax = (400.-2.*p.textsize)/((x1-x0)*(1.+2.*sep));
+ # bx = 506.-ax*x1;
+ # (y0,y1) = datarange(p.ymin,p.ymax,margin);
+ # ay = (400.-2.*p.textsize)/((y1-y0)*(1.+2.*sep));
+ # by = 596.-ay*y1;
+
+ # convert from fraction of box to PostScript units
+ ticksize *= ax*(x1-x0);
+ sep *= ax*(x1-x0);
+
+ # revert to original drawing order
+ log := p.op;
+ oplist = nil;
+ while(log!=nil){
+ oplist = hd log :: oplist;
+ log = tl log;
+ }
+ p.op = oplist;
+
+ toplevel = p.t;
+ #------------send display list to Tk-----------------
+ while(oplist!=nil){
+ op := hd oplist;
+ n := op.n;
+ case op.code{
+ GRAPH =>
+ if(Pen == DASHED){
+ dashed(17, op.x, op.y);
+ }else if(Pen == DOTTED){
+ dashed(85, op.x, op.y);
+ }else{
+ for(i:=0; i<n; i++){
+ xx := op.x[i];
+ yy := op.y[i];
+ if(Pen == CIRCLE){
+ circle(xx, yy, radius/(ax+ay));
+ }else if(Pen == CROSS){
+ mv(xx-radius/ax, yy);
+ vec(xx+radius/ax, yy);
+ stroke();
+ mv(xx, yy-radius/ay);
+ vec(xx, yy+radius/ay);
+ stroke();
+ }else if(Pen == INVIS){
+ }else{
+ if(i==0){
+ mv(xx, yy);
+ }else{
+ vec(xx, yy);
+ }
+ }
+ }
+ stroke();
+ }
+ TEXT =>
+ angle := 0.;
+ if(op.n&UP) angle = 90.;
+ text(op.t,op.x[0],op.y[0],xoff[n&7],yoff[(n>>3)&7]);
+ PEN =>
+ Pen = n;
+ if( Pen==SOLID && width!=SOLID ){
+ linewidth = 2.;
+ width=SOLID;
+ }else if( Pen==REFERENCE && width!=REFERENCE ){
+ linewidth = 0.8;
+ width=REFERENCE;
+ }
+ }
+ oplist = tl oplist;
+ }
+
+ #--------------------now add axes-----------------------
+ mv(x0,y0);
+ vec(x1,y0);
+ vec(x1,y1);
+ vec(x0,y1);
+ vec(x0,y0);
+ stroke();
+
+ # x ticks
+ (lab1,labn,labinc,k,u,s) := mytic(x0,x1);
+ for (i := lab1; i <= labn; i += labinc){
+ r := real i*s*u;
+ mv(r,y0);
+ vec(r,y0+ticksize/ay);
+ stroke();
+ mv(r,y1);
+ vec(r,y1-ticksize/ay);
+ stroke();
+ text(labfmt(real i*s),r,y0-sep/ay,"","n");
+ }
+ yy := y0-(2.*sep+p.textsize)/ay;
+ labelstr := "";
+ if(xlabel!=nil)
+ labelstr = xlabel;
+ if(k!=0||xunit!=nil)
+ labelstr += " /";
+ if(k!=0)
+ labelstr += " ₁₀"+ string k;
+ if(xunit!=nil)
+ labelstr += " " + xunit;
+ text(labelstr,(x0+x1)/2.,yy,"","n");
+
+ # y ticks
+ (lab1,labn,labinc,k,u,s) = mytic(y0,y1);
+ for (i = lab1; i <= labn; i += labinc){
+ r := real i*s*u;
+ mv(x0,r);
+ vec(x0+ticksize/ax,r);
+ stroke();
+ mv(x1,r);
+ vec(x1-ticksize/ax,r);
+ stroke();
+ text(labfmt(real i*s),x0-sep/ax,r,"e","");
+ }
+ xx := x0-(4.*sep+p.textsize)/ax;
+ labelstr = "";
+ if(ylabel!=nil)
+ labelstr = ylabel;
+ if(k!=0||yunit!=nil)
+ labelstr += " /";
+ if(k!=0)
+ labelstr += " ₁₀"+ string k;
+ if(yunit!=nil)
+ labelstr += " " + yunit;
+ text(labelstr,xx,(y0+y1)/2.,"e","");
+
+ TkCmd(p.t, "update");
+}
+
+
+
+# automatic tic choice Eric Grosse 9 Dec 84
+# Input: low and high endpoints of expanded data range
+# Output: lab1, labn, labinc, k, u, s where the tics are
+# (lab1*s, (lab1+labinc)*s, ..., labn*s) * 10^k
+# and u = 10^k. k is metric, i.e. k=0 mod 3.
+
+max3(a, b, c: real): real
+{
+ if(a<b) a=b;
+ if(a<c) a=c;
+ return(a);
+}
+
+my_mod(i, n: int): int
+{
+ while(i< 0) i+=n;
+ while(i>=n) i-=n;
+ return(i);
+}
+
+mytic(l, h: real): (int,int,int,int,real,real)
+{
+ lab1, labn, labinc, k, nlab, j, ndig, t1, tn: int;
+ u, s: real;
+ eps := .0001;
+ k = int floor( log10((h-l)/(3.+eps)) );
+ u = pow10(k);
+ t1 = int ceil(l/u-eps);
+ tn = int floor(h/u+eps);
+ lab1 = t1;
+ labn = tn;
+ labinc = 1;
+ nlab = labn - lab1 + 1;
+ if( nlab>5 ){
+ lab1 = t1 + my_mod(-t1,2);
+ labn = tn - my_mod( tn,2);
+ labinc = 2;
+ nlab = (labn-lab1)/labinc + 1;
+ if( nlab>5 ){
+ lab1 = t1 + my_mod(-t1,5);
+ labn = tn - my_mod( tn,5);
+ labinc = 5;
+ nlab = (labn-lab1)/labinc + 1;
+ if( nlab>5 ){
+ u *= 10.;
+ k++;
+ lab1 = int ceil(l/u-eps);
+ labn = int floor(h/u+eps);
+ nlab = labn - lab1 + 1;
+ labinc = 1;
+ } else if( nlab<3 ){
+ lab1 = t1 + my_mod(-t1,4);
+ labn = tn - my_mod( tn,4);
+ labinc = 4;
+ nlab = (labn-lab1)/labinc + 1;
+ }
+ }
+ }
+ ndig = int(1.+floor(log10(max3(fabs(real lab1),fabs(real labn),1.e-30))));
+ if( ((k<=0)&&(k>=-ndig)) # no zeros have to be added
+ || ((k<0)&&(k>=-3))
+ || ((k>0)&&(ndig+k<=4)) ){ # even with zeros, label is small
+ s = u;
+ k = 0;
+ u = 1.;
+ }else if(k>0){
+ s = 1.;
+ j = ndig;
+ while(k%3!=0){
+ k--;
+ u/=10.;
+ s*=10.;
+ j++;
+ }
+ if(j-3>0){
+ k+=3;
+ u*=1000.;
+ s/=1000.;
+ }
+ }else{ # k<0
+ s = 1.;
+ j = ndig;
+ while(k%3!=0){
+ k++;
+ u*=10.;
+ s/=10.;
+ j--;
+ }
+ if(j<0){
+ k-=3;
+ u/=1000.;
+ s*=1000.;
+ }
+ }
+ return (lab1, labn, labinc, k, u, s);
+}
diff --git a/appl/math/graph0.b b/appl/math/graph0.b
new file mode 100644
index 00000000..1e016b51
--- /dev/null
+++ b/appl/math/graph0.b
@@ -0,0 +1,84 @@
+implement Graph0;
+
+include "sys.m";
+ sys: Sys;
+ print: import sys;
+
+include "draw.m";
+include "tk.m";
+ tk: Tk;
+
+include "bufio.m";
+ bufio: Bufio;
+ Iobuf: import bufio;
+
+include "gr.m";
+ gr: GR;
+ Plot: import gr;
+
+Graph0: module{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+
+plotfile(ctxt: ref Draw->Context, nil: list of string, filename: string){
+ p := gr->open(ctxt,filename);
+ input := bufio->open(filename,bufio->OREAD);
+ if(input==nil){
+ print("can't read %s",filename);
+ exit;
+ }
+
+ n := 0;
+ maxn := 100;
+ x := array[maxn] of real;
+ y := array[maxn] of real;
+ while(1){
+ xn := input.gett(" \t\n\r");
+ if(xn==nil)
+ break;
+ yn := input.gett(" \t\n\r");
+ if(yn==nil){
+ print("after reading %d pairs, saw singleton\n",n);
+ exit;
+ }
+ if(n>=maxn){
+ maxn *= 2;
+ newx := array[maxn] of real;
+ newy := array[maxn] of real;
+ for(i:=0; i<n; i++){
+ newx[i] = x[i];
+ newy[i] = y[i];
+ }
+ x = newx;
+ y = newy;
+ }
+ x[n] = real xn;
+ y[n] = real yn;
+ n++;
+ }
+ if(n==0){
+ print("empty input\n");
+ exit;
+ }
+
+ p.graph(x[0:n],y[0:n]);
+ p.pen(GR->CIRCLE);
+ p.graph(x[0:n],y[0:n]);
+ p.paint("",nil,"",nil);
+ p.bye();
+}
+
+init(ctxt: ref Draw->Context, argv: list of string){
+ sys = load Sys Sys->PATH;
+ tk = load Tk Tk->PATH;
+ bufio = load Bufio Bufio->PATH;
+ if((gr = load GR GR->PATH) == nil){
+ sys->print("%s: Can't load gr\n",hd argv);
+ exit;
+ }
+
+ argv = tl argv;
+ if(argv!=nil)
+ plotfile(ctxt,argv,hd argv);
+}
diff --git a/appl/math/hist0.b b/appl/math/hist0.b
new file mode 100644
index 00000000..0484c98b
--- /dev/null
+++ b/appl/math/hist0.b
@@ -0,0 +1,99 @@
+implement Hist0;
+
+include "sys.m";
+ sys: Sys;
+ print: import sys;
+include "math.m";
+ math: Math;
+ fmin: import math;
+
+include "draw.m";
+include "tk.m";
+ tk: Tk;
+
+include "bufio.m";
+ bufio: Bufio;
+ Iobuf: import bufio;
+
+include "gr.m";
+ gr: GR;
+ Plot: import gr;
+
+Hist0: module{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+
+plotfile(ctxt: ref Draw->Context, nil: list of string, filename: string){
+ p := gr->open(ctxt,filename);
+ input := bufio->open(filename,bufio->OREAD);
+ if(input==nil){
+ print("can't read %s",filename);
+ exit;
+ }
+
+ n := 0;
+ maxn := 100;
+ x := array[maxn] of real;
+ y := array[maxn] of real;
+ while(1){
+ xn := input.gett(" \t\n\r");
+ if(xn==nil)
+ break;
+ yn := input.gett(" \t\n\r");
+ if(yn==nil){
+ print("after reading %d pairs, saw singleton\n",n);
+ exit;
+ }
+ if(n>=maxn){
+ maxn *= 2;
+ newx := array[maxn] of real;
+ newy := array[maxn] of real;
+ for(i:=0; i<n; i++){
+ newx[i] = x[i];
+ newy[i] = y[i];
+ }
+ x = newx;
+ y = newy;
+ }
+ x[n] = real xn;
+ y[n] = real yn;
+ n++;
+ }
+ if(n==0){
+ print("empty input\n");
+ exit;
+ }
+
+ for(i:=0; i<n; i++){
+ h := 0.2;
+ if(i==0) h *= x[i+1]-x[i];
+ else if(i==n-1) h *= x[i]-x[i-1];
+ else h *= fmin(x[i+1]-x[i],x[i]-x[i-1]);
+ barx := array[] of{x[i]-h,x[i]-h,x[i]+h,x[i]+h,x[i]-h};
+ bary := array[] of{0.,y[i],y[i],0.,0.};
+ p.graph(barx,bary);
+ }
+ p.paint("",nil,"",nil);
+ p.bye();
+}
+
+init(ctxt: ref Draw->Context, argv: list of string){
+ sys = load Sys Sys->PATH;
+ math = load Math Math->PATH;
+ tk = load Tk Tk->PATH;
+ bufio = load Bufio Bufio->PATH;
+ if((gr = load GR GR->PATH) == nil){
+ sys->print("%s: Can't load gr\n",hd argv);
+ exit;
+ }
+
+ tkargs := "";
+ argv = tl argv;
+ if(argv != nil) {
+ tkargs = hd argv;
+ argv = tl argv;
+ }
+ if(argv!=nil)
+ plotfile(ctxt,argv,hd argv);
+}
diff --git a/appl/math/linalg.b b/appl/math/linalg.b
new file mode 100644
index 00000000..5477c4ca
--- /dev/null
+++ b/appl/math/linalg.b
@@ -0,0 +1,197 @@
+implement LinAlg;
+
+include "sys.m";
+sys: Sys;
+print: import sys;
+
+include "math.m";
+math: Math;
+ceil, fabs, floor, Infinity, log10, pow10, sqrt: import math;
+dot, gemm, iamax: import math;
+
+include "linalg.m";
+
+# print a matrix in MATLAB-compatible format
+printmat(label:string, a:array of real, lda, m, n:int)
+{
+ if(m>30 || n>10)
+ return;
+ if(sys==nil){
+ sys = load Sys Sys->PATH;
+ math = load Math Math->PATH;
+ }
+ print("%% %d by %d matrix\n",m,n);
+ print("%s = [",label);
+ for(i:=0; i<m; i++){
+ print("%.4g",a[i]);
+ for(j:=1; j<n; j++)
+ print(", %.4g",a[i+lda*j]);
+ if(i==m-1)
+ print("]\n");
+ else
+ print(";\n");
+ }
+}
+
+
+# Constant times a vector plus a vector.
+daxpy(da:real, dx:array of real, dy:array of real)
+{
+ n := len dx;
+ gemm('N','N',n,1,n,da,nil,0,dx,n,1.,dy,n);
+}
+
+# Scales a vector by a constant.
+dscal(da:real, dx:array of real)
+{
+ n := len dx;
+ gemm('N','N',n,1,n,0.,nil,0,nil,0,da,dx,n);
+}
+
+# gaussian elimination with partial pivoting
+# dgefa factors a double precision matrix by gaussian elimination.
+# dgefa is usually called by dgeco, but it can be called
+# directly with a saving in time if rcond is not needed.
+# (time for dgeco) = (1 + 9/n)*(time for dgefa) .
+# on entry
+# a REAL precision[n][lda]
+# the matrix to be factored.
+# lda integer
+# the leading dimension of the array a .
+# n integer
+# the order of the matrix a .
+# on return
+# a an upper triangular matrix and the multipliers
+# which were used to obtain it.
+# the factorization can be written a = l*u where
+# l is a product of permutation and unit lower
+# triangular matrices and u is upper triangular.
+# ipvt integer[n]
+# an integer vector of pivot indices.
+# info integer
+# = 0 normal value.
+# = k if u[k][k] .eq. 0.0 . this is not an error
+# condition for this subroutine, but it does
+# indicate that dgesl or dgedi will divide by zero
+# if called. use rcond in dgeco for a reliable
+# indication of singularity.
+dgefa(a:array of real, lda, n:int, ipvt:array of int): int
+{
+ if(sys==nil){
+ sys = load Sys Sys->PATH;
+ math = load Math Math->PATH;
+ }
+ info := 0;
+ nm1 := n - 1;
+ if(nm1 >= 0)
+ for(k := 0; k < nm1; k++){
+ kp1 := k + 1;
+ ldak := lda*k;
+
+ # find l = pivot index
+ l := iamax(a[ldak+k:ldak+n]) + k;
+ ipvt[k] = l;
+
+ # zero pivot implies this column already triangularized
+ if(a[ldak+l]!=0.){
+
+ # interchange if necessary
+ if(l!=k){
+ t := a[ldak+l];
+ a[ldak+l] = a[ldak+k];
+ a[ldak+k] = t;
+ }
+
+ # compute multipliers
+ t := -1./a[ldak+k];
+ dscal(t,a[ldak+k+1:ldak+n]);
+
+ # row elimination with column indexing
+ for(j := kp1; j < n; j++){
+ ldaj := lda*j;
+ t = a[ldaj+l];
+ if(l!=k){
+ a[ldaj+l] = a[ldaj+k];
+ a[ldaj+k] = t;
+ }
+ daxpy(t,a[ldak+k+1:ldak+n],a[ldaj+k+1:ldaj+n]);
+ }
+ }else
+ info = k;
+ }
+ ipvt[n-1] = n-1;
+ if(a[lda*(n-1)+(n-1)] == 0.)
+ info = n-1;
+ return info;
+}
+
+
+# dgesl solves the double precision system
+# a * x = b or trans(a) * x = b
+# using the factors computed by dgeco or dgefa.
+# on entry
+# a double precision[n][lda]
+# the output from dgeco or dgefa.
+# lda integer
+# the leading dimension of the array a .
+# n integer
+# the order of the matrix a .
+# ipvt integer[n]
+# the pivot vector from dgeco or dgefa.
+# b double precision[n]
+# the right hand side vector.
+# job integer
+# = 0 to solve a*x = b ,
+# = nonzero to solve trans(a)*x = b where
+# trans(a) is the transpose.
+# on return
+# b the solution vector x .
+# error condition
+# a division by zero will occur if the input factor contains a
+# zero on the diagonal. technically this indicates singularity
+# but it is often caused by improper arguments or improper
+# setting of lda.
+dgesl(a:array of real, lda, n:int, ipvt:array of int, b:array of real, job:int)
+{
+ nm1 := n - 1;
+ if(job == 0){ # job = 0 , solve a * x = b
+ # first solve l*y = b
+ if(nm1 >= 1)
+ for(k := 0; k < nm1; k++){
+ l := ipvt[k];
+ t := b[l];
+ if(l!=k){
+ b[l] = b[k];
+ b[k] = t;
+ }
+ daxpy(t,a[lda*k+k+1:lda*k+n],b[k+1:n]);
+ }
+
+ # now solve u*x = y
+ for(kb := 0; kb < n; kb++){
+ k = n - (kb + 1);
+ b[k] = b[k]/a[lda*k+k];
+ t := -b[k];
+ daxpy(t,a[lda*k:lda*k+k],b[0:k]);
+ }
+ }else{ # job = nonzero, solve trans(a) * x = b
+ # first solve trans(u)*y = b
+ for(k := 0; k < n; k++){
+ t := dot(a[lda*k:lda*k+k],b[0:k]);
+ b[k] = (b[k] - t)/a[lda*k+k];
+ }
+
+ # now solve trans(l)*x = y
+ if(nm1 >= 1)
+ for(kb := 1; kb < nm1; kb++){
+ k = n - (kb+1);
+ b[k] += dot(a[lda*k+k+1:lda*k+n],b[k+1:n]);
+ l := ipvt[k];
+ if(l!=k){
+ t := b[l];
+ b[l] = b[k];
+ b[k] = t;
+ }
+ }
+ }
+}
diff --git a/appl/math/linbench.b b/appl/math/linbench.b
new file mode 100644
index 00000000..1857093a
--- /dev/null
+++ b/appl/math/linbench.b
@@ -0,0 +1,197 @@
+# Translated to Limbo by Eric Grosse <ehg@netlib.bell-labs.com> 3/96
+# Translated to Java by Reed Wade (wade@cs.utk.edu) 2/96
+# Translated to C by Bonnie Toy 5/88
+# Will Menninger, 10/93
+# Jack Dongarra, linpack, 3/11/78.
+# Cleve Moler, linpack, 08/14/78
+
+implement linbench;
+
+include "sys.m";
+ sys: Sys;
+ print: import sys;
+include "math.m";
+ math: Math;
+ dot, fabs, gemm, iamax: import math;
+include "draw.m";
+ draw: Draw;
+ Rect, Screen, Display, Image: import draw;
+include "tk.m";
+ tk: Tk;
+ Toplevel: import tk;
+include "tkclient.m";
+ tkclient: Tkclient;
+include "bufio.m";
+ bufio: Bufio;
+ Iobuf: import bufio;
+include "linalg.m";
+ linalg: LinAlg;
+ dgefa, dgesl, printmat: import linalg;
+
+ctxt: ref Draw->Context;
+top: ref Tk->Toplevel;
+buttonlist: string;
+
+
+linbench: module
+{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+tkcmd(arg: string): string{
+ rv := tk->cmd(top,arg);
+ if(rv!=nil && rv[0]=='!')
+ print("tk->cmd(%s): %s\n",arg,rv);
+ return rv;
+}
+
+init(xctxt: ref Draw->Context, nil: list of string)
+{
+ sys = load Sys Sys->PATH;
+ math = load Math Math->PATH;
+ draw = load Draw Draw->PATH;
+ tk = load Tk Tk->PATH;
+ tkclient = load Tkclient Tkclient->PATH;
+ linalg = load LinAlg LinAlg->PATH;
+ if(linalg==nil) print("couldn't load LinAlg\n");
+ sys->pctl(Sys->NEWPGRP, nil);
+ tkclient->init();
+ ctxt = xctxt;
+ menubut: chan of string;
+ (top, menubut) = tkclient->toplevel(ctxt, "",
+ "Linpack in Limbo", Tkclient->Appl);
+ cmd := chan of string;
+ tk->namechan(top, cmd, "cmd");
+
+ tkcmd("pack .Wm_t -fill x");
+ tkcmd("frame .b");
+ tkcmd("button .b.Run -text Run -command {send cmd run}");
+ tkcmd("entry .b.N -width 10w"); tkcmd(".b.N insert 0 200");
+ tkcmd("pack .b.Run .b.N -side left");
+ tkcmd("pack .b -anchor w");
+ tkcmd("frame .d");
+ tkcmd("listbox .d.f -width 35w -height 150 -selectmode single -yscrollcommand {.d.fscr set}");
+ tkcmd("scrollbar .d.fscr -command {.d.f yview}");
+ tkcmd("pack .d.f .d.fscr -expand 1 -fill y -side right");
+ tkcmd("pack .d -side top");
+ tkcmd("focus .b.N");
+ tkcmd("pack propagate . 0;update");
+ tkclient->onscreen(top, nil);
+ tkclient->startinput(top, "kbd"::"ptr"::nil);
+
+ for(;;) alt {
+ s := <-top.ctxt.kbd =>
+ tk->keyboard(top, s);
+ s := <-top.ctxt.ptr =>
+ tk->pointer(top, *s);
+ s := <-top.ctxt.ctl or
+ s = <-top.wreq or
+ s = <-menubut =>
+ tkclient->wmctl(top, s);
+ press := <-cmd =>
+ case press {
+ "run" =>
+ tkcmd("cursor -bitmap cursor.wait; update");
+ nstr := tkcmd(".b.N get");
+ n := int nstr;
+ (mflops,secs) := benchmark(n);
+ result := sys->sprint("%8.2f Mflops %8.1f secs",mflops,secs);
+ tkcmd("cursor -default");
+ tkcmd(".d.f insert end {" + result + "}");
+ tkcmd(".d.f yview moveto 1; update");
+ }
+ }
+}
+
+benchmark(n: int): (real,real)
+{
+ math = load Math Math->PATH;
+
+ time := array [2] of real;
+ lda := 201;
+ if(n>lda) lda = n;
+ a := array [lda*n] of real;
+ b := array [n] of real;
+ x := array [n] of real;
+ ipvt := array [n] of int;
+ ops := (2*n*n*n)/3 + 2*n*n;
+
+ norma := matgen(a,lda,n,b);
+ printmat("a",a,lda,n,n);
+ printmat("b",b,lda,n,1);
+ t1 := second();
+ dgefa(a,lda,n,ipvt);
+ time[0] = second() - t1;
+ printmat("a",a,lda,n,n);
+ t1 = second();
+ dgesl(a,lda,n,ipvt,b,0);
+ time[1] = second() - t1;
+ total := time[0] + time[1];
+
+ for(i := 0; i < n; i++) {
+ x[i] = b[i];
+ }
+ printmat("x",x,lda,n,1);
+ norma = matgen(a,lda,n,b);
+ for(i = 0; i < n; i++) {
+ b[i] = -b[i];
+ }
+ dmxpy(b,x,a,lda);
+ resid := 0.;
+ normx := 0.;
+ for(i = 0; i < n; i++){
+ if(resid<fabs(b[i])) resid = fabs(b[i]);
+ if(normx<fabs(x[i])) normx = fabs(x[i]);
+ }
+
+ eps_result := math->MachEps;
+ residn_result := (real n)*norma*normx*eps_result;
+ if(residn_result!=0.)
+ residn_result = resid/residn_result;
+ else
+ print("can't scale residual.");
+ if(residn_result>math->sqrt(real n))
+ print("resid/MachEps=%.3g\n",residn_result);
+ time_result := total;
+ mflops_result := 0.;
+ if(total!=0.)
+ mflops_result = real ops/(1e6*total);
+ else
+ print("can't measure time\n");
+ return (mflops_result,time_result);
+}
+
+
+# multiply matrix m times vector x and add the r_result to vector y.
+dmxpy(y, x, m:array of real, ldm: int)
+{
+ n1 := len y;
+ n2 := len x;
+ gemm('N','N',n1,1,n2,1.,m,ldm,x,n2,1.,y,n1);
+}
+
+second(): real
+{
+ return(real sys->millisec()/1000.);
+}
+
+
+# generate a (fixed) random matrix and right hand side
+# a[i][j] => a[lda*i+j]
+matgen(a: array of real, lda, n: int, b: array of real): real
+{
+ seed := 1325;
+ norma := 0.;
+ for(j := 0; j < n; j++)
+ for(i := 0; i < n; i++){
+ seed = 3125*seed % 65536;
+ a[lda*j+i] = (real seed - 32768.0)/16384.0;
+ if(norma < a[lda*j+i]) norma = a[lda*j+i];
+ }
+ for (i = 0; i < n; i++)
+ b[i] = 0.;
+ for (j = 0; j < n; j++)
+ for (i = 0; i < n; i++)
+ b[i] += a[lda*j+i];
+ return norma;
+}
diff --git a/appl/math/mersenne.b b/appl/math/mersenne.b
new file mode 100644
index 00000000..c90254b0
--- /dev/null
+++ b/appl/math/mersenne.b
@@ -0,0 +1,99 @@
+implement Mersenne;
+
+include "sys.m";
+ sys : Sys;
+include "draw.m";
+include "keyring.m";
+ keyring: Keyring;
+ IPint: import keyring;
+
+# Test primality of Mersenne numbers
+
+Mersenne: module
+{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+init(nil: ref Draw->Context, argv: list of string)
+{
+ sys = load Sys Sys->PATH;
+ keyring = load Keyring Keyring->PATH;
+ p := 3;
+ if(tl argv != nil)
+ p = int hd tl argv;
+ if(isprime(p) && (p == 2 || lucas(p)))
+ s := "";
+ else
+ s = "not ";
+ sys->print("2^%d-1 is %sprime\n", p, s);
+}
+
+# s such that s^2 <= n
+sqrt(n: int): int
+{
+ v := n;
+ r := 0;
+ for(t := 1<<30; t; t >>= 2){
+ if(t+r <= v){
+ v -= t+r;
+ r = (r>>1)|t;
+ }
+ else
+ r = r>>1;
+ }
+ return r;
+}
+
+isprime(n: int): int
+{
+ if(n < 2)
+ return 0;
+ if(n == 2)
+ return 1;
+ if((n&1) == 0)
+ return 0;
+ s := sqrt(n);
+ for(i := 3; i <= s; i += 2)
+ if(n%i == 0)
+ return 0;
+ return 1;
+}
+
+pow(b : ref IPint, n : int): ref IPint
+{
+ zero := IPint.inttoip(0);
+ one := IPint.inttoip(1);
+ if((b.cmp(zero) == 0 && n != 0) || b.cmp(one) == 0 || n == 1)
+ return b;
+ if(n == 0)
+ return one;
+ c := b;
+ b = one;
+ while(n){
+ while(!(n & 1)){
+ n >>= 1;
+ c = c.mul(c);
+ }
+ n--;
+ b = c.mul(b);
+ }
+ return b;
+}
+
+lucas(p: int): int
+{
+ zero := IPint.inttoip(0);
+ one := IPint.inttoip(1);
+ two := IPint.inttoip(2);
+ bigp := pow(two, p).sub(one);
+ u := IPint.inttoip(4);
+ for(i := 2; i < p; i++){
+ u = u.mul(u);
+ if(u.cmp(two) <= 0)
+ u = two.sub(u);
+ else
+ u = u.sub(two).expmod(one, bigp);
+ }
+ return u.cmp(zero) == 0;
+}
+
diff --git a/appl/math/mkfile b/appl/math/mkfile
new file mode 100644
index 00000000..7d21050a
--- /dev/null
+++ b/appl/math/mkfile
@@ -0,0 +1,43 @@
+<../../mkconfig
+
+TARG=\
+ ack.dis\
+ crackerbarrel.dis\
+ factor.dis\
+ ffts.dis\
+ fibonacci.dis\
+ fit.dis\
+ genprimes.dis\
+ geodesy.dis\
+ graph0.dis\
+ gr.dis\
+ hist0.dis\
+ linalg.dis\
+ linbench.dis\
+ mersenne.dis\
+ parts.dis\
+ perms.dis\
+ pi.dis\
+ polyfill.dis\
+ polyhedra.dis\
+ powers.dis\
+ primes.dis\
+ sieve.dis\
+
+MODULES=
+
+SYSMODULES=\
+ math.m\
+ gr.m\
+ linalg.m\
+ draw.m\
+ sys.m\
+ tk.m\
+ wmlib.m\
+ bufio.m\
+ math/polyfill.m\
+ math/polyhedra.m\
+
+DISBIN=$ROOT/dis/math
+
+<$ROOT/mkfiles/mkdis
diff --git a/appl/math/parts.b b/appl/math/parts.b
new file mode 100644
index 00000000..36c36d8e
--- /dev/null
+++ b/appl/math/parts.b
@@ -0,0 +1,125 @@
+implement Partitions;
+
+include "sys.m";
+ sys : Sys;
+include "draw.m";
+include "keyring.m";
+ keyring: Keyring;
+ IPint: import keyring;
+
+#
+# the number p(n) of partitions of n
+# based upon the formula :-
+# p(n) = p(n-1)+p(n-2)-p(n-5)-p(n-7)+p(n-12)+p(n-15)-p(n-22)-p(n-26)+.....
+# where p[0] = 1 and p[m] = 0 for m < 0
+#
+
+aflag := 0;
+cflag := 0;
+
+Partitions: module
+{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+init(nil: ref Draw->Context, argv: list of string)
+{
+ sys = load Sys Sys->PATH;
+ keyring = load Keyring Keyring->PATH;
+ argv = tl argv;
+ while(argv != nil){
+ s := hd argv;
+ if(s != nil && s[0] == '-'){
+ for(i := 1; i < len s; i++){
+ case s[i]{
+ 'a' => aflag = 1;
+ 'c' => cflag = 1;
+ }
+ }
+ }
+ else
+ parts(int s);
+ argv = tl argv;
+ }
+}
+
+parts(m : int)
+{
+ if (aflag)
+ sys->print("n p(n)\n");
+ if (m <= 0) {
+ p := 0;
+ if (m == 0)
+ p = 1;
+ if (aflag)
+ sys->print("%d %d\n", m, p);
+ else
+ sys->print("p[%d] = %d\n", m, p);
+ return;
+ }
+ p := array[m+1] of ref IPint;
+ if (p == nil)
+ return;
+ p[0] = IPint.inttoip(1);
+ for (i := 1; i <= m; i++) {
+ k := i;
+ s := 1;
+ n := IPint.inttoip(0);
+ for (j := 1; ; j++) {
+ k -= 2*j-1;
+ if (k < 0)
+ break;
+ if (s == 1)
+ n = n.add(p[k]);
+ else
+ n = n.sub(p[k]);
+ k -= j;
+ if (k < 0)
+ break;
+ if (s == 1)
+ n = n.add(p[k]);
+ else
+ n = n.sub(p[k]);
+ s = -s;
+ }
+ if (aflag)
+ sys->print("%d %s\n", i, n.iptostr(10));
+ p[i] = n;
+ }
+ if (!aflag)
+ sys->print("p[%d] = %s\n", m, p[m].iptostr(10));
+ if (cflag)
+ check(m, p);
+}
+
+#
+# given p[0]..p[m], search for congruences of the form
+# p[ni+j] = r mod i
+#
+check(m : int, p : array of ref IPint)
+{
+ one := IPint.inttoip(1);
+ for (i := 2; i < m/3; i++) {
+ ip := IPint.inttoip(i);
+ for (j := 0; j < i; j++) {
+ k := j;
+ r := p[k].expmod(one, ip).iptoint();
+ s := 1;
+ for (;;) {
+ k += i;
+ if (k > m)
+ break;
+ if (p[k].expmod(one, ip).iptoint() != r) {
+ r = -1;
+ break;
+ }
+ s++;
+ }
+ if (r >= 0)
+ if (j == 0)
+ sys->print("p(%dm) = %d mod %d ?\n", i, r, i);
+ else
+ sys->print("p(%dm+%d) = %d mod %d ?\n", i, j, r, i);
+ }
+ }
+}
diff --git a/appl/math/perms.b b/appl/math/perms.b
new file mode 100644
index 00000000..23c5e52b
--- /dev/null
+++ b/appl/math/perms.b
@@ -0,0 +1,131 @@
+#
+# initially generated by c2l
+#
+
+implement Perms;
+
+include "draw.m";
+
+Perms: module
+{
+ init: fn(nil: ref Draw->Context, argl: list of string);
+};
+
+include "sys.m";
+ sys: Sys;
+
+init(nil: ref Draw->Context, argl: list of string)
+{
+ sys = load Sys Sys->PATH;
+ main(len argl, argl);
+}
+
+#
+# * generate permutations of N elements
+# * from ``On Programming, an interim report on the SETL project'',
+# * Jacob T Schwartz (ed), New York University
+#
+Seq: adt{
+ nel: int;
+ el: array of int;
+};
+
+origin: int = 1;
+
+main(argc: int, argv: list of string): int
+{
+ n: int;
+
+ if(argc > 1 && (as := hd tl argv)[0] == '-'){
+ origin = int (as[1: ]);
+ argc--;
+ argv = tl argv;
+ }
+ if(argc != 2){
+ sys->fprint(sys->fildes(2), "Usage: perms #elements\n");
+ exit;
+ }
+ n = int hd tl argv;
+ if(n > 0)
+ perms(n);
+ exit;
+}
+
+
+perms(n: int)
+{
+ seq: ref Seq;
+
+ seq = newseq(n);
+ do
+ putseq(seq);
+ while(eperm(seq) != nil);
+}
+
+putseq(seq: ref Seq)
+{
+ k: int;
+
+ for(k = 0; k < seq.nel; k++)
+ sys->print(" %d", seq.el[k]+origin);
+ sys->print("\n");
+}
+
+eperm(seq: ref Seq): ref Seq
+{
+ j, k, n: int;
+
+ n = seq.nel;
+ # if sequence is monotone decreasing, there are no more
+ # permutations. Otherwise, find last point of increase
+ hit := 0;
+ for(j = n-2; j >= 0; j--)
+ if(seq.el[j] < seq.el[j+1]){
+ hit = 1;
+ break;
+ }
+ if(!hit)
+ return nil;
+ # then find the last seq[k] which exceeds seq[j] and swop
+ for(k = seq.nel-1; k > j; k--)
+ if(seq.el[k] > seq.el[j]){
+ {
+ t: int;
+
+ t = seq.el[k];
+ seq.el[k] = seq.el[j];
+ seq.el[j] = t;
+ }
+ ;
+ break;
+ }
+ # then re-arrange all the elements from seq[j+1] into
+ # increasing order
+ for(k = j+1; k < (n+j+1)/2; k++){
+ kk: int;
+
+ kk = n-k+j;
+ {
+ t: int;
+
+ t = seq.el[k];
+ seq.el[k] = seq.el[kk];
+ seq.el[kk] = t;
+ }
+ ;
+ }
+ return seq;
+}
+
+newseq(n: int): ref Seq
+{
+ seq: ref Seq;
+ k: int;
+
+ seq = ref Seq;
+ seq.nel = n;
+ seq.el = array[n] of int;
+ for(k = 0; k < n; k++)
+ seq.el[k] = k;
+ return seq;
+}
diff --git a/appl/math/pi.b b/appl/math/pi.b
new file mode 100644
index 00000000..2ddd1a14
--- /dev/null
+++ b/appl/math/pi.b
@@ -0,0 +1,405 @@
+implement Pi;
+
+include "sys.m";
+ sys: Sys;
+include "draw.m";
+include "math.m";
+ math: Math;
+ log: import math;
+include "daytime.m";
+ daytime: Daytime;
+
+LBASE: con 3; # 4
+BASE: con 1000; # 10000
+
+stderr: ref Sys->FD;
+
+# spawn process for each series ?
+
+Pi: module
+{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+init(nil: ref Draw->Context, argv: list of string)
+{
+ sys = load Sys Sys->PATH;
+ math = load Math Math->PATH;
+ daytime = load Daytime Daytime->PATH;
+
+ stderr = sys->fildes(2);
+ dp := 1000;
+ argv = tl argv;
+ if(argv != nil){
+ if(tl argv != nil){
+ picmp(hd argv, hd tl argv);
+ exit;
+ }
+ dp = int hd argv;
+ }
+ if(dp <= 0)
+ exit;
+ # t1 := daytime->now();
+ p2 := pi2(dp+1);
+ # t2 := daytime->now();
+ prpi(p2);
+ p1 := pi1(dp+1);
+ # t3 := daytime->now();
+ # sys->print("%d %d\n", t2-t1, t3-t2);
+ if(p1 == nil && p2 == nil)
+ fatal("too many dp: reduce dp or source base");
+ else if(p1 == nil)
+ p1 = p2;
+ else if(p2 == nil)
+ p2 = p1;
+ n1 := len p1;
+ n2 := len p2;
+ if(n1 != n2)
+ fatal(sys->sprint("lens differ %d %d", n1, n2));
+ f := array[10] of { * => 0 };
+ for(i := 0; i < n1; i++){
+ if(p1[i] != p2[i])
+ fatal(sys->sprint("arrays differ %d/%d: %d %d", i, n1, p1[i], p2[i]));
+ if(p1[i] < 0 || p1[i] >= BASE)
+ fatal(sys->sprint("bad array element %d: %d", i, p1[i]));
+ if(0){
+ p := p1[i];
+ for(j := 0; j < LBASE; j++){
+ f[p%10]++;
+ p /= 10;
+ }
+ }
+ }
+ # prpi(p1);
+ if(0){
+ t := 0;
+ for(i = 0; i < 10; i++){
+ sys->print("%d %d\n", i, f[i]);
+ t += f[i];
+ }
+ sys->print("T %d\n", t);
+ }
+}
+
+terms(dp: int, f: int, v: int): (int, int)
+{
+ p := dp;
+ t := 0;
+ for(;;){
+ t = 2 + int ((real p*log(real 10)+log(real v))/log(real f));
+ if(!(t&1))
+ t++;
+ e := int (log(real (v*(t+1)/2))/log(real 10))+1;
+ if(dp <= p-e)
+ break;
+ p += e;
+ }
+ # sys->fprint(stderr, "dp=%d p=%d f=%d v=%d terms=%d\n", dp, p, f, v, t);
+ if(t < f*f)
+ k := f*f;
+ else
+ k = t;
+ m := BASE*k;
+ if(m < 0 || m < BASE || m < k || m/BASE != k || m/k != BASE)
+ return (-1, -1);
+ return (t, p);
+}
+
+prpi(p: array of int)
+{
+ n := len p;
+ sys->print("π ≅ ");
+ m := BASE/10;
+ sys->print("%d.%.*d", p[0]/m, LBASE-1, p[0]%m);
+ for(i := 1; i < n; i++)
+ sys->print("%.*d", LBASE, p[i]);
+ sys->print("\n");
+}
+
+memcmp(b1: array of byte, b2: array of byte, n: int): (int, int, int)
+{
+ for(i := 0; i < n; i++)
+ if(b1[i] != b2[i])
+ return (i, int b1[i], int b2[i]);
+ return (-1, 0, 0);
+}
+
+picmp(f1: string, f2: string)
+{
+ fd1 := sys->open(f1, Sys->OREAD);
+ fd2 := sys->open(f2, Sys->OREAD);
+ if(fd1 == nil || fd2 == nil)
+ fatal(sys->sprint("cannot open %s or %s", f1, f2));
+ b1 := array[Sys->ATOMICIO] of byte;
+ b2 := array[Sys->ATOMICIO] of byte;
+ t := 0;
+ shouldexit := 0;
+ for(;;){
+ n1 := sys->read(fd1, b1, len b1);
+ n2 := sys->read(fd2, b2, len b2);
+ if(n1 <= 0 || n2 <= 0)
+ return;
+ if(shouldexit)
+ fatal("bad picmp");
+ if(n1 < n2)
+ (d, v1, v2) := memcmp(b1, b2, n1);
+ else
+ (d, v1, v2) = memcmp(b1, b2, n2);
+ if(d >= 0){
+ if(v1 == '\n' || v2 == '\n')
+ shouldexit = 1;
+ else
+ fatal(sys->sprint("%s %s differ at byte %d(%c %c)", f1, f2, t+d, v1, v2));
+ }
+ t += n1;
+ if(n1 != n2)
+ shouldexit = 1;
+ }
+}
+
+roundup(n: int, m: int): (int, int)
+{
+ r := m*((n+m-1)/m);
+ return (r, r/m);
+}
+
+pi1(dp: int): array of int
+{
+ fs := array[2] of { 5, 239 };
+ vs := array[2] of { 16, 4 };
+ ss := array[2] of { 1, -1 };
+ # sys->fprint(stderr, "π1\n");
+ return pi(dp, fs, vs, ss);
+}
+
+pi2(dp: int): array of int
+{
+ fs := array[3] of { 18, 57, 239 };
+ vs := array[3] of { 48, 32, 20 };
+ ss := array[3] of { 1, 1, -1 };
+ # sys->fprint(stderr, "π2\n");
+ return pi(dp, fs, vs, ss);
+}
+
+pi3(dp: int): array of int
+{
+ fs := array[4] of { 57, 239, 682, 12943 };
+ vs := array[4] of { 176, 28, 48, 96 };
+ ss := array[4] of { 1, 1, -1, 1 };
+ # sys->fprint(stderr, "π3\n");
+ return pi(dp, fs, vs, ss);
+}
+
+pi(dp: int, fs: array of int, vs: array of int, ss: array of int): array of int
+{
+ k := len fs;
+ n := cn := adp := 0;
+ (dp, n) = roundup(dp, LBASE);
+ cdp := dp;
+ m := array[k] of int;
+ for(i := 0; i < k; i++){
+ (m[i], adp) = terms(dp+1, fs[i], vs[i]);
+ if(m[i] < 0)
+ return nil;
+ if(adp > cdp)
+ cdp = adp;
+ }
+ (cdp, cn) = roundup(cdp, LBASE);
+ a := array[cn] of int;
+ p := array[cn] of int;
+ for(i = 0; i < cn; i++)
+ p[i] = 0;
+ for(i = 0; i < k; i++){
+ series(m[i], cn, fs[i], (vs[i]*BASE)/10, ss[i], a, p);
+ # sys->fprint(stderr, "term %d done\n", i+1);
+ }
+ return p[0: n];
+}
+
+series(m: int, n: int, f: int, v: int, s: int, a: array of int, p: array of int)
+{
+ i, j, k, q, r, r1, r2, n0: int;
+
+ v *= f;
+ f *= f;
+ for(j = 0; j < n; j++)
+ a[j] = 0;
+ a[0] = v;
+
+ if(s == 1)
+ series1(m, n, f, v, a, p);
+ else
+ series2(m, n, f, v, a, p);
+ return;
+
+ # following code now split
+ n0 = 0; # reaches n when very close to m so no check needed
+ for(i = 1; i <= m; i += 2){
+ r1 = r2 = 0;
+ for(j = n0; j < n; j++){
+ v = a[j]+r1;
+ q = v/f;
+ r1 = (v-q*f)*BASE;
+ a[j] = q;
+ v = q+r2;
+ q = v/i;
+ r2 = (v-q*i)*BASE;
+ for(k = j; q > 0; k--){
+ r = p[k]+s*q;
+ if(r >= BASE){
+ p[k] = r-BASE;
+ q = 1;
+ }
+ else if(r < 0){
+ p[k] = r+BASE;
+ q = 1;
+ }
+ else{
+ p[k] = r;
+ q = 0;
+ }
+ }
+ }
+ for(j = n0; j < n; j++){
+ if(a[j] == 0)
+ n0++;
+ else
+ break;
+ }
+ s = -s;
+ }
+}
+
+series1(m: int, n: int, f: int, v: int, a: array of int, p: array of int)
+{
+ i, j, k, q, r, r1, r2, n0: int;
+
+ n0 = 0;
+ for(i = 1; i <= m; i += 2){
+ r1 = r2 = 0;
+ for(j = n0; j < n; j++){
+ v = a[j]+r1;
+ q = v/f;
+ r1 = (v-q*f)*BASE;
+ a[j] = q;
+ v = q+r2;
+ q = v/i;
+ r2 = (v-q*i)*BASE;
+ for(k = j; q > 0; k--){
+ r = p[k]+q;
+ if(r >= BASE){
+ p[k] = r-BASE;
+ q = 1;
+ }
+ else{
+ p[k] = r;
+ q = 0;
+ }
+ }
+ }
+ for(j = n0; j < n; j++){
+ if(a[j] == 0)
+ n0++;
+ else
+ break;
+ }
+ i += 2;
+ r1 = r2 = 0;
+ for(j = n0; j < n; j++){
+ v = a[j]+r1;
+ q = v/f;
+ r1 = (v-q*f)*BASE;
+ a[j] = q;
+ v = q+r2;
+ q = v/i;
+ r2 = (v-q*i)*BASE;
+ for(k = j; q > 0; k--){
+ r = p[k]-q;
+ if(r < 0){
+ p[k] = r+BASE;
+ q = 1;
+ }
+ else{
+ p[k] = r;
+ q = 0;
+ }
+ }
+ }
+ for(j = n0; j < n; j++){
+ if(a[j] == 0)
+ n0++;
+ else
+ break;
+ }
+ }
+}
+
+series2(m: int, n: int, f: int, v: int, a: array of int, p: array of int)
+{
+ i, j, k, q, r, r1, r2, n0: int;
+
+ n0 = 0;
+ for(i = 1; i <= m; i += 2){
+ r1 = r2 = 0;
+ for(j = n0; j < n; j++){
+ v = a[j]+r1;
+ q = v/f;
+ r1 = (v-q*f)*BASE;
+ a[j] = q;
+ v = q+r2;
+ q = v/i;
+ r2 = (v-q*i)*BASE;
+ for(k = j; q > 0; k--){
+ r = p[k]-q;
+ if(r < 0){
+ p[k] = r+BASE;
+ q = 1;
+ }
+ else{
+ p[k] = r;
+ q = 0;
+ }
+ }
+ }
+ for(j = n0; j < n; j++){
+ if(a[j] == 0)
+ n0++;
+ else
+ break;
+ }
+ i += 2;
+ r1 = r2 = 0;
+ for(j = n0; j < n; j++){
+ v = a[j]+r1;
+ q = v/f;
+ r1 = (v-q*f)*BASE;
+ a[j] = q;
+ v = q+r2;
+ q = v/i;
+ r2 = (v-q*i)*BASE;
+ for(k = j; q > 0; k--){
+ r = p[k]+q;
+ if(r >= BASE){
+ p[k] = r-BASE;
+ q = 1;
+ }
+ else{
+ p[k] = r;
+ q = 0;
+ }
+ }
+ }
+ for(j = n0; j < n; j++){
+ if(a[j] == 0)
+ n0++;
+ else
+ break;
+ }
+ }
+}
+
+fatal(e: string)
+{
+ sys->print("%s\n", e);
+ exit;
+}
diff --git a/appl/math/polyfill.b b/appl/math/polyfill.b
new file mode 100644
index 00000000..c799d3db
--- /dev/null
+++ b/appl/math/polyfill.b
@@ -0,0 +1,440 @@
+implement Polyfill;
+
+include "sys.m";
+ sys: Sys;
+include "draw.m";
+ draw: Draw;
+ Point, Rect, Image, Endsquare: import draw;
+include "math/polyfill.m";
+
+∞: con 16r7fffffff;
+
+init()
+{
+ sys = load Sys Sys->PATH;
+ draw = load Draw Draw->PATH;
+}
+
+initzbuf(r: Rect): ref Zstate
+{
+ if(sys == nil)
+ init();
+ s := ref Zstate;
+ s.r = r;
+ s.xlen = r.dx();
+ s.ylen = r.dy();
+ s.xylen = s.xlen*s.ylen;
+ s.zbuf0 = array[s.xylen] of int;
+ s.zbuf1 = array[s.xylen] of int;
+ return s;
+}
+
+clearzbuf(s: ref Zstate)
+{
+ b0 := s.zbuf0;
+ b1 := s.zbuf1;
+ n := s.xylen;
+ for(i := 0; i < n; i++)
+ b0[i] = b1[i] = ∞;
+}
+
+setzbuf(s: ref Zstate, zd: int)
+{
+ b0 := s.zbuf0;
+ b1 := s.zbuf1;
+ n := s.xylen;
+ for(i := 0; i < n; i++)
+ b0[i] = b1[i] = zd;
+}
+
+Seg: adt
+{
+ p0: Point;
+ p1: Point;
+ num: int;
+ den: int;
+ dz: int;
+ dzrem: int;
+ z: int;
+ zerr: int;
+ d: int;
+};
+
+fillline(dst: ref Image, left: int, right: int, y: int, src: ref Image, p: Point)
+{
+ p.x += left;
+ p.y += y;
+ dst.line((left, y), (right, y), Endsquare, Endsquare, 0, src, p);
+}
+
+filllinez(dst: ref Image, left: int, right: int, y: int, z: int, e: int, dx: int, k: int, zbuf0: array of int, zbuf1: array of int, src: ref Image, p: Point)
+{
+ prevx := ∞;
+ for(x := left; x <= right; x++){
+ if(z+e < zbuf0[k] || (z-e <= zbuf1[k] && x != right && prevx != ∞)){
+ zbuf0[k] = z-e;
+ zbuf1[k] = z+e;
+ if(prevx == ∞)
+ prevx = x;
+ }
+ else if(prevx != ∞){
+ fillline(dst, prevx, x-1, y, src, p);
+ prevx = ∞;
+ }
+ z += dx;
+ k++;
+ }
+ if(prevx != ∞)
+ fillline(dst, prevx, right, y, src, p);
+}
+
+fillpoly(dst: ref Image, vert: array of Point, w: int, src: ref Image, sp: Point, zstate: ref Zstate, dc: int, dx: int, dy: int)
+{
+ p0: Point;
+ i: int;
+
+ nvert := len vert;
+ if(nvert == 0)
+ return;
+ fixshift := 0;
+ seg := array[nvert+2] of ref Seg;
+ if(seg == nil)
+ return;
+ segtab := array[nvert+1] of ref Seg;
+ if(segtab == nil)
+ return;
+
+ sp.x = (sp.x - vert[0].x) >> fixshift;
+ sp.y = (sp.y - vert[0].y) >> fixshift;
+ p0 = vert[nvert-1];
+ if(!fixshift) {
+ p0.x <<= 1;
+ p0.y <<= 1;
+ }
+ for(i = 0; i < nvert; i++) {
+ segtab[i] = ref Seg;
+ segtab[i].p0 = p0;
+ p0 = vert[i];
+ if(!fixshift) {
+ p0.x <<= 1;
+ p0.y <<= 1;
+ }
+ segtab[i].p1 = p0;
+ segtab[i].d = 1;
+ }
+ if(!fixshift)
+ fixshift = 1;
+
+ xscan(dst, seg, segtab, nvert, w, src, sp, zstate, dc, dx, dy, fixshift);
+}
+
+mod(x: int, y: int): int
+{
+ z: int;
+
+ z = x%y;
+ if((z^y) > 0 || z == 0)
+ return z;
+ return z + y;
+}
+
+sdiv(x: int, y: int): int
+{
+ if((x^y) >= 0 || x == 0)
+ return x/y;
+ return (x+((y>>30)|1))/y-1;
+}
+
+smuldivmod(x: int, y: int, z: int): (int, int)
+{
+ mod: int;
+ vx: int;
+
+ if(x == 0 || y == 0)
+ return (0, 0);
+ vx = x;
+ vx *= y;
+ mod = vx % z;
+ if(mod < 0)
+ mod += z;
+ if((vx < 0) == (z < 0))
+ return (vx/z, mod);
+ return (-((-vx)/z), mod);
+}
+
+xscan(dst: ref Image, seg: array of ref Seg, segtab: array of ref Seg, nseg: int, wind: int, src: ref Image, spt: Point, zstate: ref Zstate, dc: int, dx: int, dy: int, fixshift: int)
+{
+ y, maxy, x, x2, onehalf: int;
+ ep, next, p, q, s: int;
+ n, i, iy, cnt, ix, ix2, minx, maxx, zinc, k, zv: int;
+ pt: Point;
+ sp: ref Seg;
+
+ er := (abs(dx)+abs(dy)+1)/2;
+ zr := zstate.r;
+ xlen := zstate.xlen;
+ zbuf0 := zstate.zbuf0;
+ zbuf1 := zstate.zbuf1;
+ s = 0;
+ p = 0;
+ for(i=0; i<nseg; i++) {
+ sp = seg[p] = segtab[s];
+ if(sp.p0.y == sp.p1.y){
+ s++;
+ continue;
+ }
+ if(sp.p0.y > sp.p1.y) {
+ pt = sp.p0;
+ sp.p0 = sp.p1;
+ sp.p1 = pt;
+ sp.d = -sp.d;
+ }
+ sp.num = sp.p1.x - sp.p0.x;
+ sp.den = sp.p1.y - sp.p0.y;
+ sp.dz = sdiv(sp.num, sp.den) << fixshift;
+ sp.dzrem = mod(sp.num, sp.den) << fixshift;
+ sp.dz += sdiv(sp.dzrem, sp.den);
+ sp.dzrem = mod(sp.dzrem, sp.den);
+ p++;
+ s++;
+ }
+ n = p;
+ if(n == 0)
+ return;
+ seg[p] = nil;
+ qsortycompare(seg, p);
+
+ onehalf = 0;
+ if(fixshift)
+ onehalf = 1 << (fixshift-1);
+
+ minx = dst.clipr.min.x;
+ maxx = dst.clipr.max.x;
+
+ y = seg[0].p0.y;
+ if(y < (dst.clipr.min.y << fixshift))
+ y = dst.clipr.min.y << fixshift;
+ iy = (y + onehalf) >> fixshift;
+ y = (iy << fixshift) + onehalf;
+ maxy = dst.clipr.max.y << fixshift;
+ k = (iy-zr.min.y)*xlen;
+ zv = dc+iy*dy;
+
+ ep = next = 0;
+
+ while(y<maxy) {
+ for(q = p = 0; p < ep; p++) {
+ sp = seg[p];
+ if(sp.p1.y < y)
+ continue;
+ sp.z += sp.dz;
+ sp.zerr += sp.dzrem;
+ if(sp.zerr >= sp.den) {
+ sp.z++;
+ sp.zerr -= sp.den;
+ if(sp.zerr < 0 || sp.zerr >= sp.den)
+ sys->print("bad ratzerr1: %d den %d dzrem %d\n", sp.zerr, sp.den, sp.dzrem);
+ }
+ seg[q] = sp;
+ q++;
+ }
+
+ for(p = next; seg[p] != nil; p++) {
+ sp = seg[p];
+ if(sp.p0.y >= y)
+ break;
+ if(sp.p1.y < y)
+ continue;
+ sp.z = sp.p0.x;
+ (zinc, sp.zerr) = smuldivmod(y - sp.p0.y, sp.num, sp.den);
+ sp.z += zinc;
+ if(sp.zerr < 0 || sp.zerr >= sp.den)
+ sys->print("bad ratzerr2: %d den %d ratdzrem %d\n", sp.zerr, sp.den, sp.dzrem);
+ seg[q] = sp;
+ q++;
+ }
+ ep = q;
+ next = p;
+
+ if(ep == 0) {
+ if(seg[next] == nil)
+ break;
+ iy = (seg[next].p0.y + onehalf) >> fixshift;
+ y = (iy << fixshift) + onehalf;
+ k = (iy-zr.min.y)*xlen;
+ zv = dc+iy*dy;
+ continue;
+ }
+
+ zsort(seg, ep);
+
+ for(p = 0; p < ep; p++) {
+ sp = seg[p];
+ cnt = 0;
+ x = sp.z;
+ ix = (x + onehalf) >> fixshift;
+ if(ix >= maxx)
+ break;
+ if(ix < minx)
+ ix = minx;
+ cnt += sp.d;
+ p++;
+ sp = seg[p];
+ for(;;) {
+ if(p == ep) {
+ sys->print("xscan: fill to infinity");
+ return;
+ }
+ cnt += sp.d;
+ if((cnt&wind) == 0)
+ break;
+ p++;
+ sp = seg[p];
+ }
+ x2 = sp.z;
+ ix2 = (x2 + onehalf) >> fixshift;
+ if(ix2 <= minx)
+ continue;
+ if(ix2 > maxx)
+ ix2 = maxx;
+ filllinez(dst, ix, ix2, iy, zv+ix*dx, er, dx, k+ix-zr.min.x, zbuf0, zbuf1, src, spt);
+ }
+ y += (1<<fixshift);
+ iy++;
+ k += xlen;
+ zv += dy;
+ }
+}
+
+zsort(seg: array of ref Seg, ep: int)
+{
+ done: int;
+ s: ref Seg;
+ q, p: int;
+
+ if(ep < 20) {
+ # bubble sort by z - they should be almost sorted already
+ q = ep;
+ do {
+ done = 1;
+ q--;
+ for(p = 0; p < q; p++) {
+ if(seg[p].z > seg[p+1].z) {
+ s = seg[p];
+ seg[p] = seg[p+1];
+ seg[p+1] = s;
+ done = 0;
+ }
+ }
+ } while(!done);
+ } else {
+ q = ep-1;
+ for(p = 0; p < q; p++) {
+ if(seg[p].z > seg[p+1].z) {
+ qsortzcompare(seg, ep);
+ break;
+ }
+ }
+ }
+}
+
+ycompare(s0: ref Seg, s1: ref Seg): int
+{
+ y0, y1: int;
+
+ y0 = s0.p0.y;
+ y1 = s1.p0.y;
+
+ if(y0 < y1)
+ return -1;
+ if(y0 == y1)
+ return 0;
+ return 1;
+}
+
+zcompare(s0: ref Seg, s1: ref Seg): int
+{
+ z0, z1: int;
+
+ z0 = s0.z;
+ z1 = s1.z;
+
+ if(z0 < z1)
+ return -1;
+ if(z0 == z1)
+ return 0;
+ return 1;
+}
+
+qsortycompare(a : array of ref Seg, n : int)
+{
+ i, j : int;
+ t : ref Seg;
+
+ while(n > 1) {
+ i = n>>1;
+ t = a[0]; a[0] = a[i]; a[i] = t;
+ i = 0;
+ j = n;
+ for(;;) {
+ do
+ i++;
+ while(i < n && ycompare(a[i], a[0]) < 0);
+ do
+ j--;
+ while(j > 0 && ycompare(a[j], a[0]) > 0);
+ if(j < i)
+ break;
+ t = a[i]; a[i] = a[j]; a[j] = t;
+ }
+ t = a[0]; a[0] = a[j]; a[j] = t;
+ n = n-j-1;
+ if(j >= n) {
+ qsortycompare(a, j);
+ a = a[j+1:];
+ } else {
+ qsortycompare(a[j+1:], n);
+ n = j;
+ }
+ }
+}
+
+qsortzcompare(a : array of ref Seg, n : int)
+{
+ i, j : int;
+ t : ref Seg;
+
+ while(n > 1) {
+ i = n>>1;
+ t = a[0]; a[0] = a[i]; a[i] = t;
+ i = 0;
+ j = n;
+ for(;;) {
+ do
+ i++;
+ while(i < n && zcompare(a[i], a[0]) < 0);
+ do
+ j--;
+ while(j > 0 && zcompare(a[j], a[0]) > 0);
+ if(j < i)
+ break;
+ t = a[i]; a[i] = a[j]; a[j] = t;
+ }
+ t = a[0]; a[0] = a[j]; a[j] = t;
+ n = n-j-1;
+ if(j >= n) {
+ qsortzcompare(a, j);
+ a = a[j+1:];
+ } else {
+ qsortzcompare(a[j+1:], n);
+ n = j;
+ }
+ }
+}
+
+abs(n: int): int
+{
+ if(n < 0)
+ return -n;
+ return n;
+}
diff --git a/appl/math/polyhedra.b b/appl/math/polyhedra.b
new file mode 100644
index 00000000..e7a51d6c
--- /dev/null
+++ b/appl/math/polyhedra.b
@@ -0,0 +1,195 @@
+implement Polyhedra;
+
+include "sys.m";
+ sys: Sys;
+include "bufio.m";
+ bufio: Bufio;
+ Iobuf: import bufio;
+include "math/polyhedra.m";
+
+scanpolyhedra(f: string): (int, ref Polyhedron, ref Iobuf)
+{
+ first, last: ref Polyhedron;
+ D: int;
+
+ if(sys == nil)
+ sys = load Sys Sys->PATH;
+ if(bufio == nil)
+ bufio = load Bufio Bufio->PATH;
+ b := bufio->open(f, Sys->OREAD);
+ if(b == nil)
+ return (0, nil, nil);
+ n := 0;
+ for(;;){
+ s := getstring(b);
+ if(s == nil)
+ break;
+ n++;
+ p := ref Polyhedron;
+ if(first == nil)
+ first = p;
+ else{
+ last.nxt = p;
+ p.prv = last;
+ }
+ last = p;
+ p.name = s;
+ p.dname = getstring(b);
+ b.gets('\n');
+ (p.allf, p.adj) = scanvc(getstring(b));
+ b.gets('\n');
+ b.gets('\n');
+ b.gets('\n');
+ l := getstring(b);
+ (p.indx, l) = getint(l);
+ (p.V, l) = getint(l);
+ (p.E, l) = getint(l);
+ (p.F, l) = getint(l);
+ (nil, l) = getint(l);
+ (D, l) = getint(l);
+ (p.anti, l) = getint(l);
+ p.concave = D != 1 || p.allf;
+ p.offset = b.offset();
+ tot := 2*p.V+2*p.F;
+ for(i := 0; i < tot; i++)
+ b.gets('\n');
+ if(p.indx < 58 || p.indx == 59 || p.indx == 66 || p.indx == 67)
+ p.inc = 0.1;
+ else
+ p.inc = 0.0;
+ # sys->print("%d: %d %d %d %d %s\n", p.indx, p.allf, D != 1, p.anti, p.concave, vc);
+ }
+ first.prv = last;
+ last.nxt = first;
+ return (n, first, b);
+}
+
+getpolyhedra(p: ref Polyhedron, b: ref Iobuf)
+{
+ q := p;
+ do{
+ getpolyhedron(q, b);
+ q = q.nxt;
+ }while(q != p);
+}
+
+getpolyhedron(p: ref Polyhedron, b: ref Iobuf)
+{
+ if(p.v != nil)
+ return;
+ b.seek(p.offset, Bufio->SEEKSTART);
+ p.v = array[p.V] of Vector;
+ for(i := 0; i < p.V; i++)
+ p.v[i] = getvector(b);
+ p.f = array[p.F] of Vector;
+ for(i = 0; i < p.F; i++)
+ p.f[i] = getvector(b);
+ p.fv = array[p.F] of array of int;
+ for(i = 0; i < p.F; i++)
+ p.fv[i] = getarray(b, p.adj);
+ p.vf = array[p.V] of array of int;
+ for(i = 0; i < p.V; i++)
+ p.vf[i] = getarray(b, p.adj);
+}
+
+getstring(b: ref Iobuf): string
+{
+ s := b.gets('\n');
+ if(s == nil)
+ return nil;
+ if(s[0] == '#')
+ return getstring(b);
+ if(s[len s - 1] == '\n')
+ return s[0: len s - 1];
+ return s;
+}
+
+getvector(b: ref Iobuf): Vector
+{
+ v: Vector;
+
+ s := getstring(b);
+ (v.x, s) = getreal(s);
+ (v.y, s) = getreal(s);
+ (v.z, s) = getreal(s);
+ return v;
+}
+
+getarray(b: ref Iobuf, adj: int): array of int
+{
+ n, d: int;
+
+ s := getstring(b);
+ (n, s) = getint(s);
+ a := array[n+2] of int;
+ a[0] = n;
+ for(i := 1; i <= n; i++)
+ (a[i], s) = getint(s);
+ (d, s) = getint(s);
+ if(d == 0 || d == n-1 || adj)
+ d = 1;
+ a[n+1] = d;
+ return a;
+}
+
+getint(s: string): (int, string)
+{
+ n := int s;
+ for(i := 0; i < len s && s[i] == ' '; i++)
+ ;
+ for( ; i < len s; i++)
+ if(s[i] == ' ')
+ return (n, s[i+1:]);
+ return (n, nil);
+}
+
+getreal(s: string): (real, string)
+{
+ r := real s;
+ for(i := 0; i < len s && s[i] == ' '; i++)
+ ;
+ for( ; i < len s; i++)
+ if(s[i] == ' ')
+ return (r, s[i+1:]);
+ return (r, nil);
+}
+
+vftab := array[] of { 0, 0, 0, 2, 3, 3, 5, 0, 3, 0, 3 };
+
+scanvc(s: string): (int, int)
+{
+ af := 0;
+ ad := 0;
+ fd := ld := 1;
+ ln := len s;
+ if(ln > 0 && s[0] == '('){
+ s = s[1:];
+ ln--;
+ }
+ while(ln > 0 && s[ln-1] != ')'){
+ s = s[0: ln-1];
+ ln--;
+ }
+ (m, lst) := sys->tokenize(s, ".");
+ for(l := lst ; l != nil; l = tl l){
+ (m, lst) = sys->tokenize(hd l, "/");
+ if(m == 1)
+ (n, d) := (int hd lst, 1);
+ else if(m == 2)
+ (n, d) = (int hd lst, int hd tl lst);
+ else
+ sys->print("vc error\n");
+ if(d != 1 && d == vftab[n])
+ af = 1;
+ if(d == n-1)
+ d = 1;
+ if(l == lst)
+ fd = d;
+ else if(ld != 1 && d != 1)
+ ad = 1;
+ ld = d;
+ }
+ if(ld != 1 && fd != 1)
+ ad = 1;
+ return (af, ad);
+}
diff --git a/appl/math/powers.b b/appl/math/powers.b
new file mode 100644
index 00000000..d9c5259c
--- /dev/null
+++ b/appl/math/powers.b
@@ -0,0 +1,672 @@
+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;
+ }
+}
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";
+}
+
diff --git a/appl/math/sieve.b b/appl/math/sieve.b
new file mode 100644
index 00000000..25d4b0f1
--- /dev/null
+++ b/appl/math/sieve.b
@@ -0,0 +1,220 @@
+implement Sieve;
+
+include "sys.m";
+ sys: Sys;
+include "draw.m";
+include "arg.m";
+ arg: Arg;
+
+M: con 16*1024*1024;
+N: con 8*M;
+T: con 2*1024*1024;
+
+limit := array[5] of { M, N, 2*N, 3*N, 15*(N/4) };
+
+Sieve: module
+{
+ init: fn(nil: ref Draw->Context, argv: list of string);
+};
+
+init(nil: ref Draw->Context, argv: list of string)
+{
+ sys = load Sys Sys->PATH;
+ arg = load Arg Arg->PATH;
+
+ np := 0;
+ alg := 3;
+ arg->init(argv);
+ while((c := arg->opt()) != 0){
+ case (c){
+ 'a' =>
+ alg = int arg->arg();
+ }
+ }
+ if(alg < 0 || alg > 4)
+ alg = 3;
+ lim := limit[alg];
+ argv = arg->argv();
+ if(argv != nil)
+ lim = int hd argv;
+ if(lim < 0 || lim > limit[alg])
+ lim = limit[alg];
+ if(lim < 6){
+ if(lim > 2){
+ sys->print("2\n");
+ np++;
+ }
+ if(lim > 3){
+ sys->print("3\n");
+ np++;
+ }
+ }
+ else{
+ case (alg){
+ 0 => np = init0(lim);
+ 1 => np = init1(lim);
+ 2 => np = init2(lim);
+ 3 => np = init3(lim);
+ 4 => np = init4(lim);
+ }
+ }
+ sys->print("%d primes < %d\n", np, lim);
+}
+
+init0(lim: int): int
+{
+ p := array[lim] of byte;
+ for(i := 0; i < lim; i++)
+ p[i] = byte 1;
+ p[0] = p[1] = byte 0;
+ np := 0;
+ for(i = 0; i < lim; i++){
+ if(p[i] == byte 1){
+ np++;
+ sys->print("%d\n", i);
+ for(j := i+i; j < lim; j += i)
+ p[j] = byte 0;
+ }
+ }
+ return np;
+}
+
+init1(lim: int): int
+{
+ n := (lim+31)/32;
+ p := array[n] of int;
+ for(i := 0; i < n; i++)
+ p[i] = int 16rffffffff;
+ p[0] = int 16rfffffffc;
+ np := 0;
+ for(i = 0; i < lim; i++){
+ if(p[i>>5] & (1<<(i&31))){
+ np++;
+ sys->print("%d\n", i);
+ for(j := i+i; j < lim; j += i)
+ p[j>>5] &= ~(1<<(j&31));
+ }
+ }
+ return np;
+}
+
+init2(lim: int): int
+{
+ n := ((lim+1)/2+31)/32;
+ p := array[n] of int;
+ for(i := 0; i < n; i++)
+ p[i] = int 16rffffffff;
+ p[0] = int 16rfffffffe;
+ np := 1;
+ sys->print("%d\n", 2);
+ for(i = 1; i < lim; i += 2){
+ k := (i-1)>>1;
+ if(p[k>>5] & (1<<(k&31))){
+ np++;
+ sys->print("%d\n", i);
+ inc := i+i;
+ for(j := i+i+i; j < lim; j += inc){
+ k = (j-1)>>1;
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ }
+ }
+ return np;
+}
+
+init3(lim: int): int
+{
+ n := ((lim+2)/3+31)/32;
+ p := array[n] of int;
+ for(i := 0; i < n; i++)
+ p[i] = int 16rffffffff;
+ p[0] = int 16rfffffffe;
+ np := 2;
+ sys->print("%d\n", 2);
+ sys->print("%d\n", 3);
+ d := 2;
+ for(i = 1; i < lim; i += d){
+ k := (i-1)/3;
+ if(p[k>>5] & (1<<(k&31))){
+ np++;
+ sys->print("%d\n", i);
+ inc := 6*i;
+ for(j := 5*i; j > 0 && j < lim; j += inc){
+ k = (j-1)/3;
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ for(j = 7*i; j > 0 && j < lim; j += inc){
+ k = (j-1)/3;
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ }
+ d = 6-d;
+ }
+ return np;
+}
+
+init4(lim: int): int
+{
+ n := (4*((lim+14)/15)+31)/32;
+ p := array[n] of int;
+ for(i := 0; i < n; i++)
+ p[i] = int 16rffffffff;
+ p[0] = int 16rfffffffe;
+ np := 3;
+ sys->print("%d\n", 2);
+ sys->print("%d\n", 3);
+ sys->print("%d\n", 5);
+ m := -1;
+ d := array[8] of { 6, 4, 2, 4, 2, 4, 6, 2 };
+ for(i = 1; i < lim; i += d[m]){
+ k := (17*(i%30-1))/60+8*(i/30);
+ if(p[k>>5] & (1<<(k&31))){
+ np++;
+ sys->print("%d\n", i);
+ inc := 30*i;
+ for(j := 7*i; j > 0 && j < lim; j += inc){
+ k = (17*(j%30-1))/60+8*(j/30);
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ for(j = 11*i; j > 0 && j < lim; j += inc){
+ k = (17*(j%30-1))/60+8*(j/30);
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ for(j = 13*i; j > 0 && j < lim; j += inc){
+ k = (17*(j%30-1))/60+8*(j/30);
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ for(j = 17*i; j > 0 && j < lim; j += inc){
+ k = (17*(j%30-1))/60+8*(j/30);
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ for(j = 19*i; j > 0 && j < lim; j += inc){
+ k = (17*(j%30-1))/60+8*(j/30);
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ for(j = 23*i; j > 0 && j < lim; j += inc){
+ k = (17*(j%30-1))/60+8*(j/30);
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ for(j = 29*i; j > 0 && j < lim; j += inc){
+ k = (17*(j%30-1))/60+8*(j/30);
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ for(j = 31*i; j > 0 && j < lim; j += inc){
+ k = (17*(j%30-1))/60+8*(j/30);
+ p[k>>5] &= ~(1<<(k&31));
+ }
+ }
+ m++;
+ if(m == 8)
+ m = 0;
+ }
+ return np;
+}
+
+init5(lim: int): int
+{
+ # you must be joking
+ lim = 0;
+ return 0;
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-12-10 22:40:29 +0000