use Array: all; use SimplePrint: all; use SimpleFibre: all; use Math: all; use String:{tochar}; import Char:all; use Bits:all; /* % This is the APEX stdlib.sis include file. % Standard equates and constants for APL compiler % Also standard coercion functions */ #define toB(x) tob((x)) #define toI(x) toi((x)) #define toD(x) tod((x)) #define toC(x) (x) #define toc(x) ((x)) void APEXERROR(char[.] msg) { /* Error function. This needs work, as it should kill * the running task, too. */ /*print(msg); */ return(); } /* Structural function utility functions */ /* Ravel utility */ #define APEXRavel(y) (reshape([prod(shape(y))],y)) inline int VectorRotateAmount(int x, int y) { /* Normalize x rotate for array of shape y on selected axis */ /* normalize rotation count */ if (x>0) z = _mod_(x,y); else z = y - _mod_(_abs_(x),y); return(z); } #define APEXRESHAPE(TYPE) \ inline TYPE[*] APEXReshape(int[.] x, TYPE[*] y) \ { /* APEX vector x reshape, with item reuse */ \ ry = APEXRavel(y); \ zxrho = prod(x); /* THIS NEEDS XRHO FOR CODE SAFETY!! */ \ yxrho = shape(ry)[0]; \ if( zxrho <= yxrho) { /* No element resuse case */ \ z = take([zxrho],ry); \ } else { \ ncopies = zxrho/yxrho; /* # complete copies of y. */ \ /* FIXME: y empty case !*/ \ z = with(. <= [i] <= .) \ genarray( [ncopies], y,y); \ /* Now append the leftover bits */ \ z = APEXRavel(z) ++ take([zxrho-(ncopies*yxrho)],ry); \ } \ return(reshape(x,z)); \ } APEXRESHAPE(bool) APEXRESHAPE(int) APEXRESHAPE(double) APEXRESHAPE(char) inline bool SameShape(int[*] x, int[*] y) { /* Predicate for two shape vectors having same shape */ if (dim(x) != dim(y)) z = false; else z = with ([0] <= i < [dim(y)]) fold(andBBB, true, (shape(x))[i] == (shape(y))[i]); return(z); } #define APEXmiota(N) with (. <= [iv] <= .) genarray([N], iv) #define APEXPI 3.1415926535897932d #define APEXE 2.718281828d /* APEXPINFINITYI largest integer */ #define APEXPINFINITYD 1.7976931348623156D308 /* APEXMINFINITYI smallest integer */ #define APEXMINFINITYD -1.7976931348623156D308 /* % Floating-point utilities % This taken from page 93 of the 1993-01-06 version % of Committee Draft 1 of the Extended ISO APL Standard % NO {quad}ct support yet! */ inline int Dsignum(double y) { /* signum double */ if (0.0 == y) z = 0; else if (0.0 > y) z = -1; else z = 1; return(z); } inline int Isignum(int y) { /* signum int */ if (0 == y) z = 0; else if (0 > y) z = -1; else z = 1; return(z); } inline double Dresidue(double x, double y, double QUADct) { /* Double residue double */ /* See Iresidue for definition */ if (0.0 == x) nx = 1.0; else nx = x; z = y - x * tod(Dfloor(y/nx, QUADct)); return(z); } inline int Iresidue(int x, int y) { /* Integer residue integer */ /* This definition is taken from SHARP APL Refman May 1991, p.6-26. * It extends the definition of residue to fractional right arguments * and to zero, negative and fractional left arguments. * r= y-x times floor y divide x+0=x */ if (0 == x) nx = 1; else nx = x; z = y - x * (y/nx); return(z); } inline int DfloorNoFuzz(double y) { /* Exact floor (no fuzz) */ return(toi(floor(y))); } inline int Dfloor(double y, double QUADct) { /* Fuzzy floor */ /* Definition taken from SHARP APL Refman May 1991, p.6-23 * floor: n <- (signum y) times nofuzzfloor 0.5+abs y) * z <- n-(QUADct times 1 max abs y)<(n-y) * * If you want a double result, write: "y - 1| y". * */ n = tod(DfloorNoFuzz(0.5+fabs(y))); if (y < 0.0) n = -n; else if (0.0 == y) n = 0.0; range = fabs(y); if (1.0 > range) range = 1.0; fuzzlim = QUADct*range; ny = n-y; if (fuzzlim < ny) z = n - 1.0; else z = n; return(toi(z)); } /* % 1996-12-07 Pearl Harbor Day. Try to fix bug in DBank infdivm % benchmark. (4 5 rho 6) + rank 1 (4 1 rho 7) introduces singleton % vectors into scalar fn calls. % Hence, we need support for singletons within scalar fns. % There is an similar, but independent, bug in rank support for % the extension case. */ /* vector-scalar simple search loops * This is origin-0 x1 iota y0 * cfn is comparator function type suffix to use, e.g, b,i,d,c, */ #define FindFirst(x,y,cfn,QUADct) \ sx = shape(x)[0]; \ z = sx; /* if not found */ \ for(i=0; iy[i]; i=shp; } return(z); } inline bool comparatorfnb(bool x, bool y,double QUADct) { /* Boolean comparator */ z = (x == y); return(z); } inline bool comparatorfni(int x, int y, double QUADct) { /* integer comparator */ z = (x == y); return(z); } inline bool comparatorfnd(double x, double y, double QUADct) { /* double comparator with fuzz */ L = _abs_(x-y); /* Tolerant equality */ R = QUADct*max(_abs_(x),_abs_(y)); z = (L <= R); return(z); } inline bool comparatorfnc(char x, char y,double QUADct) { /* char comparator */ z = (x == y); return(z); } inline bool comparatorfndnf(double x,double y, double QUADct) { /* double comparator, no fuzz */ z = (x == y); return(z); } /* Lehmer's random number generator */ #define lehmer(qrl) mod(qrl*16807,2147483647) /* Transposes */ #define APEXTRANSPOSE(TYPE,OTFILL) \ inline TYPE[*] APEXTranspose(TYPE[*] y) \ { \ z = with(iv) \ ( . <= iv <= .) : y[reverse( iv)]; \ genarray( reverse( shape(y)), OTFILL); \ return(z); \ } APEXTRANSPOSE(bool,false) APEXTRANSPOSE(int,0) APEXTRANSPOSE(double,0.0) APEXTRANSPOSE(char,' ') /* End of transposes */ /* End of boilerplate */ /* Monadic Scalar Function kernels on scalars */ inline int,int queryXII(int y, int QUADrl) {return(42,QUADrl+1); } inline bool modXBB(bool y) { return(y); } inline int modXII(int y) { return(abs(y)); } inline double modXDD(double y) { return(abs(y)); } inline bool barXBB(bool y) { return(!y); } inline int barXII(int y) { return(-y); } inline double barXDD(double y) { return(-y); } inline bool plusXBB(bool y) { return(y); } inline int plusXII(int y) { return(y); } inline double plusXDD(double y) { return(y); } inline bool minXBB(bool y) { return(y); } inline int minXII(int y) { return(y); } inline int minXDI(double y,double QUADct) { return(Dfloor(y,QUADct)); } inline bool maxXBB(bool y) { return(y); } inline int maxXII(int y) { return(y); } inline int maxXDI(double y, double QUADct) { return(-Dfloor(-y, QUADct)); } inline bool mpyXBB(bool y) { return(y); } inline int mpyXII(int y) { return(Isignum(y)); } inline int mpyXDD(double y) { return(Dsignum(y)); } inline bool notXBB(bool y) { return(!y); } inline bool notXII(int y) { return(!tob(y)); } inline bool notXDD(double y) { return(!tob(y)); } inline double divXBB(bool y) { return(1.0/tod(y)); } inline double divXII(int y) { return(1.0/tod(y)); } inline double divXDD(double y) { return(1.0/y); } inline double starXBD(bool y) { return(pow(APEXE,tod(y))); } inline double starXID(int y) { return(pow(APEXE,tod(y))); } inline double starXDD(double y) { return(pow(APEXE,y)); } inline double logXBB(bool y) { return(log(tod(y))); } inline double logXII(int y) { return(log(tod(y))); } inline double logXDD(double y) { return(log(y)); } inline double circXBD(bool y) { return(APEXPI*tod(y)); } inline double circXID(int y) { return(APEXPI*tod(y)); } inline double circXDD(double y) { return(APEXPI*y); } /* Dyadic Scalar Function kernel macro definitions */ /* x plus y */ inline int plusBBI(bool x, bool y) { return(toi(x)+toi(y)); } inline int plusBII(bool x, int y) { return(toi(x)+y); } inline int plusIBI(int x, bool y) { return(x+toi(y)); } inline int plusIII(int x, int y) { return(x+y); } inline double plusDDD(double x, double y) { return(x+y); } inline double plusBDD(bool x, double y) { return(tod(x)+y); } inline double plusDBD(double x, bool y) { return(x+tod(y)); } inline double plusIDD(int x, double y) { return(tod(x)+y); } inline double plusDID(double x, int y) { return(x+tod(y)); } /* x minus y */ inline int barBBI(bool x, bool y) { return(toi(x)-toi(y)); } inline int barBII(bool x, int y) { return(toi(x)-y); } inline int barIBI(int x, bool y) { return(x-toi(y)); } inline int barIII(int x, int y) { return(x-y); } inline double barDDD(double x, double y) { return(x-y); } inline double barBDD(bool x, double y) { return(tod(x)-y); } inline double barDBD(double x, bool y) { return(x-tod(y)); } inline double barIDD(int x, double y) { return(tod(x)-y); } inline double barDID(double x, int y) { return(x-tod(y)); } /* x times y */ inline bool mpyBBB(bool x, bool y) { return(x & y); } inline int mpyBII(bool x, int y) { return(toi(x)*y); } inline int mpyIBI(int x, bool y) { return(x*toi(y)); } inline int mpyIII(int x, int y) { return(x*y); } inline double mpyDDD(double x, double y) { return(x*y); } inline double mpyBDD(bool x, double y) { return(tod(x)*y); } inline double mpyDBD(double x, bool y) { return(x*tod(y)); } inline double mpyIDD(int x, double y) { return(tod(x)*y); } inline double mpyDID(double x, int y) { return(x*tod(y)); } /* x divided by y */ inline double divDDD(double x, double y) { if (x == y) z = 1.0d; else z = x/y; return(z); } inline double divBDD(bool x, double y) { if (tod(x) == y) z = 1.0d; else z = tod(x)/y; return(z); } inline double divIDD(int x, double y) { if (tod(x) == y) z = 1.0d; else z = tod(x)/y; return(z); } inline double divBID(bool x, int y) { if (tod(x) == tod(y)) z = 1.0d; else z = tod(x)/tod(y); return(z); } inline double divIID(int x, int y) { if (tod(x) == tod(y)) z = 1.0d; else z = tod(x)/tod(y); return(z); } inline double divDID(double x, int y) { if (x == tod(y)) z = 1.0d; else z = x/tod(y); return(z); } inline double divBBD(bool x, bool y) { if (x == y) z = 1.0d; else z = tod(x)/tod(y); return(z); } inline double divIBD(int x, bool y) { if (tod(x) == tod(y)) z = 1.0d; else z = tod(x)/tod(y); return(z); } inline double diviDBD(double x, bool y) { if (x == tod(y)) z = 1.0d; else z = x/tod(y); return(z); } /* x min y */ inline bool minBBB(bool x, bool y) { return (x&y); } inline int minIII(int x, int y) { if (x <= y) z = x; else z = y; return (z); } inline double minDDD(double x, double y) { if (x <= y) z = x; else z = y; return (z); } inline char minCCC(char x, char y) { if (x <= y) z = x; else z = y; return (z); } inline int minBII(bool x, int y) { if (toi(x) <= y) z = toi(x); else z = y; return (z); } inline int minIBI(int x, bool y) { if (x <= toi(y)) z = x; else z = toi(y); return (z); } inline double minBDD(bool x, double y) { if (tod(x) <= y) z = tod(x); else z = y; return (z); } inline double minDBD(double x, bool y) { if (x <= tod(y)) z = x; else z = tod(y); return (z); } inline double minIDD(int x, double y) { if (tod(x) <= y) z = tod(x); else z = y; return (z); } inline double minDID(double x, int y) { if (x <= tod(y)) z = x; else z = tod(y); return (z); } /* max */ inline bool maxBBB(bool x, bool y) { return (x&y); } inline int maxIII(int x, int y) { if (x <= y) z = y; else z = x; return (z); } inline double maxDDD(double x, double y) { if (x <= y) z = y; else z = x; return (z); } inline char maxCCC(char x, char y) { if (x <= y) z = y; else z = x; return (z); } inline int maxBII(bool x, int y) { if (toi(x) <= y) z = y; else z = toi(x); return (z); } inline int maxIBI(int x, bool y) { if (x <= toi(y)) z = toi(y); else z = x; return (z); } inline double maxBDD(bool x, double y) { if (tod(x) <= y) z = y; else z = tod(x); return (z); } inline double maxDBD(double x, bool y) { if (x <= tod(y)) z = tod(y); else z = x; return (z); } inline double maxIDD(int x, double y) { if (tod(x) <= y) z = y; else z = tod(x); return (z); } inline double maxDID(double x, int y) { if (x <= tod(y)) z = tod(y); else z = x; return (z); } /* x mod y */ inline bool modBBB(bool x, bool y) { return((!x)&y); } inline int modIII(int x, int y) { return(Iresidue(x,y)); } inline double modIII(double x, double y, double QUADct) { return(Dresidue(x,y,QUADct)); } /* x star y */ inline bool starBBB(bool x, bool y) { return(x| !y); } inline int starIBI(int x, bool y) { /* SxS int to Boolean power */ if (y) z = x; else z = 1; return(z); } inline double starDBD(double x, bool y) { /* SxS double to Boolean power */ if (y) z = x; else z = 1.0; return(z); } inline double starIDD(int x, double y) { return(pow(tod(x),y)); } inline double starBID(bool x, int y) { return(pow(tod(x),tod(y))); } inline double starBDD(bool x, double y) { return(pow(tod(x),y)); } inline double starIID(int x, int y) { /* SxS integer power integer */ return(pow(tod(x),tod(y))); } inline double starDID(double x, int y) { return(pow(x,tod(y))); } inline double starDDD(double x, double y) { return(pow(x,y)); } inline double logBDD(bool x, double y) { return(log(tod(x))/log(y)); } inline double logDBD(double x, bool y) { return(log(tod(x))/log(tod(y))); } inline double logBBD(bool x, bool y) { return(log(tod(x))/log(tod(y))); } inline double logIBD(int x, bool y) { return(log(tod(x))/log(tod(y))); } inline double logBID(bool x, int y) { return(log(tod(x))/log(tod(y))); } inline double logIID(int x, int y) { return(log(tod(x))/log(tod(y))); } inline double logDID(int x, double y) { return(log(tod(x))/log(y)); } inline double logIDD(int x, double y) { return(log(tod(x))/log(y)); } inline double logDDD(double x, double y) { return(log(x)/log(y)); } /* NB. APEX Extension of ISO APL to allow comparison of characters */ /* relationals */ inline bool ltBBB(bool x, bool y) { return((!x)&y); } inline bool ltIBB(int x, bool y) { return(x y); } inline bool gtDDB(double x, double y, double QUADct) { return((x > y) && !APEXFUZZEQ(x,y,QUADct)); } inline bool gtCCB(char x, char y) { return(x > y); } inline bool geBBB(bool x, bool y) { return(x | !y); } inline bool geIIB(int x, int y) { return(x >= y); } inline bool geDDB(double x, double y, double QUADct) { return((x >= y) || APEXFUZZEQ(x,y,QUADct)); } inline bool geCCB(char x, char y) { return(x >= y); } /* Boolean functions */ inline bool orBBB(bool x, bool y) { return(x | y); } inline bool andBBB(bool x, bool y) { return(x&y); } /* As of 2004-09-16, we don't support lcm/gcd. Needs work in code generator and dfa to support side effects in scan, etc. rbe */ /* Euclids algorithm for lcm */ #define dandII(XV,YV) (for initial ax := abs(XV); ay := abs(YV); \ u = min(ax,ay); v := max (ax,ay); \ while (v ~= 0) repeat \ v := mod(old u,old v); \ u := old v; \ returns value of (ax*ay)/u \ end for) /* Euclids algorithm for lcm */ #define dandDD(XV,YV) (for initial ax := abs(XV); ay := abs(YV); \ u = min(ax,ay); v := max (ax,ay); \ while (v ~= 0) repeat \ v := mod(old u,old v); \ u := old v; \ returns value of (ax*ay)/u \ end for) /* Euclids algorithm for gcd */ #define dorII(XV,YV) (for initial \ ax := abs(XV); ay := abs(YV); \ u = min(ax,ay); v := max (ax,ay); \ while (v ~= 0) repeat \ v := mod(old u,old v); \ u := old v; \ returns value of u \ end for) /* Euclids algorithm for gcd */ #define dorDDD(XV,YV) (for initial \ ax := abs(XV); ay := abs(YV); \ u = min(ax,ay); v := max (ax,ay); \ while (v ~= 0) repeat \ v := mod(old u,old v); \ u := old v; \ returns value of u \ end for) inline bool nandBBB(bool x, bool y) { return(!(x&y)); } inline bool notBBB(bool x, bool y) { return(!(x|y)); } #define dcircDDD(XV,YV) (if (XV = 1.0d) then sin(YV) \ elseif (XV = 2.0d) then cos(YV) \ elseif (XV = 3.0d) then tan(YV) \ elseif (XV = 4.0d) then pow((1.0d+YV*YV),0.5d) \ else error[double_real] end if) /* domain error check above */ /* 1 circle */ #define dcirc1DDD(XV,YV) (sin(YV)) /* 2 circle */ #define dcirc2DDD(XV,YV) (cos(YV)) /* 3 circle */ #define dcirc3DDD(XV,YV) (tan(YV)) /* 3 circle */ #define dcirc4DDD(XV,YV) (pow((1.0d+YV*YV),0.5d)) inline int[.] rotrXII(int[.] y) { /* Vector reverse */ n = shape(y); cell = 0; z = with( . <= iv <= .) genarray(n, y[(n[0]-1)-iv],cell); return(z); } inline int[+] rotrXII(int[+] y) {/* Last axis reverse on rank>1 */ n = shape(y)[dim(y)-1]; cell = genarray([n], 0); frame = drop([-1],shape(y)); z = with( . <= iv <= .) genarray(frame, rotrXII(y[iv]),cell); return(z); } inline int[*] rhoIII(int[.] x, int[*] y) { z = APEXReshape(toi(x),y); return(z); } inline int[+] rhoIII(int[.] x, int y) { /* Vector reshape scalar to matrix) */ zxrho = prod(toi(x)); /* Result element count */ z = genarray([zxrho], y); /* allocate result */ z = reshape(toi(x),z); return(z); } inline int[.] iotaXII(int shp, int QUADio) { /* Index generator on scalar */ /* HELP! Needs domain check for negative shp */ res = with( . <= [i] <= .) genarray( [toI(shp)], i+QUADio); return( res); } inline int[.] iotaXIINonNegsy(int shp, int QUADio) { /* Index generator on ScalarN when N is non-negative integer */ res = with( . <= [i] <= .) genarray( [toI(shp)], i+QUADio); return( res); } int quadXBB(bool[*] y, int QUADpp, int QUADpw) { /* {quad}{<-} anything */ r=display(y); /* KLUDGE!!! print doesn't support booleans! */ /* MORE kludges - sac tends to evaporate side effects, so we'll help things along abit until that is repaired. */ return(r); } int quadXII(int[*] y, int QUADpp, int QUADpw) { /* {quad}{<-} anything */ r=display(y); /* KLUDGE!!! print doesn't support booleans! */ /* MORE kludges - sac tends to evaporate side effects, so we'll help things along abit until that is repaired. */ return(r); } inline int[2] comaIII(int x, int y) {/* SxS catenate first (or last) axis */ return([toI(x)]++[toI(y)]); } inline int[.] comaIII(int x, int[.] y) {/* SxV catenate first (or last) axis */ return([toI(x)]++toI(y)); } inline bool[2] comaBBB(bool x, bool y) {/* SxS catenate first (or last) axis */ return([toB(x)]++[toB(y)]); } inline bool[.] comaBBB(bool[.] x, bool y) {/* VxS catenate first (or last) axis */ return(toB(x)++[toB(y)]); } inline bool sameIIB(int[+] x, int[+] y) { /* Non-scalar match non-scalar */ if (dim(x) != dim(y)) z = false; else if (! SameShape(x,y)) z = false; else z = with(_mul_SxA_(0,shape(y)) <= iv < shape(y)) fold(andBBB, true, eqIIB(toI(x[iv]), toI(y[iv]))); return(z); } inline bool sameIIB(int[.] x, int[.] y) { /* vector-vector match */ if (dim(x) != dim(y)) z = false; else z = with(_mul_SxA_(0,shape(y)) <= iv < shape(y)) fold(andBBB, true, eqIIB(toI(x[iv]), toI(y[iv]))); return(z); } inline int[*] indrfr(int fr, int i, int[+] X) { /* X[;;;scalari;;;], where i has fr semicolons to its left */ /* Indexing is origin-0. Caller will correct this */ /* This could stand some optimization, perhaps, for boolean i, * unless SAC avoids building an array-valued temp of toI(i). */ frameshape = take([fr],shape(X)); cellshape = drop([fr+1], shape(X)); cell = genarray(cellshape,0); z = with (. <= iv <= .) genarray(frameshape,(X[iv])[i], cell); return(z); } inline int[*] indrfr(int fr, int[+] i, int[+] X) { /* X[;;;i;;;], where i has fr semicolons to its left */ /* Indexing is origin-0. Caller will correct this */ /* This could stand some optimization, perhaps, for boolean i, * unless SAC avoids building an array-valued temp of toI(i). */ frameshape = take([fr],shape(X)); cellshape = shape(i)++drop([fr+1],shape(X)); cell = genarray(cellshape,0); /* not used, but SAC needs help */ z = with (. <= iv <= .) genarray(frameshape,indrfr0(i,X[iv]),cell); return(z); } inline int[*] indrfr0(int i, int[*] X) { /* X[i;;;] i is scalar */ z = X[[i]]; return(z); } inline int[*] indrfr0(int[+] i, int[*] X) { /* X[im;;;] im is array */ cellshape = drop([1],shape(X)); cell = genarray(cellshape, 0); z = with (. <= iv <= .) genarray(shape(i), X[i[iv]],cell); return(z); } inline int[+] indsfr(int fr, int[*] i, int[+] X, int[+] Y) { /* X[;;;i;;;]<- nonscalar Y, where i has fr semicolons to its left */ cellshape = shape(i)++drop([fr],shape(X)); cell = genarray(cellshape,0); /* not used, but SAC needs help */ frameshape = take([fr],shape(X)); z = with (. <= iv <= .) genarray(frameshape,indsfr0(i,X[iv], Y[iv]),cell); zshape = frameshape++cellshape; return(reshape(zshape,z)); } inline int[+] indsfr(int fr, int[+] i, int[+] X, int Y) { /* X[;;;i;;;]<- scalar Y, where i has fr semicolons to its left */ cellshape = drop([fr+1],shape(X)); cell = genarray(cellshape,0); /* not used, but SAC needs help */ frameshape = take([fr],shape(X)); z = with (. <= iv <= .) genarray(frameshape,indsfr0(i,X[iv], Y),cell); return(z); } inline int[+] indsfr(int fr, int i, int[+] X, int Y) { /* X[;;;i;;;]<- scalar Y, where i has fr semicolons to its left */ cellshape = drop([fr+1],shape(X)); cell = genarray(cellshape,0); frameshape = take([fr],shape(X)); z = with (. <= iv <= .) genarray(frameshape,indsfr0(i,X[iv], Y),cell); zshape = frameshape++cellshape; return(z); } inline int[+] indsfr0(int i, int[+] X, int Y) { /* Case 1. X[scalarI;;]<- scalarY NB. Leading axis */ cell = genarray(drop([1],shape(X)),toi(Y)); z = toi(X); z[[i]] = cell; return(z); } inline int[+] indsfr0(int[.] iv, int[+] X, int Y) { /* 2. X[non-scalarIV;;]<- scalarY NB. Leading axis */ /* This would almost work under a with-loop, but the potential * for duplicates in iv scuppers that. Ergo, FOR loop. */ z = toi(X); cellshape = drop([1],shape(X)); cell = genarray (cellshape, toi(Y)); raveli = APEXRavel(iv); for(i=0; i