From b69bef8f327a3aa71cb85d342857dcdea694d1a7 Mon Sep 17 00:00:00 2001 From: AJ Rossini Date: Thu, 4 Oct 2012 15:45:53 +0200 Subject: [PATCH] removing liblispstat. We don't need this anymore, code is worn out --- lib/Makefile | 39 -- lib/Makefile.macosx | 41 -- lib/betabase.c | 282 ------------- lib/bivnor.c | 131 ------ lib/cbayes.c | 241 ----------- lib/cdists.c | 859 --------------------------------------- lib/cffi-glue.c | 477 ---------------------- lib/cfft.c | 827 -------------------------------------- lib/cholesky.c | 79 ---- lib/clib.make | 134 ------- lib/clinalg.c | 218 ---------- lib/cmpinclude.h | 695 -------------------------------- lib/complex.c | 131 ------ lib/complex.h | 18 - lib/derivatives.c | 71 ---- lib/eigen.c | 625 ----------------------------- lib/functions.c | 1115 --------------------------------------------------- lib/gamln.c | 35 -- lib/gammabase.c | 289 ------------- lib/kernel.c | 91 ----- lib/linalg.h | 27 -- lib/linalgdata.c | 92 ----- lib/lowess.c | 170 -------- lib/ludecomp.c | 176 -------- lib/makerotation.c | 54 --- lib/minimize.c | 868 --------------------------------------- lib/nor.c | 100 ----- lib/num_sfun.c | 655 ------------------------------ lib/ppnd.c | 60 --- lib/qrdecomp.c | 186 --------- lib/rcondest.c | 102 ----- lib/splines.c | 112 ------ lib/studentbase.c | 125 ------ lib/svdecomp.c | 263 ------------ lib/xlisp.h | 9 - lib/xmath.h | 5 - 36 files changed, 9402 deletions(-) delete mode 100644 lib/Makefile delete mode 100644 lib/Makefile.macosx delete mode 100644 lib/betabase.c delete mode 100644 lib/bivnor.c delete mode 100644 lib/cbayes.c delete mode 100644 lib/cdists.c delete mode 100644 lib/cffi-glue.c delete mode 100644 lib/cfft.c delete mode 100644 lib/cholesky.c delete mode 100644 lib/clib.make delete mode 100644 lib/clinalg.c delete mode 100644 lib/cmpinclude.h delete mode 100644 lib/complex.c delete mode 100644 lib/complex.h delete mode 100644 lib/derivatives.c delete mode 100644 lib/eigen.c delete mode 100644 lib/functions.c delete mode 100644 lib/gamln.c delete mode 100644 lib/gammabase.c delete mode 100644 lib/kernel.c delete mode 100644 lib/linalg.h delete mode 100644 lib/linalgdata.c delete mode 100644 lib/lowess.c delete mode 100644 lib/ludecomp.c delete mode 100644 lib/makerotation.c delete mode 100644 lib/minimize.c delete mode 100644 lib/nor.c delete mode 100644 lib/num_sfun.c delete mode 100644 lib/ppnd.c delete mode 100644 lib/qrdecomp.c delete mode 100644 lib/rcondest.c delete mode 100644 lib/splines.c delete mode 100644 lib/studentbase.c delete mode 100644 lib/svdecomp.c delete mode 100644 lib/xlisp.h delete mode 100644 lib/xmath.h diff --git a/lib/Makefile b/lib/Makefile deleted file mode 100644 index fcc0ee6..0000000 --- a/lib/Makefile +++ /dev/null @@ -1,39 +0,0 @@ -## -*- mode: Makefile -*- - -LISPSTAT-SO = liblispstat.so - -CDISTOBJS = betabase.o bivnor.o cdists.o gamln.o gammabase.o nor.o ppnd.o \ - studentbase.o - -CLINALGOBJS = cholesky.o linalgdata.o ludecomp.o complex.o svdecomp.o \ - qrdecomp.o cfft.o clinalg.o rcondest.o kernel.o lowess.o splines.o \ - makerotation.o eigen.o - -CBAYESOBJS = cbayes.o functions.o derivatives.o minimize.o - -CGLUE = cffi-glue.o - -## CFLAGS = -g -DINTPTR -fPIC -Wall -CFLAGS = -g -Wall - -cffi : $(LISPSTAT-SO) - -sbcl : clib.a - -kcl: clib.a - -excl: clib.a exclglue.o - -clib.a: ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - ar rcv clib.a ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - ranlib clib.a - -$(LISPSTAT-SO) :${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - gcc -shared -o $(LISPSTAT-SO) ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - -clean: - rm -f ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) *~ - -cleanall: - make clean - rm -f clib.a exclglue.o diff --git a/lib/Makefile.macosx b/lib/Makefile.macosx deleted file mode 100644 index 771a67b..0000000 --- a/lib/Makefile.macosx +++ /dev/null @@ -1,41 +0,0 @@ -## -*- mode: Makefile -*- -ARCHS = -arch i386 -arch x86_64 -SPEC_WARNINGS = -Wconversion -Wformat -Wshorten-64-to-32 - -LISPSTAT-SO = liblispstat.dylib - -CDISTOBJS = betabase.o bivnor.o cdists.o gamln.o gammabase.o nor.o ppnd.o \ - studentbase.o - -CLINALGOBJS = cholesky.o linalgdata.o ludecomp.o complex.o svdecomp.o \ - qrdecomp.o cfft.o clinalg.o rcondest.o kernel.o lowess.o splines.o \ - makerotation.o eigen.o - -CBAYESOBJS = cbayes.o functions.o derivatives.o minimize.o - -CGLUE = cffi-glue.o - -## CFLAGS = -g -DINTPTR -fPIC -Wall -CFLAGS = -g -Wall $(ARCHS) $(SPEC_WARNINGS) - -cffi : $(LISPSTAT-SO) - -sbcl : clib.a - -kcl: clib.a - -excl: clib.a exclglue.o - -clib.a: ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - ar rcv clib.a ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - ranlib clib.a - -$(LISPSTAT-SO) :${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - gcc $(ARCHS) -dynamiclib -single_module -undefined dynamic_lookup -o $(LISPSTAT-SO) ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - -clean: - rm -f ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) *~ - -cleanall: - make clean - rm -f clib.a exclglue.o diff --git a/lib/betabase.c b/lib/betabase.c deleted file mode 100644 index 7e5e6d4..0000000 --- a/lib/betabase.c +++ /dev/null @@ -1,282 +0,0 @@ -#include "xmath.h" - -#define TRUE 1 -#define FALSE 0 - -/* external function */ -extern double macheps(), gamma(); - -/* forward declarations */ -static double logbeta(), betai(), xinbta(); - -void -betabase(double *x, double *a, double *b, long *gia, long *gib, double *cdf) -{ - *cdf = betai(*x, *a, *b); -} - -double ppbeta(double p, double a, double b, int *ifault) -{ - double lbeta; - - lbeta = gamma(a) + gamma(b) - gamma(a + b); - return(xinbta(&a, &b, &lbeta, &p, ifault)); -} - -/* - Static routines -*/ - -static double logbeta(double p, double q) -{ - return(gamma(p) + gamma(q) - gamma(p + q)); -} - -#define Min(x,y) (((x) < (y)) ? (x) : (y)) -#define Max(x,y) (((x) > (y)) ? (x) : (y)) - -static double betai(double x, double pin, double qin) -{ - /* Translated from FORTRAN - july 1977 edition. w. fullerton, c3, los alamos scientific lab. - based on bosten and battiste, remark on algorithm 179, comm. acm, - v 17, p 153, (1974). - - input arguments -- - x upper limit of integration. x must be in (0,1) inclusive. - p first beta distribution parameter. p must be gt 0.0. - q second beta distribution parameter. q must be gt 0.0. - betai the incomplete beta function ratio is the probability that a - random variable from a beta distribution having parameters - p and q will be less than or equal to x. - */ - double c, finsum, p, ps, q, term, xb, xi, y, dbetai, p1; - static double eps = 0.0, alneps = 0.0, sml = 0.0, alnsml = 0.0; - long i, n, ib; - - /* I'm not sure these tolerances are optimal */ - if (eps == 0.0) { - eps = macheps(); - alneps = log(eps); - sml = macheps(); - alnsml = log(sml); - } - - y = x; - p = pin; - q = qin; - if (q > p || x >= 0.8) - if (x >= 0.2) { - y = 1.0 - y; - p = qin; - q = pin; - } - - if ((p + q) * y / (p + 1.0) < eps) { - dbetai = 0.0; - xb = p * log(Max(y, sml)) - log(p) - logbeta(p, q); - if (xb > alnsml && y != 0.0) dbetai = exp(xb); - if (y != x || p != pin) dbetai = 1.0 - dbetai; - } - else { - /* - * evaluate the infinite sum first. term will equal - * y**p/beta(ps,p) * (1.-ps)-sub-i * y**i / fac(i) . - */ - ps = q - floor(q); - if (ps == 0.0) ps = 1.0; - xb = p * log(y) - logbeta(ps, p) - log(p); - dbetai = 0.0; - if (xb >= alnsml) { - - dbetai = exp(xb); - term = dbetai * p; - if (ps != 1.0) { - n = Max(alneps / log(y), 4.0); - for (i = 1; i <= n; i++) { - xi = i; - term = term * (xi - ps) * y / xi; - dbetai = dbetai + term / (p + xi); - } - } - } - /* - * now evaluate the finite sum, maybe. - */ - if (q > 1.0) { - - xb = p * log(y) + q * log(1.0 - y) - logbeta(p,q) - log(q); - ib = Max(xb / alnsml, 0.0); - term = exp(xb - ((float) ib) * alnsml); - c = 1.0 / (1.0 - y); - p1 = q * c / (p + q - 1.0); - - finsum = 0.0; - n = q; - if (q == (float) n) n = n - 1; - for (i = 1; i <= n; i++) { - if (p1 <= 1.0 && term / eps <= finsum) break; - xi = i; - term = (q - xi + 1.0) * c * term / (p + q - xi); - - if (term > 1.0) ib = ib - 1; - if (term > 1.0) term = term * sml; - - if (ib==0) finsum = finsum + term; - } - - dbetai = dbetai + finsum; - } - if (y != x || p != pin) dbetai = 1.0 - dbetai; - dbetai = Max(Min(dbetai, 1.0), 0.0); - } - return(dbetai); -} - -/* - xinbta.f -- translated by f2c and modified - - algorithm as 109 appl. statist. (1977), vol.26, no.1 - (replacing algorithm as 64 appl. statist. (1973), vol.22, no.3) - - Remark AS R83 has been incorporated in this version. - - Computes inverse of the incomplete beta function - ratio for given positive values of the arguments - p and q, alpha between zero and one. - log of complete beta function, beta, is assumed to be known. - - Auxiliary function required: betai - - SAE below is the most negative decimal exponent which does not - cause an underflow; a value of -308 or thereabouts will often be -*/ - -static double xinbta(double *p, double *q, double *beta, double *alpha, int *ifault) -{ - /* Initialized data */ - static double sae = -30.0; /* this should be sufficient */ - static double zero = 0.0; - static double one = 1.0; - static double two = 2.0; - static double three = 3.0; - static double four = 4.0; - static double five = 5.0; - static double six = 6.0; - - /* System generated locals */ - double ret_val, d_1, d_2; - - /* Local variables */ - static int indx; - static double prev, a, g, h, r, s, t, w, y, yprev, pp, qq; - static double sq, tx, adj, acu; - static long iex; - static double fpu, xin; - - /* Define accuracy and initialise. */ - fpu = sae * 10.; - ret_val = *alpha; - - /* test for admissibility of parameters */ - *ifault = 1; - if (*p <= zero || *q <= zero) return ret_val; - *ifault = 2; - if (*alpha < zero || *alpha > one) return ret_val; - *ifault = 0; - if (*alpha == zero || *alpha == one) return ret_val; - - /* change tail if necessary */ - if (*alpha <= .5) { - a = *alpha; - pp = *p; - qq = *q; - indx = FALSE; - } - else { - a = one - *alpha; - pp = *q; - qq = *p; - indx = TRUE; - } - - /* calculate the initial approximation */ - r = sqrt(-log(a * a)); - y = r - (r * .27061 + 2.30753) / (one + (r * .04481 + .99229) * r); - if (pp > one && qq > one) { - r = (y * y - three) / six; - s = one / (pp + pp - one); - t = one / (qq + qq - one); - h = two / (s + t); - d_1 = y * sqrt(h + r) / h; - d_2 = (t - s) * (r + five / six - two / (three * h)); - w = d_1 - d_2; - ret_val = pp / (pp + qq * exp(w + w)); - } - else { - r = qq + qq; - t = one / (qq * 9.); - /* Computing 3rd power */ - d_1 = one - t + y * sqrt(t); - t = r * (d_1 * d_1 * d_1); - if (t <= zero) { - ret_val = one - exp((log((one - a) * qq) + *beta) / qq); - } - else { - t = (four * pp + r - two) / t; - if (t <= one) ret_val = exp((log(a * pp) + *beta) / pp); - else ret_val = one - two / (t + one); - } - } - - - /* - solve for x by a modified newton-raphson method, using the function betai - */ - r = one - pp; - t = one - qq; - yprev = zero; - sq = one; - prev = one; - if (ret_val < 1e-4) ret_val = 1e-4; - if (ret_val > .9999) ret_val = .9999; - /* Computing MAX, two 2nd powers */ - d_1 = -5.0 / (pp * pp) - 1.0 / (a * a) - 13.0; - iex = (sae > d_1) ? sae : d_1; - acu = pow(10.0, (double) iex); - do { - y = betai(ret_val, pp, qq); - if (*ifault != 0) { - *ifault = 3; - return ret_val; - } - xin = ret_val; - y = (y - a) * exp(*beta + r * log(xin) + t * log(one - xin)); - if (y * yprev <= zero) { - prev = (sq > fpu) ? sq : fpu; - } - g = one; - do { - adj = g * y; - sq = adj * adj; - if (sq < prev) { - tx = ret_val - adj; - if (tx >= zero && tx <= one) { - if (prev <= acu || y * y <= acu) { - if (indx) ret_val = one - ret_val; - return ret_val; - } - if (tx != zero && tx != one) break; - } - } - g /= three; - } while (TRUE); - if (tx == ret_val) { - if (indx) ret_val = one - ret_val; - return ret_val; - } - ret_val = tx; - yprev = y; - } while (TRUE); - return ret_val; -} diff --git a/lib/bivnor.c b/lib/bivnor.c deleted file mode 100644 index d10f573..0000000 --- a/lib/bivnor.c +++ /dev/null @@ -1,131 +0,0 @@ -#include "xmath.h" - -#define twopi 6.283195307179587 -#define con (twopi / 2.0) * 10.0e-10 - -extern void normbase(double *, double *); - -double bivnor(double ah, double ak, double r) -{ - /* - based on alg 4628 comm. acm oct 73 - gives the probability that a bivariate normal exceeds (ah,ak). - gh and gk are .5 times the right tail areas of ah, ak under a n(0,1) - - Tranlated from FORTRAN to ratfor using struct; from ratfor to C by hand. - */ - double a2, ap, b, cn, conex, ex, g2, gh, gk, gw, h2, h4, rr, s1, s2, - sgn, sn, sp, sqr, t, temp, w2, wh, wk; - int is; - - temp = -ah; - normbase(&temp, &gh); - gh = gh / 2.0; - temp = -ak; - normbase(&temp, &gk); - gk = gk / 2.0; - - b = 0; - - if (r==0) - b = 4*gh*gk; - else { - rr = 1-r*r; - if (rr<0) - return(-1.0); - if (rr!=0) { - sqr = sqrt(rr); - if (ah!=0) { - b = gh; - if (ah*ak<0) - b = b-.5; - else if (ah*ak==0) - goto label10; - } - else if (ak==0) { - b = atan(r/sqr)/twopi+.25; - goto label50; - } - b = b+gk; - if (ah==0) - goto label20; - label10: - wh = -ah; - wk = (ak/ah-r)/sqr; - gw = 2*gh; - is = -1; - goto label30; - label20: - do { - wh = -ak; - wk = (ah/ak-r)/sqr; - gw = 2*gk; - is = 1; - label30: - sgn = -1; - t = 0; - if (wk!=0) { - if (fabs(wk)>=1) { - if (fabs(wk)==1) { - t = wk*gw*(1-gw)/2; - goto label40; - } else { - sgn = -sgn; - wh = wh*wk; - normbase(&wh, &g2); - wk = 1/wk; - if (wk<0) { - b = b+.5; - } - b = b-(gw+g2)/2+gw*g2; - } - } - h2 = wh*wh; - a2 = wk*wk; - h4 = h2*.5; - ex = 0; - if (h4<150.0) { - ex = exp(-h4); - } - w2 = h4*ex; - ap = 1; - s2 = ap-ex; - sp = ap; - s1 = 0; - sn = s1; - conex = fabs(con/wk); - do { - cn = ap*s2/(sn+sp); - s1 = s1+cn; - if (fabs(cn)<=conex) - break; - sn = sp; - sp = sp+1; - s2 = s2-w2; - w2 = w2*h4/sp; - ap = -ap*a2; - } while (1); - t = (atan(wk)-wk*s1)/twopi; - label40: - b = b+sgn*t; - } - if (is>=0) - break; - } while(ak!=0); - } - else if (r>=0) - if (ah>=ak) - b = 2*gh; - else - b = 2*gk; - else if (ah+ak<0) - b = 2*(gh+gk)-1; - } - label50: - if (b<0) - b = 0; - if (b>1) - b = 1; - - return(b); -} diff --git a/lib/cbayes.c b/lib/cbayes.c deleted file mode 100644 index 8a842a4..0000000 --- a/lib/cbayes.c +++ /dev/null @@ -1,241 +0,0 @@ -/* cbayes - Lisp interface to laplace approximation stuff */ -/* Copyright (c) 1990, by Luke Tierney */ - -#include /* for calloc/realloc */ -#include "linalg.h" - -#ifdef INTPTR -typedef int PTR; -#else -typedef char *PTR; -#endif - -extern void maximize_callback(size_t, PTR, PTR, PTR, PTR, PTR); -extern void evalfront(char **, size_t *, double *, double *, double *, - double *, double *, double *); - -extern void maxfront(); -extern void bufputstr(char *); - -/************************************************************************/ -/** **/ -/** Definitions and Globals **/ -/** **/ -/************************************************************************/ - -#define MAXALLOC 20 - -static char *mem[MAXALLOC], memcount; - -typedef struct { - int n, m, k, itnlimit, backtrack, verbose, vals_suppl, exptilt; - int count, termcode; -} MaxIPars; - -typedef struct { - double typf, h, gradtol, steptol, maxstep, dflt, tilt, newtilt, hessadd; -} MaxDPars; - -typedef struct { - MaxIPars max; - int full, covar; -} MomIPars; - -typedef struct { - MaxDPars max; - double mgfdel; -} MomDPars; - -/************************************************************************/ -/** **/ -/** Fake Replacements for S Interface **/ -/** **/ -/************************************************************************/ - -static void -meminit() -{ - static int inited = FALSE; - int i; - - if (! inited) { - for (i = 0; i < MAXALLOC; i++) mem[i] = nil; - inited = TRUE; - } - - memcount = 0; -} - -static int -makespace(void **pptr, size_t size) /* why are we using **char? */ -{ - if (size <= 0) { - return(1); /* we've done, by default, what we asked for. */ - } - if (*pptr == nil) { - *pptr = calloc(size, 1); - } else { - *pptr = realloc(*pptr, size); - } - if (size > 0 && *pptr == nil) { - return(0); /* xlfail("memory allocation failed"); FIXME:AJR xlfail redef. */ - } - return(1); -} - - -/************************************************************************/ -/** **/ -/** Callback Function **/ -/** **/ -/************************************************************************/ - -static void -callLminfun(size_t n, - double *x, double *fval, double *grad, double *hess, - long *derivs) -{ - maximize_callback(n, (PTR) x, - (PTR) fval, (PTR) grad, (PTR) hess, (PTR) derivs); -} - -void -call_S(char *fun, long narg, char **args, char **mode, long *length,char **names, - long nvals, char **values) -{ - long n = length[0], derivs; - static double *fval = nil, *grad = nil, *hess = nil; - - makespace((void **)&fval, 1 * sizeof(double)); /* probably should test the - result of this and the next - 2 to make sure that they are - appropriate */ - makespace((void **)&grad, n * sizeof(double)); - makespace((void **)&hess, n * n * sizeof(double)); - - callLminfun(n,(double *)args[0], fval, grad, hess, &derivs); - - values[0] = (char *) fval; - values[1] = (derivs > 0) ? (char *) grad : nil; - values[2] = (derivs > 1) ? (char *) hess : nil; -} - -void -Recover(s, w) - char *s, *w; -{ - /* FIXME:AJR: xlfail(s); */ - return; -} - -/************************************************************************/ -/** **/ -/** Numerical Derivatives **/ -/** **/ -/************************************************************************/ - -void -numgrad_front(size_t n, PTR px, PTR pgrad, double h, PTR pscale) -{ - LVAL f = nil; - double fval; - - evalfront((char **)&f, &n, (double *) px, - &fval, (double *) pgrad, nil, &h, (double *) pscale); -} - -void -numhess_front(size_t n, PTR px, PTR pf, PTR pgrad, PTR phess, double h, PTR pscale) -{ - LVAL f = nil; - - evalfront((char **)&f, &n, (double *) px, - (double *) pf, (double *) pgrad, (double *) phess, - &h, (double *) pscale); -} - -/************************************************************************/ -/** **/ -/** Maximization Interface **/ -/** **/ -/************************************************************************/ - -/* internals array information */ -#define INTLEN 12 -#define F_POS 0 -#define G_POS 1 -#define C_POS 2 -#define X_POS 3 -#define SCALE_POS 4 -#define FVALS_POS 5 -#define CVALS_POS 6 -#define CTARG_POS 7 -#define IPARS_POS 8 -#define DPARS_POS 9 -#define TSCAL_POS 10 -#define MULT_POS 11 - -static MaxIPars getMaxIPars(int *ipars) -{ - MaxIPars ip; - - ip.n = ipars[0]; - ip.m = ipars[1]; - ip.k = ipars[2]; - ip.itnlimit = ipars[3]; - ip.backtrack = ipars[4]; - ip.verbose = ipars[5]; - ip.vals_suppl = ipars[6]; - ip.exptilt = ipars[7]; - ip.count = ipars[8]; - ip.termcode = ipars[9]; - - return(ip); -} - -static MaxDPars getMaxDPars(dpars) - double *dpars; -{ - MaxDPars dp; - - dp.typf = dpars[0]; - dp.h = dpars[1]; - dp.gradtol = dpars[2]; - dp.steptol = dpars[3]; - dp.maxstep = dpars[4]; - dp.dflt = dpars[5]; - dp.tilt = dpars[6]; - dp.newtilt = dpars[7]; - dp.hessadd = dpars[8]; - - return(dp); -} - -void -minfo_maximize(PTR px, PTR pfvals, PTR pscale, PTR pip, PTR pdp, int verbose) -{ - LVAL f = nil; - MaxIPars ip; - MaxDPars dp; - int n, m, k; - static double *dx, *typx, *fvals; - char *msg; - - dx = (double *) px; - typx = (double *) pscale; - fvals = (double *) pfvals; - - ip = getMaxIPars((int *) pip); - dp = getMaxDPars((double *) pdp); - - m = 0; - k = 0; - n = ip.n; - if (verbose >= 0) ip.verbose = verbose; - - meminit(); - maxfront(&f, nil, nil, dx, typx, fvals, nil, nil, nil, - &ip, &dp, nil, &msg); - - bufputstr(msg); -} diff --git a/lib/cdists.c b/lib/cdists.c deleted file mode 100644 index e8cad37..0000000 --- a/lib/cdists.c +++ /dev/null @@ -1,859 +0,0 @@ -#include "xmath.h" - -extern void normbase(double *,double *); -extern void gammabase(double *,double *, double *); -extern void studentbase(double *,double *, double *); -extern void betabase(double *, double *,double *,long *, long *,double *); -extern long max(long, long); -extern long min(long, long); - -extern void xlfail(char *); - -extern double ppnd(), gamma(), bivnor(), uni(), ppgamma(), ppbeta(), - ppstudent(); - -/* forward declaration */ -extern double tdens(); - -#ifndef PI -#define PI 3.14159265358979323846 -#endif /* PI*/ -#ifdef __LP64__ -#define TRUE 1L -#define FALSE 0L -#else -#define TRUE 1 -#define FALSE 0 -#endif - -/* - * Under ULTRIX 3.1 (the cc1.31 compilers in particular) the _G0 math - * library does not really exist! You either need to figure out how - * to get tan() and floor() in at load time for kcl, or use the ones - * here. - */ - -#ifdef mips -#ifdef UX31 -double tan(double x) -{ - return(sin(x) / cos(x)); -} - -double floor(double x) -{ - long ix = x; - double dx = ix; - return((dx <= x) ? dx : dx - 1.0); -} -#endif -#endif - -static void -checkflag(int flag) -{ - /* do nothing for now */ -} - -static void -checkexp(double a) -{ - if (a <= 0.0) { - xlfail("non-positive gamma or beta exponent"); - } -} - -static void -checkdf(double df) -{ - if (df <= 0.0) { - xlfail("non-positive degrees of freedom"); - } -} - -static void -checkprob(double p, long zerostrict, long onestrict) -{ - if (zerostrict) { - if (p <= 0.0) xlfail("non-positive probability argument"); - } else { - if (p < 0.0) xlfail("negative probability argument"); - } - if (onestrict) { - if (p >= 1.0) xlfail("probability argument not less than one"); - } else { - if (p > 1.0) xlfail("probability argument greater than one"); - } -} - -static void -checkrho(double r) -{ - if (r < -1 || r > 1) { - xlfail("correlation out of range"); - } -} - -static void -checkpoisson(double L) -{ - if (L < 0.0) { - xlfail("negative Poisson mean"); - } -} - -static void -checkbinomial(long n, double p) -{ - if (p < 0.0 || p > 1.0) xlfail("binomial p out of range"); - if (n < 1) xlfail("non-positive binomial n"); -} - -/*****************************************************************************/ -/*****************************************************************************/ -/** **/ -/** Uniform Distribution **/ -/** **/ -/*****************************************************************************/ -/*****************************************************************************/ - -/* uniform generator - avoids zero and one */ -double unirand() -{ - double u; - do { - u = uni(); - } while ((u <= 0.0) || (u >= 1.0)); - return(u); -} - -/*****************************************************************************/ -/*****************************************************************************/ -/** **/ -/** Normal Distribution **/ -/** **/ -/*****************************************************************************/ -/*****************************************************************************/ - -/* standard normal cdf */ -double normalcdf(double x) -{ - double p; - normbase(&x, &p); - return(p); -} - -/* standard normal quantile function */ -double normalquant(double p) -{ - int flag; - double x; - - checkprob(p, TRUE, TRUE); - x = ppnd(p, &flag); - checkflag(flag); - return(x); -} - -/* standard normal density */ -double normaldens(double x) -{ - return(exp(- 0.5 * x * x) / sqrt(2.0 * PI)); -} - -/* standard normal generator */ -double normalrand() -{ - double x, y, u, u1, v; - static double c = -1.0; - - if (c < 0.0) c = sqrt(2.0 / exp(1.0)); - - /* ratio of uniforms with linear pretest */ - do { - u = unirand(); - u1 = unirand(); - v = c * (2 * u1 - 1); - x = v / u; - y = x * x / 4.0; - } while(y > (1 - u) && y > - log(u)); - return(x); -} - -/* bivariate normal cdf */ -double bnormcdf(double x, double y, double r) -{ - checkrho(r); - return(bivnor(-x, -y, r)); -} - -/*****************************************************************************/ -/*****************************************************************************/ -/** **/ -/** Cauchy Distribution **/ -/** **/ -/*****************************************************************************/ -/*****************************************************************************/ - -double -cauchycdf(double dx) -{ - return((atan(dx) + PI / 2) / PI); -} - -double -cauchyquant(double p) -{ - checkprob(p, TRUE, TRUE); - return(tan(PI * (p - 0.5))); -} - -double -cauchydens(double dx) -{ - return(tdens(dx, 1.0)); -} - -/* cauchy generator */ -double -cauchyrand() -{ - double u1, u2, v1, v2; - - /* ratio of uniforms on half disk */ - do { - u1 = unirand(); - u2 = unirand(); - v1 = 2.0 * u1 - 1.0; - v2 = u2; - } while(v1 * v1 + v2 * v2 > 1.0); - return(v1 / v2); -} - -/*****************************************************************************/ -/** **/ -/** Gamma Distribution **/ -/** **/ -/*****************************************************************************/ - -double -gammacdf(double x, double a) -{ - double p; - - checkexp(a); - if (x <= 0.0) p = 0.0; - else gammabase(&x, &a, &p); - return(p); -} - -double -gammaquant(double p, double a) -{ - int flag; - double x; - - checkexp(a); - checkprob(p, FALSE, TRUE); - x = ppgamma(p, a, &flag); - checkflag(flag); - return(x); -} - -/* gamma density */ -double gammadens(double x, double a) -{ - double dens; - - checkexp(a); - if (x <= 0.0) { - dens = 0.0; - } else { - dens = exp(log(x) * (a - 1) - x - gamma(a)); - } - return(dens); -} - -/* gamma generator */ -double -gammarand(double a) -{ - double x, u0, u1, u2, v, w, c, c1, c2, c3, c4, c5; - static double e = -1.0; - long done; - - checkexp(a); - if (e < 0.0) e = exp(1.0); - - if (a < 1.0) { - /* Ahrens and Dieter algorithm */ - done = FALSE; - c = (a + e) / e; - do { - u0 = unirand(); - u1 = unirand(); - v = c * u0; - if (v <= 1.0) { - x = exp(log(v) / a); - if (u1 <= exp(-x)) done = TRUE; - } - else { - x = -log((c - v) / a); - if (x > 0.0 && u1 < exp((a - 1.0) * log(x))) done = TRUE; - } - } while(! done); - } - else if (a == 1.0) x = -log(unirand()); - else { - /* Cheng and Feast algorithm */ - c1 = a - 1.0; - c2 = (a - 1.0 / (6.0 * a)) / c1; - c3 = 2.0 / c1; - c4 = 2.0 / (a - 1.0) + 2.0; - c5 = 1.0 / sqrt(a); - do { - do { - u1 = unirand(); - u2 = unirand(); - if (a > 2.5) u1 = u2 + c5 * (1.0 - 1.86 * u1); - } while (u1 <= 0.0 || u1 >= 1.0); - w = c2 * u2 / u1; - } while ((c3 * u1 + w + 1.0/w) > c4 && (c3 * log(u1) - log(w) + w) > 1.0); - x = c1 * w; - } - return(x); -} - -/*****************************************************************************/ -/** **/ -/** Chi-Square Distribution **/ -/** **/ -/*****************************************************************************/ - -double chisqcdf(double x, double df) -{ - double p, a; - - checkdf(df); - a = 0.5 * df; - x = 0.5 * x; - if (x <= 0.0) { - p = 0.0; - } else { - gammabase(&x, &a, &p); - } - return(p); -} - -double -chisqquant(double p, double df) -{ - double x, a; - int flag; - - checkdf(df); - checkprob(p, FALSE, TRUE); - a = 0.5 * df; - x = 2.0 * ppgamma(p, a, &flag); - checkflag(flag); - return(x); -} - -double -chisqdens(double dx, double da) -{ - checkdf(da); - da = 0.5 * da; - dx = 0.5 * dx; - return(0.5 * gammadens(dx, da)); -} - -double -chisqrand(double df) -{ - checkdf(df); - return(2.0 * gammarand(df / 2.0)); -} - -/*****************************************************************************/ -/** **/ -/** Beta Distribution **/ -/** **/ -/*****************************************************************************/ - -double -betacdf(double x, double a, double b) -{ - double p; - long ia, ib; - - checkexp(a); checkexp(b); - ia = a; ib = b; - if (x <= 0.0) { - p = 0.0; - } else { - if (x >= 1.0) { - p = 1.0; - } else { - betabase(&x, &a, &b, &ia, &ib, &p); - } - } - return(p); -} - -double -betaquant(double p, double a, double b) -{ - double x; - int flag; - - checkexp(a); - checkexp(b); - checkprob(p, FALSE, FALSE); - x = ppbeta(p, a, b, &flag); - checkflag(flag); - return(x); -} - -static double -logbeta(double a, double b) -{ - static double da = 0.0, db = 0.0, lbeta = 0.0; - - if (a != da || b != db) { /* cache most recent call */ - da = a; db = b; - lbeta = gamma(da) + gamma(db) - gamma(da + db); - } - return(lbeta); -} - -/* beta density */ -double -betadens(double x, double a, double b) -{ - double dens; - - checkexp(a); - checkexp(b); - if (x <= 0.0 || x >= 1.0) { - dens = 0.0; - } else { - dens = exp(log(x) * (a - 1) + log(1 - x) * (b - 1) - logbeta(a, b)); - } - return(dens); -} - -/* beta generator */ -double -betarand(double a, double b) -{ - double x, y; - - checkexp(a); - checkexp(b); - x = gammarand(a); - y = gammarand(b); - return(x / (x + y)); -} - -/*****************************************************************************/ -/** **/ -/** t Distribution **/ -/** **/ -/*****************************************************************************/ - -/* t cdf */ -double -tcdf(double x, double df) -{ - double p; - - checkdf(df); - studentbase(&x, &df, &p); - return(p); -} - -double -tquant(double p, double df) -{ - double x; - int flag; - - checkdf(df); - checkprob(p, TRUE, TRUE); - x = ppstudent(p, df, &flag); - checkflag(flag); - return(x); -} - -double -tdens(double x, double a) -{ - double dens; - - checkdf(a); - dens = (1.0 / sqrt(a * PI)) - * exp(gamma(0.5 * (a + 1)) - gamma(0.5 * a) - - 0.5 * (a + 1) * log(1.0 + x * x / a)); - return(dens); -} - -double trand(double df) -{ - checkdf(df); - return(normalrand() / sqrt(chisqrand(df) / df)); -} - -/*****************************************************************************/ -/** **/ -/** F Distribution **/ -/** **/ -/*****************************************************************************/ - -double fcdf(double x, double ndf, double ddf) -{ - double p, a, b; - - checkdf(ndf); checkdf(ddf); - a = 0.5 * ddf; - b = 0.5 * ndf; - if (x <= 0.0) { - p = 0.0; - } else { - x = a / (a + b * x); - p = 1.0 - betacdf(x, a, b); - } - return(p); -} - -double fquant(double p, double ndf, double ddf) -{ - double x, a, b; - int flag; - - checkdf(ndf); checkdf(ddf); - checkprob(p, FALSE, TRUE); - a = 0.5 * ddf; - b = 0.5 * ndf; - if (p == 0.0) { - x = 0.0; - } else { - p = 1.0 - p; - x = ppbeta(p, a, b, &flag); - checkflag(flag); - x = a * (1.0 / x - 1.0) / b; - } - return(x); -} - -double -fdens(double dx, double da, double db) -{ - double dens; - - checkdf(da); - checkdf(db); - if (dx <= 0.0) { - dens = 0.0; - } else { - dens = exp(0.5 * da * log(da) + 0.5 * db *log(db) - + (0.5 * da - 1.0) * log(dx) - - logbeta(0.5 * da, 0.5 * db) - - 0.5 * (da + db) * log(db + da * dx)); - } - return(dens); -} - -/* f generator */ -double frand(double ndf, double ddf) -{ - checkdf(ndf); - checkdf(ddf); - return((ddf * chisqrand(ndf)) / (ndf * chisqrand(ddf))); -} - -/*****************************************************************************/ -/** **/ -/** Poisson Distribution **/ -/** **/ -/*****************************************************************************/ - -static double -poisson_cdf(long k, double L) -{ - double dp, dx; - - if (k < 0) { - dp = 0.0; - } else { - if (L == 0.0) { - dp = (k < 0) ? 0.0 : 1.0; - } else { - dx = k + 1.0; - gammabase(&L, &dx, &dp); - dp = 1.0 - dp; - } - } - return(dp); -} - -double -poissoncdf(double k, double L) -{ - checkpoisson(L); - return(poisson_cdf((long) floor(k), L)); -} - -long -poissonquant(double x, double L) -{ - long k, k1, k2, del, ia; - double m, s, p1, p2, pk; - - checkpoisson(L); - checkprob(x, FALSE, TRUE); - m = L; - s = sqrt(L); - del = max(1, (long) (0.2 * s)); - - if (x == 0.0) { - k = 0.0; - } else { - k = m + s * ppnd(x, &ia); - } - k1 = k; - k2 = k; - - do { - k1 = k1 - del; k1 = max(0, k1); - p1 = poisson_cdf(k1, L); - } while (k1 > 0 && p1 > x); - - if (k1 == 0 && p1 >= x) { - return(k1); - } - - do { - k2 = k2 + del; - p2 = poisson_cdf(k2, L); - } while (p2 < x); - - while (k2 - k1 > 1) { - k = (k1 + k2) / 2; - pk = poisson_cdf(k, L); - if (pk < x) { - k1 = k; p1 = pk; - } else { - k2 = k; p2 = pk; - } - } - return(k2); -} - -double -poissonpmf(long k, double L) -{ - double dx, dp; - - checkpoisson(L); - dx = k; - if (L == 0.0) { - dp = (k == 0) ? 1.0 : 0.0; - } else { - if (dx < 0.0) { - dp = 0.0; - } else { - dp = exp(dx * log(L) - L - gamma(dx + 1.0)); - } - } - return(dp); -} - -/* poisson random generator from Numerical Recipes */ -long poissonrand(double xm) -{ - static double sqrt2xm, logxm, expxm, g, oldxm = -1.0; - double t, y; - long k; - - checkpoisson(xm); - if (xm < 12.0) { - if (xm != oldxm) { expxm = exp(-xm); oldxm = xm; } - k = -1; - t = 1.0; - do { - k++; - t *= uni(); - } while (t > expxm); - } else { - if (xm != oldxm) { - oldxm = xm; - sqrt2xm = sqrt(2.0 * xm); - logxm = log(xm); - g = xm * logxm - gamma(xm + 1.0); - } - do { - do { - y = tan(PI * uni()); - k = floor(sqrt2xm * y + xm); - } while (k < 0); - t = 0.9 * (1.0 + y * y) * exp(k * logxm - gamma(k + 1.0) - g); - } while (uni() > t); - } - return (k); -} - -/*****************************************************************************/ -/** **/ -/** Binomial Distribution **/ -/** **/ -/*****************************************************************************/ - -static double -binomial_cdf(long k, long n, double p) -{ - double da, db, dp; - long ia, ib; - - if (k < 0) { - dp = 0.0; - } else { - if (k >= n) { - dp = 1.0; - } else { - if (p == 0.0) { - dp = (k < 0) ? 0.0 : 1.0; - } else { - if (p == 1.0) { - dp = (k < n) ? 0.0 : 1.0; - } else { - da = k + 1.0; - db = n - k; - ia = floor(da); - ib = floor(db); - betabase(&p, &da, &db, &ia, &ib, &dp); - dp = 1.0 - dp; - } - } - } - } - return(dp); -} - -double -binomialcdf(double k, long n, double p) -{ - checkbinomial(n, p); - return(binomial_cdf((long) floor(k), n, p)); -} - -long -binomialquant(double x, long n, double p) -{ - long k, k1, k2, del, ia; - double m, s, p1, p2, pk; - - checkbinomial(n, p); - checkprob(x, FALSE, FALSE); - - m = n * p; - s = sqrt(n * p * (1 - p)); - del = max(1, (long) (0.2 * s)); - - if (x == 0.0) k = 0.0; - else if (x == 1.0) k = n; - else k = m + s * ppnd(x, &ia); - k1 = k; k2 = k; - - do { - k1 = k1 - del; k1 = max(0, k1); - p1 = binomial_cdf(k1, n, p); - } while (k1 > 0 && p1 > p); - if (k1 == 0 && p1 >= x) return(k1); - - do { - k2 = k2 + del; k2 = min(n, k2); - p2 = binomial_cdf(k2, n, p); - } while (k2 < n && p2 < x); - if (k2 == n && p2 <= x) return(k2); - - while (k2 - k1 > 1) { - k = (k1 + k2) / 2; - pk = binomial_cdf(k, n, p); - if (pk < x) { k1 = k; p1 = pk; } - else { k2 = k; p2 = pk; } - } - return(k2); -} - -double binomialpmf(long k, long n, double p) -{ - double dx, dp; - - checkbinomial(n, p); - dx = k; - if (p == 0.0) dp = (k == 0) ? 1.0 : 0.0; - else if (p == 1.0) dp = (k == n) ? 1.0 : 0.0; - else if (dx < 0.0 || dx > n) dp = 0.0; - else dp = exp(gamma(n + 1.0) - gamma(dx + 1.0) - gamma(n - dx + 1.0) - + dx * log(p) + (n - dx) * log(1.0 - p)); - return(dp); -} - -/* binomial random generator from Numerical Recipes */ -long binomialrand(long n, double pp) -{ - long j, k; - static long nold = -1; - double am, em, g, p, sq, t, y; - static double pold = -1.0, pc, plog, pclog, en, oldg; - - checkbinomial(n, pp); - - p = (pp <= 0.5) ? pp : 1.0 - pp; - - am = n * p; - if (p == 0.0) k = 0; - else if (p == 1.0) k = n; - else if (n < 50) { - k = 0; - for (j = 0; j < n; j++) if (uni() < p) k++; - } - else if (am < 1.0) { - g = exp(-am); - t = 1.0; - k = -1; - do { - k++; - t *= uni(); - } while (t > g); - if (k > n) k = n; - } - else { - if (n != nold) { - en = n; - oldg = gamma(en + 1.0); - nold = n; - } - if (p != pold) { - pc = 1.0 - p; - plog = log(p); - pclog = log(pc); - pold = p; - } - sq = sqrt(2.0 * am * pc); - do { - do { - y = tan(PI * uni()); - em = sq * y + am; - } while (em < 0.0 || em >= en + 1.0); - em = floor(em); - t = 1.2 * sq * (1.0 + y * y) - * exp(oldg - gamma(em + 1.0) - gamma(en - em + 1.0) - + em * plog + (en - em) * pclog); - } while (uni() > t); - k = em; - } - if (p != pp) k = n - k; - return(k); -} diff --git a/lib/cffi-glue.c b/lib/cffi-glue.c deleted file mode 100644 index 12f2065..0000000 --- a/lib/cffi-glue.c +++ /dev/null @@ -1,477 +0,0 @@ -/* -#ifdef __APPLE__ -#include -#else -#include -#endif -*/ -#include -#include - - - -#include - -#include "linalg.h" - -typedef char *PTR; - -extern double unirand(), gamma(); -extern double normalcdf(), normalquant(), normaldens(), normalrand(); -extern double bnormcdf(); -extern double cauchycdf(), cauchyquant(), cauchydens(), cauchyrand(); -extern double gammacdf(), gammaquant(), gammadens(), gammarand(); -extern double chisqcdf(), chisqquant(), chisqdens(), chisqrand(); -extern double betacdf(), betaquant(), betadens(), betarand(); -extern double tcdf(), tquant(), tdens(), trand(); -extern double fcdf(), fquant(), fdens(), frand(); -extern double poissoncdf(), poissonpmf(); -extern long poissonquant(), poissonrand(); -extern double binomialcdf(), binomialpmf(); -extern long binomialquant(), binomialrand(); -extern double rcondest_front(); - -void xlfail(char *); - -extern int chol_decomp_front(); -extern int lu_decomp_front(); -extern int lu_solve_front(); -extern int lu_inverse_front(); -extern int sv_decomp_front(); -extern int qr_decomp_front(); -extern int make_rotation_front(); -extern int eigen_front(); -extern int range_to_rseq(); -extern int spline_front(); -extern int kernel_dens_front(); -extern int kernel_smooth_front(); -extern int base_lowess_front(); -extern int fft_front(); -extern int numgrad_front(); -extern int numhess_front(); -extern int minfo_maximize(); - -/***********************************************************************/ -/**** Basic Utilities ****/ -/***********************************************************************/ - -/***********************************************************************/ -/** Callback Support Functions **/ -/***********************************************************************/ - -static long ccl_integer_value; -static double ccl_double_value; -static PTR ccl_ptr_value; -void ccl_store_integer(long x) { ccl_integer_value = x; } -void ccl_store_double(double x) { ccl_double_value = x; } -void ccl_store_ptr(PTR x) { ccl_ptr_value = x; } - -/***************************************************************************/ -/** Lisp-Managed Calloc/Free **/ -/***************************************************************************/ - -#ifdef DODO -static void (*new_ptr)(); -static void (*free_ptr)(); - -register_new_ptr(f) void (*f)(); { new_ptr = f; } -register_free_ptr(f) void (*f)(); { free_ptr = f; } - -char* calloc(size_t n, size_t m) -{ - size_t i, N = n * m; - char *p; - - (*new_ptr)(N); - p = (char *) ccl_ptr_value; - for (i = 0; i < N; i++) p[i] = 0; - return(p); -} - -malloc() { xlfail("malloc not available yet"); } -realloc() { xlfail("realloc not available yet"); } - -void free(char *p) -{ - (*free_ptr)(p); -} -#endif /* DODO*/ - -/***************************************************************************/ -/** **/ -/** Storage Allocation Functions **/ -/** **/ -/***************************************************************************/ - -PTR -la_base_allocate(size_t n, size_t m) -{ - char *p = calloc(n, m); - if (p == nil) xlfail("allocation failed"); - return((PTR) p); -} - -long -la_base_free_alloc(PTR p) -{ - if (p) free((char *) p); - return(0); -} - -size_t -la_mode_size(int mode) -{ - switch (mode) { - case IN: return(sizeof(long)); - case RE: return(sizeof(double)); - case CX: return(sizeof(Complex)); - } - return(0); -} - -/***************************************************************************/ -/** **/ -/** Callbacks for Internal Storage **/ -/** **/ -/***************************************************************************/ - -int (*ccl_la_allocate)(), (*ccl_la_free_alloc)(); - -void register_la_allocate(f) int (*f)(); { ccl_la_allocate = f; } -void register_la_free_alloc(f) int (*f)(); { ccl_la_free_alloc = f; } - -PTR -la_allocate(size_t n, size_t m) -{ - (*ccl_la_allocate)(n, m); - return(ccl_ptr_value); -} - -void -la_free_alloc(PTR p) -{ - (*ccl_la_free_alloc)(p); -} - -/***************************************************************************/ -/** **/ -/** Storage Access Functions **/ -/** **/ -/***************************************************************************/ - -long -la_get_integer(PTR p, size_t i) -{ - return(*(((long *) p) + i)); -} - -double -la_get_double(PTR p, size_t i) -{ - return(*(((double *) p) + i)); -} - -double -la_get_complex_real(PTR p, size_t i) -{ - Complex *c = ((Complex *) p) + i; - return(c->real); -} - -double -la_get_complex_imag(PTR p, size_t i) -{ - Complex *c = ((Complex *) p) + i; - return(c->imag); -} - -PTR -la_get_pointer(PTR p, size_t i) -{ - return(*(((PTR *) p) + i)); -} - -/***************************************************************************/ -/** **/ -/** Storage Mutation Functions **/ -/** **/ -/***************************************************************************/ - -int -la_put_integer(PTR p, size_t i, long x) -{ - *(((long *) p) + i) = x; - return(0); -} - -int la_put_double(PTR p, size_t i, double x) -{ - *(((double *) p) + i) = x; - return(0); -} - -int -la_put_complex(PTR p, size_t i, double x, double y) -{ - Complex *c = ((Complex *) p) + i; - c->real = x; - c->imag = y; - return(0); -} - -int -la_put_pointer(PTR p, size_t i, PTR x) -{ - *(((PTR *) p) + i) = x; - return(0); -} - -/***********************************************************************/ -/** **/ -/** XLISP Internal Error Message Emulation **/ -/** **/ -/***********************************************************************/ - -char buf[1000]; - -static int (*ccl_set_buf_char_fptr)(); -void register_set_buf_char(f) int (*f)(); { ccl_set_buf_char_fptr = f; } -void set_buf_char(int n, int c) { (*ccl_set_buf_char_fptr)(n, c); } - -static int (*ccl_print_buffer)(); -void register_print_buffer(f) int (*f)(); { ccl_print_buffer = f; } -void print_buffer(int n, int m) { (*ccl_print_buffer)(n, m); } - -static int bufpos = 0; - -static void -resetbuf() -{ - bufpos = 0; -} - -static void -prbuf(char *s) -{ - size_t i, n; - - n = strlen(s); - for (i = 0; i 5) { - wa[i1 - 1] = wa[i - 1]; - wa[i1] = wa[i]; - } - } - l1 = l2; - } - return 0; -} - -static int pass_(int *nac, long *ido, long *ip, long *l1, long *idl1, - double *cc, double *c1, double *c2, double *ch, - double *ch2, double *wa, int *isign) -{ - /* System generated locals */ - long cc_dim1, cc_dim2, cc_offset, ch_dim1, ch_dim2, ch_offset, - ch2_offset, c1_dim1, c1_dim2, c1_offset, ch2_dim1, c2_offset, - c2_dim1, i_1, i_2, i_3; - - /* Local variables */ - int i, j, k, l, ik; - long idij, idlj, idot, ipph, lc, jc, idj, idl, inc, idp; - double wai, war; - long ipp2; - - /* Parameter adjustments */ - cc_dim1 = *ido; - cc_dim2 = *ip; - cc_offset = cc_dim1 * (cc_dim2 + 1) + 1; - cc -= cc_offset; - c1_dim1 = *ido; - c1_dim2 = *l1; - c1_offset = c1_dim1 * (c1_dim2 + 1) + 1; - c1 -= c1_offset; - c2_dim1 = *idl1; - c2_offset = c2_dim1 + 1; - c2 -= c2_offset; - ch_dim1 = *ido; - ch_dim2 = *l1; - ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; - ch -= ch_offset; - ch2_dim1 = *idl1; - ch2_offset = ch2_dim1 + 1; - ch2 -= ch2_offset; - --wa; - - /* Function Body */ - idot = *ido / 2; - ipp2 = *ip + 2; - ipph = (*ip + 1) / 2; - idp = *ip * *ido; - - if (*ido >= *l1) { - i_1 = ipph; - for (j = 2; j <= i_1; ++j) { - jc = ipp2 - j; - i_2 = *l1; - for (k = 1; k <= i_2; ++k) { - i_3 = *ido; - for (i = 1; i <= i_3; ++i) { - ch[i + (k + j * ch_dim2) * ch_dim1] = - cc[i + (j + k * cc_dim2) * cc_dim1] - + cc[i + (jc + k * cc_dim2) * cc_dim1]; - ch[i + (k + jc * ch_dim2) * ch_dim1] = - cc[i + (j + k * cc_dim2) * cc_dim1] - - cc[i + (jc + k * cc_dim2) * cc_dim1]; - } - } - } - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - i_2 = *ido; - for (i = 1; i <= i_2; ++i) { - ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1]; - } - } - } - else { - i_1 = ipph; - for (j = 2; j <= i_1; ++j) { - jc = ipp2 - j; - i_2 = *ido; - for (i = 1; i <= i_2; ++i) { - i_3 = *l1; - for (k = 1; k <= i_3; ++k) { - ch[i + (k + j * ch_dim2) * ch_dim1] = - cc[i + (j + k * cc_dim2) * cc_dim1] - + cc[i + (jc + k * cc_dim2) * cc_dim1]; - ch[i + (k + jc * ch_dim2) * ch_dim1] = - cc[i + (j + k * cc_dim2) * cc_dim1] - - cc[i + (jc + k * cc_dim2) * cc_dim1]; - } - } - } - i_1 = *ido; - for (i = 1; i <= i_1; ++i) { - i_2 = *l1; - for (k = 1; k <= i_2; ++k) { - ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * cc_dim2 + 1) * cc_dim1]; - } - } - } - idl = 2 - *ido; - inc = 0; - i_1 = ipph; - for (l = 2; l <= i_1; ++l) { - lc = ipp2 - l; - idl += *ido; - i_2 = *idl1; - for (ik = 1; ik <= i_2; ++ik) { - c2[ik + l * c2_dim1] = ch2[ik + ch2_dim1] + wa[idl - 1] - * ch2[ik + (ch2_dim1 << 1)]; - c2[ik + lc * c2_dim1] = - *isign * wa[idl] * ch2[ik + *ip * ch2_dim1]; - } - idlj = idl; - inc += *ido; - i_2 = ipph; - for (j = 3; j <= i_2; ++j) { - jc = ipp2 - j; - idlj += inc; - if (idlj > idp) { - idlj -= idp; - } - war = wa[idlj - 1]; - wai = wa[idlj]; - i_3 = *idl1; - for (ik = 1; ik <= i_3; ++ik) { - c2[ik + l * c2_dim1] += war * ch2[ik + j * ch2_dim1]; - c2[ik + lc * c2_dim1] -= *isign * wai * ch2[ik + jc * ch2_dim1]; - } - } - } - i_1 = ipph; - for (j = 2; j <= i_1; ++j) { - i_2 = *idl1; - for (ik = 1; ik <= i_2; ++ik) { - ch2[ik + ch2_dim1] += ch2[ik + j * ch2_dim1]; - } - } - i_1 = ipph; - for (j = 2; j <= i_1; ++j) { - jc = ipp2 - j; - i_2 = *idl1; - for (ik = 2; ik <= i_2; ik += 2) { - ch2[ik - 1 + j * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - - c2[ik + jc * c2_dim1]; - ch2[ik - 1 + jc * ch2_dim1] = c2[ik - 1 + j * c2_dim1] - + c2[ik + jc * c2_dim1]; - ch2[ik + j * ch2_dim1] = c2[ik + j * c2_dim1] - + c2[ik - 1 + jc * c2_dim1]; - ch2[ik + jc * ch2_dim1] = c2[ik + j * c2_dim1] - - c2[ik - 1 + jc * c2_dim1]; - } - } - *nac = 1; - if (*ido != 2) { - *nac = 0; - i_1 = *idl1; - for (ik = 1; ik <= i_1; ++ik) { - c2[ik + c2_dim1] = ch2[ik + ch2_dim1]; - } - i_1 = *ip; - for (j = 2; j <= i_1; ++j) { - i_2 = *l1; - for (k = 1; k <= i_2; ++k) { - c1[(k + j * c1_dim2) * c1_dim1 + 1] = - ch[(k + j * ch_dim2) * ch_dim1 + 1]; - c1[(k + j * c1_dim2) * c1_dim1 + 2] = - ch[(k + j * ch_dim2) * ch_dim1 + 2]; - } - } - if (idot <= *l1) { - idij = 0; - i_1 = *ip; - for (j = 2; j <= i_1; ++j) { - idij += 2; - i_2 = *ido; - for (i = 4; i <= i_2; i += 2) { - idij += 2; - i_3 = *l1; - for (k = 1; k <= i_3; ++k) { - c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] - * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] + *isign * wa[idij] - * ch[i + (k + j * ch_dim2) * ch_dim1]; - c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] - * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij] - * ch[i - 1 + (k + j * ch_dim2) * ch_dim1]; - } - } - } - } - else { - idj = 2 - *ido; - i_1 = *ip; - for (j = 2; j <= i_1; ++j) { - idj += *ido; - i_2 = *l1; - for (k = 1; k <= i_2; ++k) { - idij = idj; - i_3 = *ido; - for (i = 4; i <= i_3; i += 2) { - idij += 2; - c1[i - 1 + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] - * ch[i - 1 + (k + j * ch_dim2) * ch_dim1] - + *isign * wa[idij] * ch[i + (k + j * ch_dim2) * ch_dim1]; - c1[i + (k + j * c1_dim2) * c1_dim1] = wa[idij - 1] - * ch[i + (k + j * ch_dim2) * ch_dim1] - *isign * wa[idij] - * ch[i - 1 + (k + j * ch_dim2) * ch_dim1]; - } - } - } - } - } - return 0; -} - -static int pass2_(ido, l1, cc, ch, wa1, isign) - int *ido, *l1; - double *cc, *ch, *wa1; - int *isign; -{ - /* System generated locals */ - int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2; - - /* Local variables */ - int i, k; - double ti2, tr2; - - /* Parameter adjustments */ - cc_dim1 = *ido; - cc_offset = cc_dim1 * 3 + 1; - cc -= cc_offset; - ch_dim1 = *ido; - ch_dim2 = *l1; - ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; - ch -= ch_offset; - --wa1; - - /* Function Body */ - if (*ido <= 2) { - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - ch[(k + ch_dim2) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] - + cc[((k << 1) + 2) * cc_dim1 + 1]; - ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cc[((k << 1) + 1) * cc_dim1 + 1] - - cc[((k << 1) + 2) * cc_dim1 + 1]; - ch[(k + ch_dim2) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] - + cc[((k << 1) + 2) * cc_dim1 + 2]; - ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = cc[((k << 1) + 1) * cc_dim1 + 2] - - cc[((k << 1) + 2) * cc_dim1 + 2]; - } - } - else { - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - i_2 = *ido; - for (i = 2; i <= i_2; i += 2) { - ch[i - 1 + (k + ch_dim2) * ch_dim1] - = cc[i - 1 + ((k << 1) + 1) * cc_dim1] - + cc[i - 1 + ((k << 1) + 2) * cc_dim1]; - tr2 = cc[i - 1 + ((k << 1) + 1) * cc_dim1] - - cc[i - 1 + ((k << 1) + 2) * cc_dim1]; - ch[i + (k + ch_dim2) * ch_dim1] = cc[i + ((k << 1) + 1) * cc_dim1] - + cc[i + ((k << 1) + 2) * cc_dim1]; - ti2 = cc[i + ((k << 1) + 1) * cc_dim1] - - cc[i + ((k << 1) + 2) * cc_dim1]; - ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ti2 - - *isign * wa1[i] * tr2; - ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * tr2 - + *isign * wa1[i] * ti2; - } - } - } - return 0; -} - -static int pass3_(ido, l1, cc, ch, wa1, wa2, isign) - int *ido, *l1; - double *cc, *ch, *wa1, *wa2; - int *isign; -{ - /* System generated locals */ - int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2; - - /* Local variables */ - double taui, taur; - int i, k; - double ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2; - - /* Parameter adjustments */ - cc_dim1 = *ido; - cc_offset = (cc_dim1 << 2) + 1; - cc -= cc_offset; - ch_dim1 = *ido; - ch_dim2 = *l1; - ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; - ch -= ch_offset; - --wa1; - --wa2; - - /* Function Body */ - taur = -.5; - taui = -(*isign) * .866025403784439; - if (*ido == 2) { - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - tr2 = cc[(k * 3 + 2) * cc_dim1 + 1] + cc[(k * 3 + 3) * cc_dim1 + 1]; - cr2 = cc[(k * 3 + 1) * cc_dim1 + 1] + taur * tr2; - ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 3 + 1) * cc_dim1 + 1] + tr2; - ti2 = cc[(k * 3 + 2) * cc_dim1 + 2] + cc[(k * 3 + 3) * cc_dim1 + 2]; - ci2 = cc[(k * 3 + 1) * cc_dim1 + 2] + taur * ti2; - ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 3 + 1) * cc_dim1 + 2] + ti2; - cr3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 1] - - cc[(k * 3 + 3) * cc_dim1 + 1]); - ci3 = taui * (cc[(k * 3 + 2) * cc_dim1 + 2] - - cc[(k * 3 + 3) * cc_dim1 + 2]); - ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci3; - ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr2 + ci3; - ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr3; - ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci2 - cr3; - } - } - else { - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - i_2 = *ido; - for (i = 2; i <= i_2; i += 2) { - tr2 = cc[i - 1 + (k * 3 + 2) * cc_dim1] - + cc[i - 1 + (k * 3 + 3) * cc_dim1]; - cr2 = cc[i - 1 + (k * 3 + 1) * cc_dim1] + taur * tr2; - ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 3 + 1) - * cc_dim1] + tr2; - ti2 = cc[i + (k * 3 + 2) * cc_dim1] + cc[i + (k * 3 + 3) * cc_dim1]; - ci2 = cc[i + (k * 3 + 1) * cc_dim1] + taur * ti2; - ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 3 + 1) * cc_dim1] + ti2; - cr3 = taui * (cc[i - 1 + (k * 3 + 2) * cc_dim1] - - cc[i - 1 + (k * 3 + 3) * cc_dim1]); - ci3 = taui * (cc[i + (k * 3 + 2) * cc_dim1] - - cc[i + (k * 3 + 3) * cc_dim1]); - dr2 = cr2 - ci3; - dr3 = cr2 + ci3; - di2 = ci2 + cr3; - di3 = ci2 - cr3; - ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 - - *isign * wa1[i] * dr2; - ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 - + *isign * wa1[i] * di2; - ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 - - *isign * wa2[i] * dr3; - ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 - + *isign * wa2[i] * di3; - } - } - } - return 0; -} - -static int pass4_(ido, l1, cc, ch, wa1, wa2, wa3, isign) - int *ido, *l1; - double *cc, *ch, *wa1, *wa2, *wa3; - int *isign; -{ - /* System generated locals */ - int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2; - - /* Local variables */ - int i, k; - double ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; - - /* Parameter adjustments */ - cc_dim1 = *ido; - cc_offset = cc_dim1 * 5 + 1; - cc -= cc_offset; - ch_dim1 = *ido; - ch_dim2 = *l1; - ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; - ch -= ch_offset; - --wa1; - --wa2; - --wa3; - - /* Function Body */ - if (*ido == 2) { - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - ti1 = cc[((k << 2) + 1) * cc_dim1 + 2] - - cc[((k << 2) + 3) * cc_dim1 + 2]; - ti2 = cc[((k << 2) + 1) * cc_dim1 + 2] - + cc[((k << 2) + 3) * cc_dim1 + 2]; - tr4 = *isign * (cc[((k << 2) + 2) * cc_dim1 + 2] - - cc[((k << 2) + 4) * cc_dim1 + 2]); - ti3 = cc[((k << 2) + 2) * cc_dim1 + 2] - + cc[((k << 2) + 4) * cc_dim1 + 2]; - tr1 = cc[((k << 2) + 1) * cc_dim1 + 1] - - cc[((k << 2) + 3) * cc_dim1 + 1]; - tr2 = cc[((k << 2) + 1) * cc_dim1 + 1] - + cc[((k << 2) + 3) * cc_dim1 + 1]; - ti4 = *isign * (cc[((k << 2) + 4) * cc_dim1 + 1] - - cc[((k << 2) + 2) * cc_dim1 + 1]); - tr3 = cc[((k << 2) + 2) * cc_dim1 + 1] - + cc[((k << 2) + 4) * cc_dim1 + 1]; - ch[(k + ch_dim2) * ch_dim1 + 1] = tr2 + tr3; - ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = tr2 - tr3; - ch[(k + ch_dim2) * ch_dim1 + 2] = ti2 + ti3; - ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ti2 - ti3; - ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = tr1 + tr4; - ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = tr1 - tr4; - ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ti1 + ti4; - ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ti1 - ti4; - } - } - else { - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - i_2 = *ido; - for (i = 2; i <= i_2; i += 2) { - ti1 = cc[i + ((k << 2) + 1) * cc_dim1] - - cc[i + ((k << 2) + 3) * cc_dim1]; - ti2 = cc[i + ((k << 2) + 1) * cc_dim1] - + cc[i + ((k << 2) + 3) * cc_dim1]; - ti3 = cc[i + ((k << 2) + 2) * cc_dim1] - + cc[i + ((k << 2) + 4) * cc_dim1]; - tr4 = *isign * (cc[i + ((k << 2) + 2) * cc_dim1] - - cc[i + ((k << 2) + 4) * cc_dim1]); - tr1 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] - - cc[i - 1 + ((k << 2) + 3) * cc_dim1]; - tr2 = cc[i - 1 + ((k << 2) + 1) * cc_dim1] - + cc[i - 1 + ((k << 2) + 3) * cc_dim1]; - ti4 = *isign * (cc[i - 1 + ((k << 2) + 4) * cc_dim1] - - cc[i - 1 + ((k << 2) + 2) * cc_dim1]); - tr3 = cc[i - 1 + ((k << 2) + 2) * cc_dim1] - + cc[i - 1 + ((k << 2) + 4) * cc_dim1]; - ch[i - 1 + (k + ch_dim2) * ch_dim1] = tr2 + tr3; - cr3 = tr2 - tr3; - ch[i + (k + ch_dim2) * ch_dim1] = ti2 + ti3; - ci3 = ti2 - ti3; - cr2 = tr1 + tr4; - cr4 = tr1 - tr4; - ci2 = ti1 + ti4; - ci4 = ti1 - ti4; - ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * cr2 - + *isign * wa1[i] * ci2; - ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * ci2 - - *isign * wa1[i] * cr2; - ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * cr3 - + *isign * wa2[i] * ci3; - ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * ci3 - - *isign * wa2[i] * cr3; - ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * cr4 - + *isign * wa3[i] * ci4; - ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * ci4 - - *isign * wa3[i] * cr4; - } - } - } - return 0; -} - -static int pass5_(ido, l1, cc, ch, wa1, wa2, wa3, wa4, isign) - int *ido, *l1; - double *cc, *ch, *wa1, *wa2, *wa3, *wa4; - int *isign; -{ - /* System generated locals */ - int cc_dim1, cc_offset, ch_dim1, ch_dim2, ch_offset, i_1, i_2; - - /* Local variables */ - int i, k; - double ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3, - ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5, ti11, ti12, tr11, - tr12; - - /* Parameter adjustments */ - cc_dim1 = *ido; - cc_offset = cc_dim1 * 6 + 1; - cc -= cc_offset; - ch_dim1 = *ido; - ch_dim2 = *l1; - ch_offset = ch_dim1 * (ch_dim2 + 1) + 1; - ch -= ch_offset; - --wa1; - --wa2; - --wa3; - --wa4; - - /* Function Body */ - tr11 = .309016994374947; - ti11 = -(*isign) * .951056516295154; - tr12 = -.809016994374947; - ti12 = -(*isign) * .587785252292473; - if (*ido == 2) { - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - ti5 = cc[(k * 5 + 2) * cc_dim1 + 2] - cc[(k * 5 + 5) * cc_dim1 + 2]; - ti2 = cc[(k * 5 + 2) * cc_dim1 + 2] + cc[(k * 5 + 5) * cc_dim1 + 2]; - ti4 = cc[(k * 5 + 3) * cc_dim1 + 2] - cc[(k * 5 + 4) * cc_dim1 + 2]; - ti3 = cc[(k * 5 + 3) * cc_dim1 + 2] + cc[(k * 5 + 4) * cc_dim1 + 2]; - tr5 = cc[(k * 5 + 2) * cc_dim1 + 1] - cc[(k * 5 + 5) * cc_dim1 + 1]; - tr2 = cc[(k * 5 + 2) * cc_dim1 + 1] + cc[(k * 5 + 5) * cc_dim1 + 1]; - tr4 = cc[(k * 5 + 3) * cc_dim1 + 1] - cc[(k * 5 + 4) * cc_dim1 + 1]; - tr3 = cc[(k * 5 + 3) * cc_dim1 + 1] + cc[(k * 5 + 4) * cc_dim1 + 1]; - ch[(k + ch_dim2) * ch_dim1 + 1] = cc[(k * 5 + 1) * cc_dim1 + 1] - + tr2 + tr3; - ch[(k + ch_dim2) * ch_dim1 + 2] = cc[(k * 5 + 1) * cc_dim1 + 2] - + ti2 + ti3; - cr2 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr11 * tr2 + tr12 * tr3; - ci2 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr11 * ti2 + tr12 * ti3; - cr3 = cc[(k * 5 + 1) * cc_dim1 + 1] + tr12 * tr2 + tr11 * tr3; - ci3 = cc[(k * 5 + 1) * cc_dim1 + 2] + tr12 * ti2 + tr11 * ti3; - cr5 = ti11 * tr5 + ti12 * tr4; - ci5 = ti11 * ti5 + ti12 * ti4; - cr4 = ti12 * tr5 - ti11 * tr4; - ci4 = ti12 * ti5 - ti11 * ti4; - ch[(k + (ch_dim2 << 1)) * ch_dim1 + 1] = cr2 - ci5; - ch[(k + ch_dim2 * 5) * ch_dim1 + 1] = cr2 + ci5; - ch[(k + (ch_dim2 << 1)) * ch_dim1 + 2] = ci2 + cr5; - ch[(k + ch_dim2 * 3) * ch_dim1 + 2] = ci3 + cr4; - ch[(k + ch_dim2 * 3) * ch_dim1 + 1] = cr3 - ci4; - ch[(k + (ch_dim2 << 2)) * ch_dim1 + 1] = cr3 + ci4; - ch[(k + (ch_dim2 << 2)) * ch_dim1 + 2] = ci3 - cr4; - ch[(k + ch_dim2 * 5) * ch_dim1 + 2] = ci2 - cr5; - } - } - else { - i_1 = *l1; - for (k = 1; k <= i_1; ++k) { - i_2 = *ido; - for (i = 2; i <= i_2; i += 2) { - ti5 = cc[i + (k * 5 + 2) * cc_dim1] - cc[i + (k * 5 + 5) * cc_dim1]; - ti2 = cc[i + (k * 5 + 2) * cc_dim1] + cc[i + (k * 5 + 5) * cc_dim1]; - ti4 = cc[i + (k * 5 + 3) * cc_dim1] - cc[i + (k * 5 + 4) * cc_dim1]; - ti3 = cc[i + (k * 5 + 3) * cc_dim1] + cc[i + (k * 5 + 4) * cc_dim1]; - tr5 = cc[i - 1 + (k * 5 + 2) * cc_dim1] - - cc[i - 1 + (k * 5 + 5) * cc_dim1]; - tr2 = cc[i - 1 + (k * 5 + 2) * cc_dim1] - + cc[i - 1 + (k * 5 + 5) * cc_dim1]; - tr4 = cc[i - 1 + (k * 5 + 3) * cc_dim1] - - cc[i - 1 + (k * 5 + 4) * cc_dim1]; - tr3 = cc[i - 1 + (k * 5 + 3) * cc_dim1] - + cc[i - 1 + (k * 5 + 4) * cc_dim1]; - ch[i - 1 + (k + ch_dim2) * ch_dim1] = cc[i - 1 + (k * 5 + 1) - * cc_dim1] + tr2 + tr3; - ch[i + (k + ch_dim2) * ch_dim1] = cc[i + (k * 5 + 1) - * cc_dim1] + ti2 + ti3; - cr2 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr11 * tr2 + tr12 * tr3; - ci2 = cc[i + (k * 5 + 1) * cc_dim1] + tr11 * ti2 + tr12 * ti3; - - cr3 = cc[i - 1 + (k * 5 + 1) * cc_dim1] + tr12 * tr2 + tr11 * tr3; - ci3 = cc[i + (k * 5 + 1) * cc_dim1] + tr12 * ti2 + tr11 * ti3; - - cr5 = ti11 * tr5 + ti12 * tr4; - ci5 = ti11 * ti5 + ti12 * ti4; - cr4 = ti12 * tr5 - ti11 * tr4; - ci4 = ti12 * ti5 - ti11 * ti4; - dr3 = cr3 - ci4; - dr4 = cr3 + ci4; - di3 = ci3 + cr4; - di4 = ci3 - cr4; - dr5 = cr2 + ci5; - dr2 = cr2 - ci5; - di5 = ci2 - cr5; - di2 = ci2 + cr5; - ch[i - 1 + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * dr2 - + *isign * wa1[i] * di2; - ch[i + (k + (ch_dim2 << 1)) * ch_dim1] = wa1[i - 1] * di2 - - *isign * wa1[i] * dr2; - ch[i - 1 + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * dr3 - + *isign * wa2[i] * di3; - ch[i + (k + ch_dim2 * 3) * ch_dim1] = wa2[i - 1] * di3 - - *isign * wa2[i] * dr3; - ch[i - 1 + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * dr4 - + *isign * wa3[i] * di4; - ch[i + (k + (ch_dim2 << 2)) * ch_dim1] = wa3[i - 1] * di4 - - *isign * wa3[i] * dr4; - ch[i - 1 + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * dr5 - + *isign * wa4[i] * di5; - ch[i + (k + ch_dim2 * 5) * ch_dim1] = wa4[i - 1] * di5 - - *isign * wa4[i] * dr5; - } - } - } - return 0; -} - - -/*static */ -int -cfft1_(int *n, double *c, double *ch, double *wa, - double /* int??*/ *ifac, int *isign) -{ - /* System generated locals */ - long i_1; - - /* Local variables */ - long idot; - int i, k1, n2, na, nac; - long l1, l2, ido, idl1, iw, ix2, ix3, ix4; - long nf, ip; - - /* Parameter adjustments */ - --c; - --ch; - --wa; - --ifac; - - /* Function Body */ - nf = ifac[2]; - na = 0; - l1 = 1; - iw = 1; - i_1 = nf; - for (k1 = 1; k1 <= i_1; ++k1) { - ip = ifac[k1 + 2]; - l2 = ip * l1; - ido = *n / l2; - idot = ido + ido; - idl1 = idot * l1; - if (ip == 4) { - ix2 = iw + idot; - ix3 = ix2 + idot; - if (na == 0) { - pass4_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], isign); - } - else { - pass4_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], isign); - } - na = 1 - na; - } - else if (ip == 2) { - if (na == 0) { - pass2_(&idot, &l1, &c[1], &ch[1], &wa[iw], isign); - } - else { - pass2_(&idot, &l1, &ch[1], &c[1], &wa[iw], isign); - } - na = 1 - na; - } - else if (ip == 3) { - ix2 = iw + idot; - if (na == 0) { - pass3_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], isign); - } - else { - pass3_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], isign); - } - na = 1 - na; - } - else if (ip == 5) { - ix2 = iw + idot; - ix3 = ix2 + idot; - ix4 = ix3 + idot; - if (na == 0) { - pass5_(&idot, &l1, &c[1], &ch[1], &wa[iw], &wa[ix2], &wa[ix3], - &wa[ix4], isign); - } - else { - pass5_(&idot, &l1, &ch[1], &c[1], &wa[iw], &wa[ix2], &wa[ix3], - &wa[ix4], isign); - } - na = 1 - na; - } - else { - if (na == 0) { - pass_(&nac, &idot, &ip, &l1, &idl1, &c[1], &c[1], &c[1], &ch[1], - &ch[1], &wa[iw], isign); - } - else { - pass_(&nac, &idot, &ip, &l1, &idl1, &ch[1], &ch[1], &ch[1], &c[1], - &c[1], &wa[iw], isign); - } - if (nac != 0) { - na = 1 - na; - } - } - l1 = l2; - iw += (ip - 1) * idot; - } - if (na != 0) { - n2 = *n + *n; - i_1 = n2; - for (i = 1; i <= i_1; ++i) { - c[i] = ch[i]; - } - } - return 0; -} - - -/* - Public Routine -*/ - -/* - * cfft computes the forward or backward complex discrete fourier transform. - * - * Input Parameters: - * - * n The length of the complex sequence c. The method is - * more efficient when n is the product of small primes. - * - * c A complex array of length n which contains the sequence - * - * wsave a real work array which must be dimensioned at least 4n+15 - * in the program that calls cfft. - * isign 1 for transform, -1 for inverse transform. - * A call of cfft with isign = 1 followed by a call of cfft with - * isign = -1 will multiply the sequence by n. - * - * Output Parameters: - * - * c For j=1,...,n - * - * c(j)=the sum from k=1,...,n of - * - * c(k)*exp(-isign*i*(j-1)*(k-1)*2*pi/n) - * - * where i=sqrt(-1) - */ - -int cfft(int n, double *c, double *wsave, int isign) -{ - int iw1, iw2; - - /* Parameter adjustments */ - --c; - --wsave; - - /* Function Body */ - if (n != 1) { - iw1 = n + n + 1; - iw2 = iw1 + n + n; - cffti1_(&n, &wsave[iw1], &wsave[iw2]); - cfft1_(&n, &c[1], &wsave[1], &wsave[iw1], &wsave[iw2], &isign); - } - return 0; -} diff --git a/lib/cholesky.c b/lib/cholesky.c deleted file mode 100644 index 2829d46..0000000 --- a/lib/cholesky.c +++ /dev/null @@ -1,79 +0,0 @@ -/* choldecomp - Cholesky decomposition routines. */ -/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */ -/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */ -/* You may give out copies of this software; for conditions see the */ -/* file COPYING included with this distribution. */ - -#include "linalg.h" - -/* -choldecomp(a, n) - RMatrix a; - int n; -{ - int i, j, k; - - for (i = 0; i < n; i++) { - for (k = 0; k < i; k++) - a[i][i] -= a[k][i] * a[k][i]; - a[i][i] = (a[i][i] > 0) ? sqrt(a[i][i]) : 0.0; - for (j = i + 1; j < n; j++) { - for (k = 0; k < i; k++) a[i][j] -= a[k][i] * a[k][j]; - a[i][j] = (a[i][i] > 0.0) ? a[i][j] / a[i][i] : 0.0; - } - } - for (i = 0; i < n; i++) - for (j = 0; j < i; j++) - a[i][j] = 0.0; -} -*/ - -static double -Max(double a, double b) -{ - return(a > b ? a : b); -} - -void -choldecomp(double **a, size_t n, double maxoffl, double *maxadd) -{ - double minl, minljj, minl2; - size_t i, j, k; - - minl = pow(macheps(), 0.25) * maxoffl; - minl2 = 0.0; - - if (maxoffl == 0.0) { - for (i = 0; i < n; i++) - maxoffl = Max(fabs(a[i][i]), maxoffl); - maxoffl = sqrt(maxoffl); - minl2 = sqrt(macheps()) * maxoffl; - } - - *maxadd = 0.0; - for (j = 0; j < n; j++) { - for (i = 0; i < j; i++) a[j][j] -= a[j][i] * a[j][i]; - - minljj = 0.0; - - for (i = j + 1; i < n; i++) { - a[i][j] = a[j][i]; - for (k = 0; k < j; k++) a[i][j] -= a[i][k] * a[j][k]; - minljj = Max(fabs(a[i][j]), minljj); - } - - minljj = Max(minljj / maxoffl, minl); - - if (a[j][j] > minljj * minljj) a[j][j] = sqrt(a[j][j]); - else { - if (minljj < minl2) minljj = minl2; - *maxadd = Max(*maxadd, minljj * minljj - a[j][j]); - a[j][j] = minljj; - } - - for (i = j + 1; i < n; i++) a[i][j] /= a[j][j]; - } - - for (i = 0; i < n; i++) - for (j = i + 1; j < n; j++) a[i][j] = 0.0; -} diff --git a/lib/clib.make b/lib/clib.make deleted file mode 100644 index c7c8350..0000000 --- a/lib/clib.make +++ /dev/null @@ -1,134 +0,0 @@ -############################################################# -## Uncomment one of the following two groups for a generic -## application or one compiled to use 68020/68881 intructions. -## -## Generic application -## -CC = C -## -## MC68020/MC68881 application -## -#CC = C -mc68881 -elems881 -mc68020 -############################################################# -############################################################# -# File: clib.make -# Target: clib -# Sources: betabase.c -# bivnor.c -# cbayes.c -# cdists.c -# cfft.c -# cholesky.c -# clinalg.c -# complex.c -# derivatives.c -# functions.c -# gamln.c -# gammabase.c -# kernel.c -# linalgdata.c -# lowess.c -# ludecomp.c -# makerotation.c -# minimize.c -# nor.c -# ppnd.c -# qrdecomp.c -# rcondest.c -# splines.c -# studentbase.c -# svdecomp.c -# utils.c -# mclglue.c -# Created: Friday, August 17, 1990 8:16:25 AM - - -OBJECTS = 6 - betabase.c.o 6 - bivnor.c.o 6 - cbayes.c.o 6 - cdists.c.o 6 - cfft.c.o 6 - cholesky.c.o 6 - clinalg.c.o 6 - complex.c.o 6 - derivatives.c.o 6 - eigen.c.o 6 - functions.c.o 6 - gamln.c.o 6 - gammabase.c.o 6 - kernel.c.o 6 - linalgdata.c.o 6 - lowess.c.o 6 - ludecomp.c.o 6 - makerotation.c.o 6 - minimize.c.o 6 - nor.c.o 6 - ppnd.c.o 6 - qrdecomp.c.o 6 - rcondest.c.o 6 - splines.c.o 6 - studentbase.c.o 6 - svdecomp.c.o 6 - - -betabase.c.o D clib.make betabase.c - {CC} -s BETABASE betabase.c -bivnor.c.o D clib.make bivnor.c - {CC} -s BIVNOR bivnor.c -cbayes.c.o D clib.make cbayes.c - {CC} -s BAYES cbayes.c -cdists.c.o D clib.make cdists.c - {CC} -s CDISTS cdists.c -cfft.c.o D clib.make cfft.c - {CC} -s CFFT cfft.c -cholesky.c.o D clib.make cholesky.c - {CC} -s CHOLESKY cholesky.c -clinalg.c.o D clib.make clinalg.c - {CC} -s LINALG clinalg.c -complex.c.o D clib.make complex.c - {CC} -s COMPLEX complex.c -derivatives.c.o D clib.make derivatives.c - {CC} -s BAYES derivatives.c -eigen.c.o D clib.make eigen.c - {CC} -s EIGEN eigen.c -functions.c.o D clib.make functions.c - {CC} -s BAYES functions.c -gamln.c.o D clib.make gamln.c - {CC} -s GAMMA gamln.c -gammabase.c.o D clib.make gammabase.c - {CC} -s GAMMA gammabase.c -kernel.c.o D clib.make kernel.c - {CC} -s SMOOTH kernel.c -linalgdata.c.o D clib.make linalgdata.c - {CC} -s LINALG linalgdata.c -lowess.c.o D clib.make lowess.c - {CC} -s SMOOTH lowess.c -ludecomp.c.o D clib.make ludecomp.c - {CC} -s LUDECOMP ludecomp.c -makerotation.c.o D clib.make makerotation.c - {CC} -s MAKEROT makerotation.c -minimize.c.o D clib.make minimize.c - {CC} -s BAYES minimize.c -nor.c.o D clib.make nor.c - {CC} -s NORMAL nor.c -ppnd.c.o D clib.make ppnd.c - {CC} -s NORMAL ppnd.c -qrdecomp.c.o D clib.make qrdecomp.c - {CC} -d SCALE=XSSCALE -s QRDECOMP qrdecomp.c # needed besause of name conflict with SCALE??? -rcondest.c.o D clib.make rcondest.c - {CC} -s RCONDEST rcondest.c -splines.c.o D clib.make splines.c - {CC} -s SMOOTH splines.c -studentbase.c.o D clib.make studentbase.c - {CC} -s STUDENT studentbase.c -svdecomp.c.o D clib.make svdecomp.c - {CC} -s SVDECOMP svdecomp.c -mclglue.c.o D clib.make mclglue.c - {CC} -s LSGLUE mclglue.c - -clib.o DD clib.make {OBJECTS} mclglue.c.o ccldists.c.o - lib -o clib.o {OBJECTS} - -clib DD clib.make clib.o mclglue.c.o - diff --git a/lib/clinalg.c b/lib/clinalg.c deleted file mode 100644 index e097ffa..0000000 --- a/lib/clinalg.c +++ /dev/null @@ -1,218 +0,0 @@ -/* - C interface to basic linear algebra routines. */ -/* Copyright (c) 1990, by Luke Tierney */ - -#include "linalg.h" - -#ifdef INTPTR -typedef int PTR; -#else -typedef char *PTR; -#endif - -extern double rcondest(); - -extern void choldecomp(); -extern int crludcmp(); -extern int crlubksb(); -extern int svdcmp(); -extern void qrdecomp(); -extern void make_rotation(); -extern void eigen(); -extern void cfft(); -extern int lowess(); -extern int fit_spline(); -extern int kernel_smooth(); - - - -long -min (long x, long y) -{ - return((x < y) ? x : y); -} - -long -max (long x, long y) -{ - return((x > y) ? x : y); -} - -/************************************************************************/ -/** **/ -/** Machine Epsilon Determination **/ -/** **/ -/************************************************************************/ - -double macheps() -{ - static int calculated = FALSE; - static double epsilon = 1.0; - - if (! calculated) - while (1.0 + epsilon / 2.0 != 1.0) epsilon = epsilon / 2.0; - calculated = TRUE; - return(epsilon); -} - -/************************************************************************/ -/** **/ -/** Lisp Interfaces to Linear Algebra Routines **/ -/** **/ -/************************************************************************/ - -void -chol_decomp_front(PTR mat, size_t n, PTR dpars) -{ - double *dp = (double *) dpars; - choldecomp((double **) mat, n, *dp, dp + 1); -} - -int -lu_decomp_front(PTR mat, size_t n, PTR iv, int mode, PTR dp) -{ - return(crludcmp((char **) mat, n, (int *) iv, mode, (double *) dp)); -} - -int -lu_solve_front(PTR a, size_t n, PTR indx, PTR b, int mode) -{ - return(crlubksb((char **) a, n, (int *) indx, (char *) b, mode)); -} - -int -lu_inverse_front(PTR pmat, size_t n, PTR piv, PTR pv, int mode, PTR pinv) -{ - Matrix mat = (Matrix) pmat, inv = (Matrix) pinv; - IVector iv = (IVector) piv; - Vector v = (Vector) pv; - CMatrix cinv; - RMatrix rinv; - CVector cv; - RVector rv; - double d; - int i, j, singular; - - singular = crludcmp(mat, n, iv, mode, &d); - - if (! singular) { - rinv = (RMatrix) inv; - cinv = (CMatrix) inv; - rv = (RVector) v; - cv = (CVector) v; - - for (j = 0; j < n; j++) { - for (i = 0; i < n; i++) { - if (mode == RE) rv[i] = rinv[i][j]; - else cv[i] = cinv[i][j]; - } - - singular = singular || crlubksb(mat, n, iv, v, mode); - - for (i = 0; i < n; i++) { - if (mode == RE) rinv[i][j] = rv[i]; - else cinv[i][j] = cv[i]; - } - } - } - return(singular); -} - -int -sv_decomp_front(PTR mat, size_t m, size_t n, PTR w, PTR v) -{ - return(svdcmp((char **) mat, m, n, (char *) w, (char **) v)); -} - -void -qr_decomp_front(PTR mat, size_t m, size_t n, PTR v, PTR jpvt, int pivot) -{ - qrdecomp((char **) mat, m, n, (char **) v, (char *) jpvt, pivot); -} - -double -rcondest_front(PTR mat, size_t n) -{ - return(rcondest((char **) mat, n)); -} - -void -make_rotation_front(size_t n, PTR rot, PTR x, PTR y, int use_alpha, double alpha) -{ - make_rotation(n, (char **) rot, (char *) x, (char *) y, use_alpha, alpha); -} - -int -eigen_front(PTR a, size_t n, PTR w, PTR z, PTR fv1) -{ - int ierr; - - eigen(&n, &n, (char *) a, (char *) w, (char *) z, (char *) fv1, &ierr); - return(ierr); -} - -void -fft_front(size_t n, PTR x, PTR work, int isign) -{ - cfft(n, (char *) x, (char *) work, isign); -} - -int -base_lowess_front(PTR x, PTR y, size_t n, double f, - size_t nsteps, double delta, PTR ys, PTR rw, PTR res) -{ - return(lowess((char *) x, (char *) y, n, f, nsteps, delta, - (char *) ys, (char *) rw, (char *) res)); -} - -void -range_to_rseq(size_t n, PTR px, long ns, PTR pxs) -{ - size_t i; - double xmin, xmax, *x, *xs; - - x = (double *) px; - xs = (double *) pxs; - for (xmax = xmin = x[0], i = 1; i < n; i++) { - if (x[i] > xmax) xmax = x[i]; - if (x[i] < xmin) xmin = x[i]; - } - for (i = 0; i < ns; i++) - xs[i] = xmin + (xmax - xmin) * ((double) i) / ((double) (ns - 1)); -} - -int -spline_front(size_t n, PTR x, PTR y, size_t ns, - PTR xs, PTR ys, PTR work) -{ - return(fit_spline(n, (char *) x, (char *) y, - ns, (char *) xs, (char *) ys, (char *) work)); -} - -int -kernel_dens_front(PTR x, size_t n, double width, PTR xs, PTR ys, size_t ns, int code) -{ - int ktype; - - if (code == 0) ktype = 'U'; - else if (code == 1) ktype = 'T'; - else if (code == 2) ktype = 'G'; - else ktype = 'B'; - - return(kernel_smooth((char *) x, nil, n, width, nil, nil, - (char *) xs, (char *) ys, ns, ktype)); -} - -int -kernel_smooth_front(PTR x, PTR y, size_t n, double width, - PTR xs, PTR ys, size_t ns, int code) -{ - int ktype; - - if (code == 0) ktype = 'U'; - else if (code == 1) ktype = 'T'; - else if (code == 2) ktype = 'G'; - else ktype = 'B'; - - return(kernel_smooth((char *) x, (char *) y, n, width, nil, nil, - (char *) xs, (char *) ys, ns, ktype)); -} diff --git a/lib/cmpinclude.h b/lib/cmpinclude.h deleted file mode 100644 index f7945de..0000000 --- a/lib/cmpinclude.h +++ /dev/null @@ -1,695 +0,0 @@ - - -/* Begin for cmpinclude */ - - -/* #define SGC */ - - -/* End for cmpinclude */ -/* -(c) Copyright Taiichi Yuasa and Masami Hagiya, 1984. All rights reserved. -Copying of this file is authorized to users who have executed the true and -proper "License Agreement for Kyoto Common LISP" with SIGLISP. -*/ -#include -#include -#include -#define TRUE 1 -#define FALSE 0 -#ifdef SGC -#define FIRSTWORD short t; char s,m -#define SGC_TOUCH(x) x->d.m=0 -#else -#define FIRSTWORD short t; short m -#define SGC_TOUCH(x) -#endif -#define STSET(type,x,i,val) do{SGC_TOUCH(x);STREF(type,x,i) = (val);} while(0) -#ifndef VOL -#define VOL -#endif -#ifndef COM_LENG -#define COM_LENG -#endif -#ifndef CHAR_SIZE -#define CHAR_SIZE 8 -#endif -typedef int bool; -typedef int fixnum; -typedef float shortfloat; -typedef double longfloat; -typedef unsigned short fatchar; -#define SIGNED_CHAR(x) (((char ) -1) < (char )0 ? (char) x \ - : (x >= (1<<(CHAR_SIZE-1)) ? \ - x - (((int)(1<<(CHAR_SIZE-1))) << 1) \ - : (char ) x)) -typedef union lispunion *object; -typedef union int_object iobject; -union int_object {int i; object o;}; - -#define OBJNULL ((object)NULL) -struct fixnum_struct { - FIRSTWORD; - fixnum FIXVAL; -}; -#define fix(x) (x)->FIX.FIXVAL -#define SMALL_FIXNUM_LIMIT 1024 -extern struct fixnum_struct small_fixnum_table[COM_LENG]; -#define small_fixnum(i) (object)(small_fixnum_table+SMALL_FIXNUM_LIMIT+(i)) - -struct bignum { - FIRSTWORD; - long *big_self; /* bignum body */ - int big_length; /* bignum length */ -}; -#define MP(x) ((GEN)(x)->big.big_self) -struct shortfloat_struct { - FIRSTWORD; - shortfloat SFVAL; -}; -#define sf(x) (x)->SF.SFVAL -struct longfloat_struct { - FIRSTWORD; - longfloat LFVAL; -}; -#define lf(x) (x)->LF.LFVAL -struct character { - FIRSTWORD; - unsigned short ch_code; - unsigned char ch_font; - unsigned char ch_bits; -}; -struct character character_table1[256+128]; -#define character_table (character_table1+128) -#define code_char(c) (object)(character_table+(c)) -#define char_code(x) (x)->ch.ch_code -#define char_font(x) (x)->ch.ch_font -#define char_bits(x) (x)->ch.ch_bits -enum stype { - stp_ordinary, - stp_constant, - stp_special -}; -struct symbol { - FIRSTWORD; - object s_dbind; - int (*s_sfdef)(); -#define s_fillp st_fillp -#define s_self st_self - int s_fillp; - char *s_self; - object s_gfdef; - object s_plist; - object s_hpack; - short s_stype; - short s_mflag; -}; -struct cons { - FIRSTWORD; - object c_cdr; - object c_car; -}; -struct array { - FIRSTWORD; - short a_rank; - short a_adjustable; - int a_dim; - int *a_dims; - object *a_self; - object a_displaced; - short a_elttype; - short a_offset; -}; - - - -struct fat_string { /* vector header */ - FIRSTWORD; - unsigned fs_raw : 24; /* tells if the things in leader are raw */ - unsigned char fs_leader_length; /* leader_Length */ - int fs_dim; /* dimension */ - int fs_fillp; /* fill pointer */ - /* For simple vectors, */ - /* fs_fillp is equal to fs_dim. */ - fatchar *fs_self; /* pointer to the vector Note the leader starts at (int *) *fs_self - fs_leader_length */ -}; - - -struct vector { - FIRSTWORD; - short v_hasfillp; - short v_adjustable; - int v_dim; - int v_fillp; - object *v_self; - object v_displaced; - short v_elttype; - short v_offset; -}; -struct string { - FIRSTWORD; - short st_hasfillp; - short st_adjustable; - int st_dim; - int st_fillp; - char *st_self; - object st_displaced; -}; -struct ustring { - FIRSTWORD; - short ust_hasfillp; - short ust_adjustable; - int ust_dim; - int ust_fillp; - unsigned char - *ust_self; - object ust_displaced; -}; -#define USHORT(x,i) (((unsigned short *)(x)->ust.ust_self)[i]) - -struct bitvector { - FIRSTWORD; - short bv_hasfillp; - short bv_adjustable; - int bv_dim; - int bv_fillp; - char *bv_self; - object bv_displaced; - short bv_elttype; - short bv_offset; -}; -struct fixarray { - FIRSTWORD; - short fixa_rank; - short fixa_adjustable; - int fixa_dim; - int *fixa_dims; - fixnum *fixa_self; - object fixa_displaced; - short fixa_elttype; - short fixa_offset; -}; -struct sfarray { - FIRSTWORD; - short sfa_rank; - short sfa_adjustable; - int sfa_dim; - int *sfa_dims; - shortfloat - *sfa_self; - object sfa_displaced; - short sfa_elttype; - short sfa_offset; -}; -struct lfarray { - FIRSTWORD; - short lfa_rank; - short lfa_adjustable; - int lfa_dim; - int *lfa_dims; - longfloat - *lfa_self; - object lfa_displaced; - short lfa_elttype; - short lfa_offset; -}; - -struct structure { /* structure header */ - FIRSTWORD; - object str_def; /* structure definition (a structure) */ - object *str_self; /* structure self */ -}; - -#define STREF(type,x,i) (*((type *)(((char *)((x)->str.str_self))+(i)))) - -struct cfun { - FIRSTWORD; - object cf_name; - int (*cf_self)(); - object cf_data; -}; - - struct dclosure { /* compiled closure header */ - FIRSTWORD; - int (*dc_self)(); /* entry address */ - object *dc_env; /* environment */ -}; - - struct cclosure { - FIRSTWORD; - - object cc_name; - int (*cc_self)(); - object cc_env; - object cc_data; - object *cc_turbo; -}; - -struct sfun { - FIRSTWORD; - object sfn_name; - int (*sfn_self)(); - object sfn_data; - int sfn_argd; - - }; -struct vfun { - FIRSTWORD; - object vfn_name; - int (*vfn_self)(); - object vfn_data; - unsigned short vfn_minargs; - unsigned short vfn_maxargs; - }; - -struct dummy { - FIRSTWORD; -}; -struct stream { - FIRSTWORD; - FILE *sm_fp; /* file pointer */ - object sm_object0; /* some object */ - object sm_object1; /* some object */ - int sm_int0; /* some int */ - int sm_int1; /* some int */ - char *sm_buffer; /* ptr to BUFSIZE block of storage */ - short sm_mode; /* stream mode */ - /* of enum smmode */ -}; -union lispunion { - struct fixnum_struct - FIX; - struct shortfloat_struct - SF; - struct stream sm; - struct longfloat_struct - LF; - struct character - ch; - struct symbol s; - struct cons c; - struct array a; - struct vector v; - struct string st; - struct ustring ust; - struct bignum big; - struct bitvector - bv; - struct structure - str; - struct cfun cf; - struct cclosure cc; - struct sfun sfn; - struct vfun vfn; - struct dummy d; - struct fat_string fs; - struct dclosure dc; - struct fixarray fixa; - struct sfarray sfa; - struct lfarray lfa; -}; -enum type { - t_cons, - t_start = 0 , /* t_cons, */ - t_fixnum, - t_bignum, - t_ratio, - t_shortfloat, - t_longfloat, - t_complex, - t_character, - t_symbol, - t_package, - t_hashtable, - t_array, - t_vector, - t_string, - t_bitvector, - t_structure, - t_stream, - t_random, - t_readtable, - t_pathname, - t_cfun, - t_cclosure, - t_sfun, - t_gfun, - t_vfun, - t_cfdata, - t_spice, - t_fat_string, - t_dclosure, - t_end, - t_contiguous, - t_relocatable, - t_other -}; -#define type_of(obje) ((enum type)(((object)(obje))->d.t)) -#define endp(obje) endp1(obje) -extern object value_stack[COM_LENG]; -#define vs_org value_stack -object *vs_limit; -object *vs_base; -object *vs_top; -#define vs_push(obje) (*vs_top++ = (obje)) -#define vs_pop (*--vs_top) -#define vs_head vs_top[-1] -#define vs_mark object *old_vs_top = vs_top -#define vs_reset vs_top = old_vs_top -#define vs_check if (vs_top >= vs_limit) \ - vs_overflow(); -#define vs_check_push(obje) \ - (vs_top >= vs_limit ? \ - (object)vs_overflow() : (*vs_top++ = (obje))) -#define check_arg(n) \ - if (vs_top - vs_base != (n)) \ - check_arg_failed(n) -#define MMcheck_arg(n) \ - if (vs_top - vs_base < (n)) \ - too_few_arguments(); \ - else if (vs_top - vs_base > (n)) \ - too_many_arguments() -#define vs_reserve(x) if(vs_base+(x) >= vs_limit) \ - vs_overflow(); -struct bds_bd { - object bds_sym; - object bds_val; -}; -extern struct bds_bd bind_stack[COM_LENG]; -typedef struct bds_bd *bds_ptr; -bds_ptr bds_org; -bds_ptr bds_limit; -bds_ptr bds_top; -#define bds_check \ - if (bds_top >= bds_limit) \ - bds_overflow() -#define bds_bind(sym, val) \ - ((++bds_top)->bds_sym = (sym), \ - bds_top->bds_val = (sym)->s.s_dbind, \ - (sym)->s.s_dbind = (val)) -#define bds_unwind1 \ - ((bds_top->bds_sym)->s.s_dbind = bds_top->bds_val, --bds_top) -typedef struct invocation_history { - object ihs_function; - object *ihs_base; -} *ihs_ptr; -extern struct invocation_history ihs_stack[COM_LENG]; -ihs_ptr ihs_org; -ihs_ptr ihs_limit; -ihs_ptr ihs_top; -#define ihs_check \ - if (ihs_top >= ihs_limit) \ - ihs_overflow() -#define ihs_push(function) \ - (++ihs_top)->ihs_function = (function); \ - ihs_top->ihs_base = vs_base -#define ihs_pop() (ihs_top--) -enum fr_class { - FRS_CATCH, - FRS_CATCHALL, - FRS_PROTECT -}; -struct frame { - jmp_buf frs_jmpbuf; - object *frs_lex; - bds_ptr frs_bds_top; - enum fr_class frs_class; - object frs_val; - ihs_ptr frs_ihs; -}; -typedef struct frame *frame_ptr; -#define alloc_frame_id() alloc_object(t_spice) -extern struct frame frame_stack[COM_LENG]; - -frame_ptr frs_org; -frame_ptr frs_limit; -frame_ptr frs_top; -#define frs_push(class, val) \ - if (++frs_top >= frs_limit) \ - frs_overflow(); \ - frs_top->frs_lex = lex_env;\ - frs_top->frs_bds_top = bds_top; \ - frs_top->frs_class = (class); \ - frs_top->frs_val = (val); \ - frs_top->frs_ihs = ihs_top; \ - setjmp(frs_top->frs_jmpbuf) -#define frs_pop() frs_top-- -bool nlj_active; -frame_ptr nlj_fr; -object nlj_tag; -object *lex_env; -object caar(); -object cadr(); -object cdar(); -object cddr(); -object caaar(); -object caadr(); -object cadar(); -object caddr(); -object cdaar(); -object cdadr(); -object cddar(); -object cdddr(); -object caaaar(); -object caaadr(); -object caadar(); -object caaddr(); -object cadaar(); -object cadadr(); -object caddar(); -object cadddr(); -object cdaaar(); -object cdaadr(); -object cdadar(); -object cdaddr(); -object cddaar(); -object cddadr(); -object cdddar(); -object cddddr(); -#define MMcons(a,d) make_cons((a),(d)) -#define MMcar(x) (x)->c.c_car -#define MMcdr(x) (x)->c.c_cdr -#define CMPcar(x) (x)->c.c_car -#define CMPcdr(x) (x)->c.c_cdr -#define CMPcaar(x) (x)->c.c_car->c.c_car -#define CMPcadr(x) (x)->c.c_cdr->c.c_car -#define CMPcdar(x) (x)->c.c_car->c.c_cdr -#define CMPcddr(x) (x)->c.c_cdr->c.c_cdr -#define CMPcaaar(x) (x)->c.c_car->c.c_car->c.c_car -#define CMPcaadr(x) (x)->c.c_cdr->c.c_car->c.c_car -#define CMPcadar(x) (x)->c.c_car->c.c_cdr->c.c_car -#define CMPcaddr(x) (x)->c.c_cdr->c.c_cdr->c.c_car -#define CMPcdaar(x) (x)->c.c_car->c.c_car->c.c_cdr -#define CMPcdadr(x) (x)->c.c_cdr->c.c_car->c.c_cdr -#define CMPcddar(x) (x)->c.c_car->c.c_cdr->c.c_cdr -#define CMPcdddr(x) (x)->c.c_cdr->c.c_cdr->c.c_cdr -#define CMPcaaaar(x) (x)->c.c_car->c.c_car->c.c_car->c.c_car -#define CMPcaaadr(x) (x)->c.c_cdr->c.c_car->c.c_car->c.c_car -#define CMPcaadar(x) (x)->c.c_car->c.c_cdr->c.c_car->c.c_car -#define CMPcaaddr(x) (x)->c.c_cdr->c.c_cdr->c.c_car->c.c_car -#define CMPcadaar(x) (x)->c.c_car->c.c_car->c.c_cdr->c.c_car -#define CMPcadadr(x) (x)->c.c_cdr->c.c_car->c.c_cdr->c.c_car -#define CMPcaddar(x) (x)->c.c_car->c.c_cdr->c.c_cdr->c.c_car -#define CMPcadddr(x) (x)->c.c_cdr->c.c_cdr->c.c_cdr->c.c_car -#define CMPcdaaar(x) (x)->c.c_car->c.c_car->c.c_car->c.c_cdr -#define CMPcdaadr(x) (x)->c.c_cdr->c.c_car->c.c_car->c.c_cdr -#define CMPcdadar(x) (x)->c.c_car->c.c_cdr->c.c_car->c.c_cdr -#define CMPcdaddr(x) (x)->c.c_cdr->c.c_cdr->c.c_car->c.c_cdr -#define CMPcddaar(x) (x)->c.c_car->c.c_car->c.c_cdr->c.c_cdr -#define CMPcddadr(x) (x)->c.c_cdr->c.c_car->c.c_cdr->c.c_cdr -#define CMPcdddar(x) (x)->c.c_car->c.c_cdr->c.c_cdr->c.c_cdr -#define CMPcddddr(x) (x)->c.c_cdr->c.c_cdr->c.c_cdr->c.c_cdr -#define CMPfuncall funcall -#define cclosure_call funcall -object simple_lispcall(); -object simple_lispcall_no_event(); -object simple_symlispcall(); -object simple_symlispcall_no_event(); -object CMPtemp; -object CMPtemp1; -object CMPtemp2; -object CMPtemp3; -#define Cnil ((object)&Cnil_body) -#define Ct ((object)&Ct_body) -struct symbol Cnil_body, Ct_body; -object MF(); -object MFnew(); -object MM(); -object Scons; -object siSfunction_documentation; -object siSvariable_documentation; -object siSpretty_print_format; -object Slist; -object keyword_package; -object alloc_object(); -object car(); -object cdr(); -object list(); -object listA(); -object coerce_to_string(); -object elt(); -object elt_set(); -frame_ptr frs_sch(); -frame_ptr frs_sch_catch(); -object make_cclosure(); -object make_cclosure_new(); -object nth(); -object nthcdr(); -object make_cons(); -object append(); -object nconc(); -object reverse(); -object nreverse(); -object number_expt(); -object number_minus(); -object number_negate(); -object number_plus(); -object number_times(); -object one_minus(); -object one_plus(); -object get(); -object getf(); -object putprop(); -object sputprop(); -object remprop(); -object string_to_object(); -object symbol_function(); -object symbol_value(); -object make_fixnum(); -object make_shortfloat(); -object make_longfloat(); -object structure_ref(); -object structure_set(); -object princ(); -object prin1(); -object print(); -object terpri(); -object aref(); -object aset(); -object aref1(); -object aset1(); -void call_or_link(); -object call_proc(); -object call_proc0(); -object call_proc1(); -object call_proc2(); -object ifuncall(); -object ifuncall1(); -object ifuncall2(); -object symbol_name(); -char object_to_char(); -int object_to_int(); -float object_to_float(); -double object_to_double(); -char *object_to_string(); -int FIXtemp; -#define CMPmake_fixnum(x) \ -((((FIXtemp=(x))+1024)&-2048)==0?small_fixnum(FIXtemp):make_fixnum(FIXtemp)) -#define Creturn(v) return((vs_top=vs,(v))) -#define Cexit return((vs_top=vs,0)) -double sin(), cos(), tan(); -object read_byte1(),read_char1(); - -#define fs_leader(ar,i) (((object *)((ar)->fs.fs_self))[-(i+1)]) -#define RPAREN ) -object make_list(); -#ifdef HAVE_ALLOCA -#ifndef alloca -char *alloca(); -#endif -char *alloca_val; -#define ALLOCA_CONS(n) (alloca_val=alloca((n)*sizeof(struct cons))) -#define ON_STACK_CONS(x,y) (alloca_val=alloca(sizeof(struct cons)), on_stack_cons(x,y)) -#define ON_STACK_LIST on_stack_list -#define ON_STACK_LIST_VECTOR on_stack_list_vector -#define ON_STACK_MAKE_LIST on_stack_make_list -object on_stack_cons(); -object on_stack_list(); -object on_stack_list_vector(); -object on_stack_make_list(); -#else -#define ALLOCA_CONS(n) 0 -#define ON_STACK_CONS(x,y) MMcons(x,y) -#define ON_STACK_LIST list -#define ON_STACK_LIST_VECTOR list_vector -#define ON_STACK_MAKE_LIST make_list -#endif - - -struct call_data { object fun; - int argd;}; -struct call_data fcall; -object fcalln(); -object list_vector(); -object MVloc[10]; -#define VARG(min,max) ((min) | (max << 8)) -#define VFUN_NARGS fcall.argd -extern object Cstd_key_defaults[]; -int vfun_wrong_number_of_args(); -int eql(),equal(),eq(); -object sublis1(); -object LVformat(),LVerror(); -#define EQ(x,y) ((x)==(y)) - - - -/* #include "../h/genpari.h"*/ -typedef unsigned long *GEN; -GEN addii(),mulii(),mulsi(),powerii(),shifti(),stoi(),dvmdii(),subii(); -int cmpii(); -#define signe(x) (((GEN)(x))[1]>>24) -#define lg(x) (((GEN)(x))[0]&0xffff) -#define setlg(x,s) (((GEN)(x))[0]=(((GEN)(x))[0]&0xffff0000)+s) -#define lgef(x) (((GEN)(x))[1]&0xffff) -#define setlgef(x,s) (((GEN)(x))[1]=(((GEN)(x))[1]&0xffff0000)+s) - -int in_saved_avma ; -#define ulong unsigned long -/* #define DEBUG_AVMA */ - -#ifdef DEBUG_AVMA -#define save_avma long lvma = (in_saved_avma = 1, avma) -#define restore_avma avma = (in_saved_avma = 0, lvma) -#else -#define save_avma long lvma = avma -#define restore_avma avma = lvma -#endif -unsigned long avma; -GEN gzero; -GEN icopy_x; - -object make_integer(); - /* copy x to y, increasing space by factor of 2 */ - - -GEN otoi(); -/* -object integ_temp; -#define otoi(x) (integ_temp = (x) , (type_of(integ_temp) == t_bignum \ - ? MP(integ_temp) :stoi(fix(integ_temp)))) -*/ - -void isetq_fix(); -#ifdef HAVE_ALLOCA -#define SETQ_II(var,alloc,val) \ - do{GEN _xx =(val) ; \ - int _n = replace_copy1(_xx,var); \ - if(_n) var = replace_copy2(_xx,alloca(_n));}while(0) - -#define SETQ_IO(var,alloc,val) {object _xx =(val) ; \ - int _n = obj_replace_copy1(_xx,var); \ - if(_n) var = obj_replace_copy2(_xx,alloca(_n));} -#define IDECL(a,b,c) ulong b[4];a =(b[0]=0x1010000 +4,b) -#else -GEN setq_io(),setq_ii(); -#define SETQ_IO(x,alloc,val) (x)=setq_io(x,&alloc,val) -#define SETQ_II(x,alloc,val) (x)=setq_ii(x,&alloc,val) -#define IDECL(a,b,c) ulong b[4];a =(b[0]=0x1010000 +4,b);object c -#endif - - -#ifdef __GNUC__ -#define alloca __builtin_alloca -#endif - - diff --git a/lib/complex.c b/lib/complex.c deleted file mode 100644 index 4056b2e..0000000 --- a/lib/complex.c +++ /dev/null @@ -1,131 +0,0 @@ -/* complex - Complex number functions */ -/* Copyright (c) 1990, by Luke Tierney */ - -/* patched up and semi-ansified, (c) 2006, AJ Rossini, blindglobe@gmail.com */ - -#include "xmath.h" -#include "complex.h" - -extern void xlfail(char *); - -/*static */ -double -phase(Complex c) -{ - double phi; - - if (c.real == 0.0) { - if (c.imag > 0.0) phi = PI / 2; - else if (c.imag == 0.0) phi = 0.0; - else phi = -PI / 2; - } else { - phi = atan(c.imag / c.real); - if (c.real < 0.0) { - if (c.imag > 0.0) { - phi += PI; - } else { - if (c.imag < 0.0) { - phi -= PI; - } else { - phi = PI; - } - } - } - } - return(phi); -} - -double -modulus(Complex c) -{ - return(sqrt(c.real * c.real + c.imag * c.imag)); -} - -Complex -cart2complex(double real, double imag) -{ - Complex val; - val.real = real; - val.imag = imag; - return(val); -} - -Complex -polar2complex(double mod, double phi) -{ - Complex val; - double cs, sn; - - if (phi == 0) { - cs = 1.0; - sn = 0.0; - } else { - if (phi == PI / 2) { - cs = 0.0; - sn = 1.0; - } else { - if (phi == PI) { - cs = -1.0; - sn = 0.0; - } else { - if (phi == -PI / 2) { - cs = 0.0; - sn = -1.0; - } else { - cs = cos(phi); - sn = sin(phi); - } - } - } - } - val.real = mod * cs; - val.imag = mod * sn; - return(val); -} - -Complex -cadd(Complex c1, - Complex c2) -{ - return(cart2complex(c1.real + c2.real, c1.imag + c2.imag)); -} - -Complex -csub(Complex c1, - Complex c2) -{ - return(cart2complex(c1.real - c2.real, c1.imag - c2.imag)); -} - -Complex -cmul(Complex c1, Complex c2) -{ - double m1, m2, p1, p2; - - m1 = modulus(c1); - p1 = phase(c1); - m2 = modulus(c2); - p2 = phase(c2); - return(polar2complex(m1 * m2, p1 + p2)); -} - -static void -checkfzero(double x) -{ - if (x == 0.0) { - xlfail("division by zero"); - } -} - -Complex -cdiv(Complex c1, Complex c2) -{ - double m1, m2, p1, p2; - - m1 = modulus(c1); - p1 = phase(c1); - m2 = modulus(c2); - p2 = phase(c2); - checkfzero(m2); - return(polar2complex(m1 / m2, p1 - p2)); -} diff --git a/lib/complex.h b/lib/complex.h deleted file mode 100644 index 3658a8f..0000000 --- a/lib/complex.h +++ /dev/null @@ -1,18 +0,0 @@ -/* I hate C programming... */ - -typedef struct { - double real, imag; -} Complex; - -#ifndef PI -#define PI 3.141592653589793 -#endif /* PI */ - -extern Complex makecomplex(); -extern Complex cart2complex(); -extern Complex polar2complex(double, double); -/* extern Complex csqrt(), cexp(), clog(), cexpt(), csin(), ccos(), ctan(); - extern Complex casin(), cacos(), catan(); */ -extern Complex cadd(), csub(), cmul(), cdiv(); -extern /* static */ double phase(Complex); -extern double modulus(Complex); diff --git a/lib/derivatives.c b/lib/derivatives.c deleted file mode 100644 index f7ff8c4..0000000 --- a/lib/derivatives.c +++ /dev/null @@ -1,71 +0,0 @@ -/* derivatives - for Bayes code in XLISP-STAT and S */ -/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */ -/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */ -/* You may give out copies of this software; for conditions see the */ -/* file COPYING included with this distribution. */ - -/* -#ifdef __APPLE__ -#include -#endif -*/ - -#include - -# include "xmath.h" -extern double macheps(); - -# define nil 0L -# define FALSE 0 -# define TRUE 1 - -void -numergrad(size_t n, double *x, double *grad, double *fsum, - int ffun(double *, double *, double *, double **), - double h, double *typx) - /* int (*ffun)();*/ -{ - size_t i; - double old_xi, f1, f2, hi; - - for (i = 0; i < n; i++) { - old_xi = x[i]; - hi = (typx != nil) ? typx[i] * h : h; - x[i] = old_xi + hi; - (*ffun)(x, &f1, nil, nil); - x[i] = old_xi - hi; - (*ffun)(x, &f2, nil, nil); - x[i] = old_xi; - grad[i] = (f1 - f2) / (2.0 * hi); - fsum[i] = f1 + f2; - } -} - -void -numerhess(size_t n, double *x, double **hess, double f, double *fsum, - int *ffun(double *, double *, double *, double *), - double h, double *typx) -{ - size_t i, j; - double old_xi, old_xj, f1, f2, hi, hj; - - for (i = 0; i < n; i++) { - hi = (typx != nil) ? typx[i] * h : h; - hess[i][i] = (fsum[i] - 2 * f) / (hi * hi); - for (j = i + 1; j < n; j++) { - hj = (typx != nil) ? typx[j] * h : h; - old_xi = x[i]; - old_xj = x[j]; - x[i] = old_xi + hi; - x[j] = old_xj + hj; - (*ffun)(x, &f1, nil, nil); - x[i] = old_xi - hi; - x[j] = old_xj - hj; - (*ffun)(x, &f2, nil, nil); - x[i] = old_xi; - x[j] = old_xj; - hess[i][j] = (2 * f + f1 + f2 - fsum[i] - fsum[j]) / (2.0 * hi * hj); - hess[j][i] = hess[i][j]; - } - } -} diff --git a/lib/eigen.c b/lib/eigen.c deleted file mode 100644 index e45d4f1..0000000 --- a/lib/eigen.c +++ /dev/null @@ -1,625 +0,0 @@ -/* eigen.f -- translated by f2c (version of 19 December 1990 16:50:21). - and modified. */ - -#include "xmath.h" - -#define integer int -#define real float -#define doublereal double -#define Min(x,y) ((x) > (y) ? (y) : (x)) -#define Max(x,y) ((x) > (y) ? (x) : (y)) -#define Abs(x) ((x) >= 0 ? (x) : -(x)) - -static /* Subroutine */ int tred2(), tql2(); -static doublereal pythag(), d_sign(); - -/* Table of constant values */ - -static doublereal c_b39 = 1.; - -/* Subroutine */ int eigen(nm, n, a, w, z, fv1, ierr) -integer *nm, *n; -doublereal *a, *w, *z, *fv1; -integer *ierr; -{ - /* System generated locals */ - integer a_dim1, a_offset, z_dim1, z_offset; - -/* THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF */ -/* SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) */ -/* TO FIND THE EIGENVALUES AND EIGENVECTORS */ -/* OF A REAL SYMMETRIC MATRIX. */ - -/* ON INPUT */ - -/* NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL */ -/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ -/* DIMENSION STATEMENT. */ - -/* N IS THE ORDER OF THE MATRIX A. */ - -/* A CONTAINS THE REAL SYMMETRIC MATRIX. */ - -/* ON OUTPUT */ - -/* W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. */ - -/* Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. */ - -/* IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR */ -/* COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT */ -/* AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. */ - -/* FV1 IS A TEMPORARY STORAGE ARRAY. */ - -/* ------------------------------------------------------------------ -*/ - - /* Parameter adjustments */ - --fv1; - z_dim1 = *nm; - z_offset = z_dim1 + 1; - z -= z_offset; - --w; - a_dim1 = *nm; - a_offset = a_dim1 + 1; - a -= a_offset; - - /* Function Body */ - if (*n <= *nm) { - goto L10; - } - *ierr = *n * 10; - goto L50; -/* .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... */ -L10: - tred2(nm, n, &a[a_offset], &w[1], &fv1[1], &z[z_offset]); - tql2(nm, n, &w[1], &fv1[1], &z[z_offset], ierr); -L50: - return 0; -} /* eigen */ - -static /* Subroutine */ int tred2(nm, n, a, d, e, z) -integer *nm, *n; -doublereal *a, *d, *e, *z; -{ - /* System generated locals */ - integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3; - doublereal d__1; - - /* Local variables */ - static doublereal f, g, h; - static integer i, j, k, l; - static doublereal scale, hh; - static integer ii, jp1; - - - -/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2, */ -/* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. */ -/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). */ - -/* THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A */ -/* SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING */ -/* ORTHOGONAL SIMILARITY TRANSFORMATIONS. */ - -/* ON INPUT */ - -/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ -/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ -/* DIMENSION STATEMENT. */ - -/* N IS THE ORDER OF THE MATRIX. */ - -/* A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE */ -/* LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. */ - -/* ON OUTPUT */ - -/* D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */ - -/* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */ -/* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */ - -/* Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX */ -/* PRODUCED IN THE REDUCTION. */ - -/* A AND Z MAY COINCIDE. IF DISTINCT, A IS UNALTERED. */ - -/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ -/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY -*/ - -/* THIS VERSION DATED AUGUST 1983. */ - -/* ------------------------------------------------------------------ -*/ - - /* Parameter adjustments */ - z_dim1 = *nm; - z_offset = z_dim1 + 1; - z -= z_offset; - --e; - --d; - a_dim1 = *nm; - a_offset = a_dim1 + 1; - a -= a_offset; - - /* Function Body */ - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - - i__2 = *n; - for (j = i; j <= i__2; ++j) { -/* L80: */ - z[j + i * z_dim1] = a[j + i * a_dim1]; - } - - d[i] = a[*n + i * a_dim1]; -/* L100: */ - } - - if (*n == 1) { - goto L510; - } -/* .......... FOR I=N STEP -1 UNTIL 2 DO -- .......... */ - i__1 = *n; - for (ii = 2; ii <= i__1; ++ii) { - i = *n + 2 - ii; - l = i - 1; - h = 0.; - scale = 0.; - if (l < 2) { - goto L130; - } -/* .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... */ - i__2 = l; - for (k = 1; k <= i__2; ++k) { -/* L120: */ - scale += (d__1 = d[k], Abs(d__1)); - } - - if (scale != 0.) { - goto L140; - } -L130: - e[i] = d[l]; - - i__2 = l; - for (j = 1; j <= i__2; ++j) { - d[j] = z[l + j * z_dim1]; - z[i + j * z_dim1] = 0.; - z[j + i * z_dim1] = 0.; -/* L135: */ - } - - goto L290; - -L140: - i__2 = l; - for (k = 1; k <= i__2; ++k) { - d[k] /= scale; - h += d[k] * d[k]; -/* L150: */ - } - - f = d[l]; - d__1 = sqrt(h); - g = -d_sign(&d__1, &f); - e[i] = scale * g; - h -= f * g; - d[l] = f - g; -/* .......... FORM A*U .......... */ - i__2 = l; - for (j = 1; j <= i__2; ++j) { -/* L170: */ - e[j] = 0.; - } - - i__2 = l; - for (j = 1; j <= i__2; ++j) { - f = d[j]; - z[j + i * z_dim1] = f; - g = e[j] + z[j + j * z_dim1] * f; - jp1 = j + 1; - if (l < jp1) { - goto L220; - } - - i__3 = l; - for (k = jp1; k <= i__3; ++k) { - g += z[k + j * z_dim1] * d[k]; - e[k] += z[k + j * z_dim1] * f; -/* L200: */ - } - -L220: - e[j] = g; -/* L240: */ - } -/* .......... FORM P .......... */ - f = 0.; - - i__2 = l; - for (j = 1; j <= i__2; ++j) { - e[j] /= h; - f += e[j] * d[j]; -/* L245: */ - } - - hh = f / (h + h); -/* .......... FORM Q .......... */ - i__2 = l; - for (j = 1; j <= i__2; ++j) { -/* L250: */ - e[j] -= hh * d[j]; - } -/* .......... FORM REDUCED A .......... */ - i__2 = l; - for (j = 1; j <= i__2; ++j) { - f = d[j]; - g = e[j]; - - i__3 = l; - for (k = j; k <= i__3; ++k) { -/* L260: */ - z[k + j * z_dim1] = z[k + j * z_dim1] - f * e[k] - g * d[k]; - } - - d[j] = z[l + j * z_dim1]; - z[i + j * z_dim1] = 0.; -/* L280: */ - } - -L290: - d[i] = h; -/* L300: */ - } -/* .......... ACCUMULATION OF TRANSFORMATION MATRICES .......... */ - i__1 = *n; - for (i = 2; i <= i__1; ++i) { - l = i - 1; - z[*n + l * z_dim1] = z[l + l * z_dim1]; - z[l + l * z_dim1] = 1.; - h = d[i]; - if (h == 0.) { - goto L380; - } - - i__2 = l; - for (k = 1; k <= i__2; ++k) { -/* L330: */ - d[k] = z[k + i * z_dim1] / h; - } - - i__2 = l; - for (j = 1; j <= i__2; ++j) { - g = 0.; - - i__3 = l; - for (k = 1; k <= i__3; ++k) { -/* L340: */ - g += z[k + i * z_dim1] * z[k + j * z_dim1]; - } - - i__3 = l; - for (k = 1; k <= i__3; ++k) { - z[k + j * z_dim1] -= g * d[k]; -/* L360: */ - } - } - -L380: - i__3 = l; - for (k = 1; k <= i__3; ++k) { -/* L400: */ - z[k + i * z_dim1] = 0.; - } - -/* L500: */ - } - -L510: - i__1 = *n; - for (i = 1; i <= i__1; ++i) { - d[i] = z[*n + i * z_dim1]; - z[*n + i * z_dim1] = 0.; -/* L520: */ - } - - z[*n + *n * z_dim1] = 1.; - e[1] = 0.; - return 0; -} /* tred2 */ - -static /* Subroutine */ int tql2(nm, n, d, e, z, ierr) -integer *nm, *n; -doublereal *d, *e, *z; -integer *ierr; -{ - /* System generated locals */ - integer z_dim1, z_offset, i__1, i__2, i__3; - doublereal d__1, d__2; - - /* Local variables */ - static doublereal c, f, g, h; - static integer i, j, k, l, m; - static doublereal p, r, s, c2, c3; - static integer l1, l2; - static doublereal s2; - static integer ii; - static doublereal dl1, el1; - static integer mml; - static doublereal tst1, tst2; - - - -/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2, */ -/* NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND */ -/* WILKINSON. */ -/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). */ - -/* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */ -/* OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD. */ -/* THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO */ -/* BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS */ -/* FULL MATRIX TO TRIDIAGONAL FORM. */ - -/* ON INPUT */ - -/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ -/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ -/* DIMENSION STATEMENT. */ - -/* N IS THE ORDER OF THE MATRIX. */ - -/* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */ - -/* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */ -/* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */ - -/* Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE */ -/* REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS */ -/* OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN */ -/* THE IDENTITY MATRIX. */ - -/* ON OUTPUT */ - -/* D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN */ -/* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT */ -/* UNORDERED FOR INDICES 1,2,...,IERR-1. */ - -/* E HAS BEEN DESTROYED. */ - -/* Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC */ -/* TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, */ -/* Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED */ -/* EIGENVALUES. */ - -/* IERR IS SET TO */ -/* ZERO FOR NORMAL RETURN, */ -/* J IF THE J-TH EIGENVALUE HAS NOT BEEN */ -/* DETERMINED AFTER 30 ITERATIONS. */ - -/* CALLS PYTHAG FOR DSQRT(A*A + B*B) . */ - -/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ -/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY -*/ - -/* THIS VERSION DATED AUGUST 1983. */ - -/* ------------------------------------------------------------------ -*/ - - /* Parameter adjustments */ - z_dim1 = *nm; - z_offset = z_dim1 + 1; - z -= z_offset; - --e; - --d; - - /* Function Body */ - *ierr = 0; - if (*n == 1) { - goto L1001; - } - - i__1 = *n; - for (i = 2; i <= i__1; ++i) { -/* L100: */ - e[i - 1] = e[i]; - } - - f = 0.; - tst1 = 0.; - e[*n] = 0.; - - i__1 = *n; - for (l = 1; l <= i__1; ++l) { - j = 0; - h = (d__1 = d[l], Abs(d__1)) + (d__2 = e[l], Abs(d__2)); - if (tst1 < h) { - tst1 = h; - } -/* .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... */ - i__2 = *n; - for (m = l; m <= i__2; ++m) { - tst2 = tst1 + (d__1 = e[m], Abs(d__1)); - if (tst2 == tst1) { - goto L120; - } -/* .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT */ -/* THROUGH THE BOTTOM OF THE LOOP .......... */ -/* L110: */ - } - -L120: - if (m == l) { - goto L220; - } -L130: - if (j == 30) { - goto L1000; - } - ++j; -/* .......... FORM SHIFT .......... */ - l1 = l + 1; - l2 = l1 + 1; - g = d[l]; - p = (d[l1] - g) / (e[l] * 2.); - r = pythag(&p, &c_b39); - d[l] = e[l] / (p + d_sign(&r, &p)); - d[l1] = e[l] * (p + d_sign(&r, &p)); - dl1 = d[l1]; - h = g - d[l]; - if (l2 > *n) { - goto L145; - } - - i__2 = *n; - for (i = l2; i <= i__2; ++i) { -/* L140: */ - d[i] -= h; - } - -L145: - f += h; -/* .......... QL TRANSFORMATION .......... */ - p = d[m]; - c = 1.; - c2 = c; - el1 = e[l1]; - s = 0.; - mml = m - l; -/* .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... */ - i__2 = mml; - for (ii = 1; ii <= i__2; ++ii) { - c3 = c2; - c2 = c; - s2 = s; - i = m - ii; - g = c * e[i]; - h = c * p; - r = pythag(&p, &e[i]); - e[i + 1] = s * r; - s = e[i] / r; - c = p / r; - p = c * d[i] - s * g; - d[i + 1] = h + s * (c * g + s * d[i]); -/* .......... FORM VECTOR .......... */ - i__3 = *n; - for (k = 1; k <= i__3; ++k) { - h = z[k + (i + 1) * z_dim1]; - z[k + (i + 1) * z_dim1] = s * z[k + i * z_dim1] + c * h; - z[k + i * z_dim1] = c * z[k + i * z_dim1] - s * h; -/* L180: */ - } - -/* L200: */ - } - - p = -s * s2 * c3 * el1 * e[l] / dl1; - e[l] = s * p; - d[l] = c * p; - tst2 = tst1 + (d__1 = e[l], Abs(d__1)); - if (tst2 > tst1) { - goto L130; - } -L220: - d[l] += f; -/* L240: */ - } -/* .......... ORDER EIGENVALUES AND EIGENVECTORS .......... */ - i__1 = *n; - for (ii = 2; ii <= i__1; ++ii) { - i = ii - 1; - k = i; - p = d[i]; - - i__2 = *n; - for (j = ii; j <= i__2; ++j) { - if (d[j] >= p) { - goto L260; - } - k = j; - p = d[j]; -L260: - ; - } - - if (k == i) { - goto L300; - } - d[k] = d[i]; - d[i] = p; - - i__2 = *n; - for (j = 1; j <= i__2; ++j) { - p = z[j + i * z_dim1]; - z[j + i * z_dim1] = z[j + k * z_dim1]; - z[j + k * z_dim1] = p; -/* L280: */ - } - -L300: - ; - } - - goto L1001; -/* .......... SET ERROR -- NO CONVERGENCE TO AN */ -/* EIGENVALUE AFTER 30 ITERATIONS .......... */ -L1000: - *ierr = l; -L1001: - return 0; -} /* tql2 */ - -static doublereal pythag(a, b) -doublereal *a, *b; -{ - /* System generated locals */ - doublereal ret_val, d__1, d__2, d__3; - - /* Local variables */ - static doublereal p, r, s, t, u; - - -/* FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW */ - - -/* Computing MAX */ - d__1 = Abs(*a), d__2 = Abs(*b); - p = Max(d__1,d__2); - if (p == 0.) { - goto L20; - } -/* Computing MIN */ - d__2 = Abs(*a), d__3 = Abs(*b); -/* Computing 2nd power */ - d__1 = Min(d__2,d__3) / p; - r = d__1 * d__1; -L10: - t = r + 4.; - if (t == 4.) { - goto L20; - } - s = r / t; - u = s * 2. + 1.; - p = u * p; -/* Computing 2nd power */ - d__1 = s / u; - r = d__1 * d__1 * r; - goto L10; -L20: - ret_val = p; - return ret_val; -} /* pythag */ - -static double d_sign(a,b) - doublereal *a, *b; -{ - double x; - x = (*a >= 0 ? *a : - *a); - return( *b >= 0 ? x : -x); -} diff --git a/lib/functions.c b/lib/functions.c deleted file mode 100644 index 6d23f8d..0000000 --- a/lib/functions.c +++ /dev/null @@ -1,1115 +0,0 @@ -/* functions - callbacks for Bayes routines in XLISP-STAT and S */ -/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */ -/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */ -/* You may give out copies of this software; for conditions see the */ -/* file COPYING included with this distribution. */ - -#include -#include "xmath.h" - -#define PRINTSTR(s) printf(s) - -extern void *S_alloc(), *calloc(), *realloc(); -extern double macheps(); -char *minresultstring(); - -extern void choldecomp(); - -extern void minsupplyvalues(double , double *, double **, - double *, double **); -extern void minimize(); -extern void minresults(double *, double *, double *, double *, double **, - double *, double **, int *, int *, double *); -extern void minsetup(size_t n, size_t k, - int *ffun(), int *gfun(), - double *x, double typf, double *typx, char *work); - /* int (*ffun)(), (*gfun)();*/ -extern void minsetoptions(double , double, double, int , int , int , int); - -extern void choldecomp(); - - - -/* next 2 from cbayes.c */ -extern void Recover(char *, char *); -extern void call_S(char *fun, long narg, char **args, char **mode, long *length,char **names, - long nvals, char **values); - -/* next 2 from derivatives.c */ -/* -extern void numergrad(size_t n, double *x, double *grad, double *fsum, - int *ffun(), double h, double *typx); -extern void numerhess(size_t n, double *x, double **hess, double f, - double *fsum, void ffun(), - double h, double *typx); -*/ -extern void numergrad(size_t , double *, double *, double *, - int ffun(double *, double *, double *, double **), - double , double *); -extern void numerhess(size_t , double *, double **, double , - double *, - int ffun(double *, double *, double *, double **), - double , double *); - -/* next from minimize */ -extern size_t minworkspacesize(size_t, size_t); - -/************************************************************************/ -/** **/ -/** Definitions and Globals **/ -/** **/ -/************************************************************************/ - -/* #define NULL 0L already defined in stddef */ - -#define nil 0L -#define TRUE 1 -#define FALSE 0 - -#define ROOT2PI 2.50662827463100050241 -#define PI_INV .31830988618379067153 - -#define GRADTOL_POWER 1.0 / 3.0 -#define H_POWER 1.0 / 6.0 - -/*typedef double **RMatrix, *RVector;*/ - -typedef struct{ - char *f, **sf, **g; - size_t n, k; - int change_sign, fderivs; - int *gderivs; - double typf, h, dflt; - double *typx, *fsum, *cvals, *ctarget; - double **gfsum; -} Fundata; - -static Fundata func, gfuncs, cfuncs; - -/************************************************************************/ -/** **/ -/** Memory Utilities **/ -/** **/ -/************************************************************************/ - -/* this function is used to maintain a statically allocated piece of */ -/* memory of a specified size. If a larger piece is needed the pointer */ -/* is realloced. This allows functions using memory allocation to be */ -/* called repeatedly (but not recursively) from within the same call */ -/* from S. It attempts to avoid the danger of dangling callocs. */ - -void static -makespace(char **pptr, size_t size) -{ - if (size <= 0) return; - if (*pptr == nil) *pptr = calloc(size, 1); - else *pptr = realloc(*pptr, size); - if (size > 0 && *pptr == nil) Recover("memory allocation failed", NULL); -} - -/************************************************************************/ -/** **/ -/** Functions Evaluation Routines **/ -/** **/ -/************************************************************************/ - -/* - * All Hessian evaluations by numerical derivatives assume the gradient is - * evaluated first at the same location. The results are cached away. - */ - -/* install log posterior function */ -static void -install_func(char *f, char **sf, size_t n, int change_sign, /*?? */ - double typf, double h, double *typx, double dflt) -{ - int i; - static int inited = FALSE; - - if (! inited) { - func.typx = nil; - func.fsum = nil; - inited = TRUE; - } - makespace(&func.typx, n * sizeof(double)); - makespace(&func.fsum, n * sizeof(double)); - - func.f = f; - func.sf = sf; - func.n = n; - func.change_sign = change_sign; - func.typf = (typf > 0.0) ? typf : 1.0; - func.h = (h > 0.0) ? h : pow(macheps(), H_POWER); - for (i = 0; i < n; i++) - func.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0; - func.dflt = dflt; - func.fderivs = 0; -} - -/* install tilt functions */ -static void -install_gfuncs(char **g, size_t n, size_t k, int change_sign, double h, double *typx) -{ - size_t i; - static int inited = FALSE; - static double *gfsumdata = nil; - - if (! inited) { - gfuncs.typx = nil; - gfuncs.gfsum = nil; - gfuncs.gderivs = nil; - inited = TRUE; - } - makespace(&gfuncs.typx, n * sizeof(double)); - makespace(&gfuncs.gfsum, k * sizeof(double *)); - makespace(&gfsumdata, k * n * sizeof(double)); - makespace(&gfuncs.gderivs, k *sizeof(int)); - - gfuncs.g = g; - gfuncs.n = n; - gfuncs.k = k; - gfuncs.change_sign = change_sign; - gfuncs.h = (h > 0.0) ? h : pow(macheps(), H_POWER); - for (i = 0; i < n; i++) - gfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0; - for (i = 0; i < k; i++) gfuncs.gfsum[i] = gfsumdata + i * n; - return; -} - -/* install constraint functions */ -static void -install_cfuncs(char **g, size_t n, size_t k, double *ctarget, double h, double *typx) -{ - size_t i; - static int inited = FALSE; - - if (! inited) { - cfuncs.typx = nil; - cfuncs.fsum = nil; - cfuncs.gderivs = nil; - inited = TRUE; - } - makespace(&cfuncs.typx, n * sizeof(double)); - makespace(&cfuncs.fsum, n * sizeof(double)); - makespace(&cfuncs.gderivs, k * sizeof(int)); - - cfuncs.g = g; - cfuncs.n = n; - cfuncs.k = k; - cfuncs.h = (h > 0.0) ? h : pow(macheps(), H_POWER); - for (i = 0; i < n; i++) - cfuncs.typx[i] = (typx != nil && typx[i] > 0.0) ? typx[i] : 1.0; - cfuncs.ctarget = ctarget; - return; -} - -static int -in_support(char **ff, size_t n, double *x) -{ - char *args[1], *values[1]; - int *result; - char *mode[1]; - long length[1]; - - if (ff == nil || ff[0] == nil) { - return(TRUE); - } else { - mode[0] = "double"; - length[0] =n; - args[0] = (char *) x; - call_S(ff[0], 1L, args, mode, length, 0L, 1L, values); - result = (int *) values[0]; - return(result[0]); - } -} - -/* callback for logposterior evaluation */ -static int -evalfunc(double *x, double *pval, double *grad, double **hess) -{ - char *args[1], *values[3]; - double *result, val; - char *mode[1]; - long length[1]; - int i, j; - - for (i = 0; i < 3; i++) values[i] = nil; - - if (in_support(func.sf, func.n, x)) { - if (pval != nil || func.fderivs > 0 || hess != nil) { - mode[0] = "double"; - length[0] = func.n; - args[0] = (char *) x; - call_S(func.f, 1L, args, mode, length, 0L, 3L, values); - result = (double *) values[0]; - val = (! func.change_sign) ? result[0] : -result[0]; - if (pval != nil) *pval = val; - if (values[2] != nil) func.fderivs = 2; - else if (values[1] != nil) func.fderivs = 1; - else func.fderivs = 0; - } - if (grad != nil) { - if (func.fderivs > 0) { - result = (double *) values[1]; - for (i = 0; i < func.n; i++) - grad[i] = (! func.change_sign) ? result[i] : -result[i]; - } - else { - numergrad(func.n, x, grad, func.fsum, evalfunc, func.h, func.typx); - } - } - if (hess != nil) { - if (func.fderivs == 2) { - result = (double *) values[2]; - for (i = 0; i < func.n; i++) - for (j = 0; j < func.n; j++) - hess[i][j] = (! func.change_sign) ? result[i + j * func.n] - : -result[i + j * func.n]; - } - else { - if (func.fderivs == 1) /* kludge to get fsum for analytic gradients */ - numergrad(func.n, x, func.fsum, func.fsum, - evalfunc, func.h, func.typx); - numerhess(func.n, x, hess, val, func.fsum, evalfunc, func.h, func.typx); - } - } - return(TRUE); - } else { - if (pval != nil) { - *pval = func.dflt; - } - return(FALSE); - } - return(TRUE); -} - - -/* callback for tilt function evaluation */ -static size_t which_gfunc; - -static int -evalgfunc(double *x, double *pval, double *grad, double **hess) -{ - char *args[1], *values[3]; - double *result, val; - char *mode[1]; - long length[1]; - int i, j; - - for (i = 0; i < 3; i++) values[i] = nil; - - if (pval != nil || gfuncs.gderivs[which_gfunc] > 0 || hess != nil) { - mode[0] = "double"; - length[0] = gfuncs.n; - args[0] = (char *) x; - call_S(gfuncs.g[which_gfunc], 1L, args, mode, length, 0L, 3L, values); - result = (double *) values[0]; - val = result[0]; - if (pval != nil) *pval = result[0]; - if (values[2] != nil) gfuncs.gderivs[which_gfunc] = 2; - else if (values[1] != nil) gfuncs.gderivs[which_gfunc] = 1; - else gfuncs.gderivs[which_gfunc] = 0; - } - if (grad != nil) { - if (gfuncs.gderivs[which_gfunc] > 0) { - result = (double *) values[1]; - for (i = 0; i < gfuncs.n; i++) grad[i] = result[i]; - } - else { - numergrad(gfuncs.n, x, grad, gfuncs.gfsum[which_gfunc], evalgfunc, - gfuncs.h, gfuncs.typx); - } - } - if (hess != nil) { - if (gfuncs.gderivs[which_gfunc] == 2) { - result = (double *) values[2]; - for (i = 0; i < gfuncs.n; i++) - for (j = 0; j < gfuncs.n; j++) - hess[i][j] = result[i + j * gfuncs.n]; - } - else { - /* kludge to get fsum if analytic gradient used */ - if (gfuncs.gderivs[which_gfunc] == 1) - numergrad(gfuncs.n, x, gfuncs.gfsum[which_gfunc], - gfuncs.gfsum[which_gfunc], evalgfunc, gfuncs.h, gfuncs.typx); - numerhess(gfuncs.n, x, hess, val, gfuncs.gfsum[which_gfunc], evalgfunc, - gfuncs.h, gfuncs.typx); - } - } - return(TRUE); -} - -/* callback for constraint function evaluation */ -static int which_cfunc; - -static int -evalcfunc(double *x, double *pval, double *grad, double **hess) -{ - char *args[1], *values[3]; - double *result, val; - char *mode[1]; - long length[1]; - int i, j; - - if (pval != nil || cfuncs.gderivs[which_cfunc] > 0 || hess != nil) { - mode[0] = "double"; - length[0] = cfuncs.n; - args[0] = (char *) x; - call_S(cfuncs.g[which_cfunc], 1L, args, mode, length, 0L, 3L, values); - result = (double *) values[0]; - val = result[0]; - if (pval != nil) { - *pval = result[0]; - if (cfuncs.ctarget != nil) *pval -= cfuncs.ctarget[which_cfunc]; - } - if (values[2] != nil) cfuncs.gderivs[which_cfunc] = 2; - else if (values[1] != nil) cfuncs.gderivs[which_cfunc] = 1; - else cfuncs.gderivs[which_cfunc] = 0; - } - if (grad != nil) { - if (cfuncs.gderivs[which_cfunc] > 0) { - result = (double *) values[1]; - for (i = 0; i 0)) - add_tilt(x, pval, grad, hess, tiltinfo.tilt, tiltinfo.exptilt); -} - -void -constfunc(double *x, double *vals, double **jac, double **hess) -{ - int i, k = cfuncs.k; - double *pvali, *jaci; - - for (i = 0; i < k; i++) { - pvali = (vals != nil) ? vals + i : nil; - jaci = (jac != nil) ? jac[i] : nil; - which_cfunc = i; - evalcfunc(x, pvali, jaci, nil); - } -} - -void -maxfront(char **ff, char **gf, char **cf, - double *x, double *typx, double *fvals, double *gvals, double *cvals, double *ctarget, - MaxIPars *ipars, MaxDPars *dpars, - double *tscale, char **msg) -{ - static char *work = nil; - static double **H = nil, **cJ = nil; - double *pf, *grad, *c; - size_t i, n, m, k; - void (*cfun)(); - - if (ipars->verbose > 0) PRINTSTR("maximizing...\n"); - - n = ipars->n; - m = ipars->m; - k = ipars->k; - if (k >= n) Recover("too many constraints", NULL); - - makespace(&H, n * sizeof(double *)); - makespace(&work, minworkspacesize(n, k)); - - pf = fvals; fvals++; - grad = fvals; fvals += n; - for (i = 0; i < n; i++, fvals += n) H[i] = fvals; - set_tilt_info(n, m, dpars->newtilt, ipars->exptilt, tscale); - - if (k == 0) { - c = nil; - cJ = nil; - cfun = nil; - } else { - c = cvals; - cvals += k; - makespace(&cJ, k * sizeof(double *)); - for (i = 0; i < k; i++) cJ[i] = cvals + i * n; - cfun = &constfunc; /* AJR: pointer to constfunc? */ - } - - install_func(ff[0], nil, n, TRUE, dpars->typf, dpars->h, typx, dpars->dflt); - install_gfuncs(gf, n, m, TRUE, dpars->h, typx); - install_cfuncs(cf, n, k, ctarget, dpars->h, typx); - - minsetup(n, k, &minfunc, &cfun, x, dpars->typf, typx, work); /* AJR: FIXME arg 3,4 */ - minsetoptions(dpars->gradtol, dpars->steptol, dpars->maxstep, - ipars->itnlimit, ipars->verbose, ipars->backtrack, TRUE); - - if (ipars->vals_suppl) { - for (i = 0; i < k; i++) c[i] -= ctarget[i]; - if (dpars->newtilt != dpars->tilt) { - add_tilt(x, pf, grad, H, dpars->newtilt - dpars->tilt, ipars->exptilt); - dpars->tilt = dpars->newtilt; - } - minsupplyvalues(*pf, grad, H, c, cJ); - } - - minimize(); - minresults(x, pf, nil, grad, H, c, cJ, &ipars->count, &ipars->termcode, - &dpars->hessadd); - msg[0] = minresultstring(ipars->termcode); - - for (i = 0; i < k; i++) c[i] += ctarget[i]; - ipars->vals_suppl = TRUE; -} - -/************************************************************************/ -/** **/ -/** Log Laplace Approximation **/ -/** **/ -/************************************************************************/ - -void -loglapdet(double *fvals, double *cvals, MaxIPars *ipars, MaxDPars *dpars, - double *val, int *detonly) -{ - int i, j, l, n = ipars->n, k = ipars->k; - double f = -fvals[0], *hessdata = fvals + n + 1, *cgraddata = cvals + k; - double ldL, ldcv, maxadd; - static double **hess = nil, **cgrad = nil; - - if (k >= n) Recover("too many constraints", NULL); - - makespace(&hess, n * sizeof(double *)); - makespace(&cgrad, k * sizeof(double *)); - - for (i = 0; i < n; i++) hess[i] = hessdata + i * n; - for (i = 0; i < k; i++) cgrad[i] = cgraddata + i * n; - - choldecomp(hess, n, 0.0, &maxadd); - /**** do something if not pos. definite ****/ - - for (i = 0, ldL = 0.0; i < n; i++) ldL += log(hess[i][i]); - - if (k > 0) { - /* forward solve for (L^-1) cgrad^T */ - for (l = 0; l < k; l++) { - for (i = 0; i < n; i++) { - if (hess[i][i] != 0.0) cgrad[l][i] /= hess[i][i]; - for (j = i + 1; j < n; j++) cgrad[l][j] -= hess[j][i] * cgrad[l][i]; - } - } - - /* compute sigma and stdev */ - for (i = 0; i < k; i++) { - for (j = i; j < k; j++) { - for (l = 0, hess[i][j] = 0.0; l < n; l++) - hess[i][j] += cgrad[i][l] * cgrad[j][l]; - hess[j][i] = hess[i][j]; - } - } - - choldecomp(hess, k, 0.0, &maxadd); - /**** do something if not pos. definite ****/ - for (i = 0, ldcv = 0.0; i < k; i++) ldcv += log(hess[i][i]); - } - else ldcv = 0.0; - - *val = (n - k) * log(ROOT2PI) - ldL - ldcv; - if (! *detonly) *val += f; -} - -#ifdef SBAYES - -loglapfront(fvals, cvals, ipars, dpars, val) - double *fvals, *cvals; - MaxIPars *ipars; - MaxDPars *dpars; - double *val; -{ - int detonly = FALSE; - - loglapdet(fvals, cvals, ipars, dpars, val, &detonly); -} - -/************************************************************************/ -/** **/ -/** First Order Moments **/ -/** **/ -/************************************************************************/ - -moms1front(char *gf, int *n, int *m, double *x, double *hessdata, double *mean, - double *stdev, double *sigmadata, double *h, double *typx) -{ - int i, j, k; - double *hess, *sigma, *delg; - double *delgdata, maxadd; - - hess = (double **) S_alloc(*n, sizeof(double *)); - sigma = (double **) S_alloc(*m, sizeof(double *)); - delg = (double **) S_alloc(*m, sizeof(double *)); - delgdata = (double *) S_alloc(*m * *n, sizeof(double)); - - for (i = 0; i < *n; i++) hess[i] = hessdata + i * *n; - for (i = 0; i < *m; i++) sigma[i] = sigmadata + i * *m; - for (i = 0; i < *m; i++) delg[i] = delgdata + i * *n; - - gevalfront(gf, n, m, x, h, typx, mean, delgdata); - - /* get the cholesky decomposition L of the hessian */ - choldecomp(hess, *n, 0.0, &maxadd); - - /* forward solve for (L^-1) delg^T */ - for (k = 0; k < *m; k++) { - for (i = 0; i < *n; i++) { - if (hess[i][i] != 0.0) delg[k][i] /= hess[i][i]; - for (j = i + 1; j < *n; j++) delg[k][j] -= hess[j][i] * delg[k][i]; - } - } - - /* compute sigma and stdev */ - for (i = 0; i < *m; i++) { - for (j = i; j < *m; j++) { - for (k = 0, sigma[i][j] = 0.0; k < *n; k++) - sigma[i][j] += delg[i][k] * delg[j][k]; - sigma[j][i] = sigma[i][j]; - } - stdev[i] = sqrt(sigma[i][i]); - } -} - -/************************************************************************/ -/** **/ -/** Second Order Moments **/ -/** **/ -/************************************************************************/ - -typedef struct { - MaxIPars max; - int full, covar; -} MomIPars; - -typedef struct { - MaxDPars max; - double mgfdel; -} MomDPars; - -moms2front(char **ff, char **gf, int *gnum, double *x, double *typx, - double *fvals, double *gvals, MomIPars *ipars, MomDPars *dpars, - double *mean, double *stdev, double *sigmadata) -{ - char *msg; - double h, loglap0, loglap1, loglap2; - double *tilts, *fvals1, *gvals1, *x1; - MomDPars dp1, *dpars1 = &dp1; - MomIPars ip1, *ipars1 = &ip1; - int i, n, m; - - n = ipars->max.n; - m = *gnum; - h = dpars->max.h; - - tilts = (double *) S_alloc(m, sizeof(double)); - x1 = (double *) S_alloc(n, sizeof(double)); - fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double)); - gvals1 = (double *) S_alloc(m + n * m, sizeof(double)); - - maxfront(ff, nil, nil, x, typx, fvals, gvals, nil, nil, - ipars, dpars, nil, &msg); - copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap0); - copypars(x, fvals, gvals, ipars, dpars, x1, fvals1, gvals1, ipars1, dpars1); - moms1front(gf, &n, &m, x1, fvals1 + n + 1, mean, stdev, sigmadata, &h, typx); - - if (ipars->full) { - for (i = 0; i < m; i++) { - copypars(x, fvals, gvals, ipars, dpars, - x1, fvals1, gvals1, ipars1, dpars1); - ipars1->max.m = 1; - ipars1->max.exptilt = FALSE; - dpars1->max.newtilt = 1.0; - maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil, - ipars1, dpars1, nil, &msg); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap1); - loglap1 -= loglap0; - - copypars(x, fvals, gvals, ipars, dpars, - x1, fvals1, gvals1, ipars1, dpars1); - ipars1->max.m = 1; - ipars1->max.exptilt = FALSE; - dpars1->max.newtilt = 2.0; - maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil, - ipars1, dpars1, nil, &msg); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap2); - loglap2 -= loglap0; - - mean[i] = exp(loglap1); - stdev[i] = sqrt(exp(loglap2) - exp(2.0 * loglap1)); - if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i]; - } - if (ipars->covar) { - char *cgf[2]; - int j; - - for (i = 0; i < m; i++) { - for (j = i + 1; j < m; j++) { - cgf[0] = gf[i]; - cgf[1] = gf[j]; - copypars(x, fvals, gvals, ipars, dpars, - x1, fvals1, gvals1, ipars1, dpars1); - ipars1->max.m = 2; - ipars1->max.exptilt = FALSE; - dpars1->max.newtilt = 1.0; - maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil, - ipars1, dpars1, nil, &msg); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap1); - loglap1 -= loglap0; - - sigmadata[i + j * m] = exp(loglap1) - mean[i] * mean[j]; - sigmadata[j + i * m] = sigmadata[i + j * m]; - } - } - } - } - else { - for (i = 0; i < m; i++) - tilts[i] = (stdev[i] > 0.0) ? dpars->mgfdel / stdev[i] : dpars->mgfdel; - - for (i = 0; i < m; i++) { - copypars(x, fvals, gvals, ipars, dpars, - x1, fvals1, gvals1, ipars1, dpars1); - ipars1->max.m = 1; - ipars1->max.exptilt = TRUE; - dpars1->max.newtilt = tilts[i]; - maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil, - ipars1, dpars1, nil, &msg); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap1); - loglap1 -= loglap0; - - copypars(x, fvals, gvals, ipars, dpars, - x1, fvals1, gvals1, ipars1, dpars1); - ipars1->max.m = 1; - ipars1->max.exptilt = TRUE; - dpars1->max.newtilt = -tilts[i]; - maxfront(ff, gf + i, nil, x1, typx, fvals1, gvals1, nil, nil, - ipars1, dpars1, nil, &msg); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap2); - loglap2 -= loglap0; - - mean[i] = (loglap1 - loglap2) / (2.0 * tilts[i]); - stdev[i] = sqrt((loglap1 + loglap2) / (tilts[i] * tilts[i])); - if (ipars->covar) sigmadata[i + i * m] = stdev[i] * stdev[i]; - } - if (ipars->covar) { - char *cgf[2]; - double ctilt, tscale[2]; - int j; - - ctilt = dpars->mgfdel; - for (i = 0; i < m; i++) { - for (j = i + 1; j < m; j++) { - cgf[0] = gf[i]; - cgf[1] = gf[j]; - tscale[0] = stdev[i] > 0 ? stdev[i] : 1.0; - tscale[1] = stdev[j] > 0 ? stdev[j] : 1.0; - - copypars(x, fvals, gvals, ipars, dpars, - x1, fvals1, gvals1, ipars1, dpars1); - ipars1->max.m = 2; - ipars1->max.exptilt = TRUE; - dpars1->max.newtilt = ctilt; - maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil, - ipars1, dpars1, tscale, &msg); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap1); - loglap1 -= loglap0; - - copypars(x, fvals, gvals, ipars, dpars, - x1, fvals1, gvals1, ipars1, dpars1); - ipars1->max.m = 2; - ipars1->max.exptilt = TRUE; - dpars1->max.newtilt = -ctilt; - maxfront(ff, cgf, nil, x1, typx, fvals1, gvals1, nil, nil, - ipars1, dpars1, tscale, &msg); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap2); - loglap2 -= loglap0; - - sigmadata[i + j * m] = stdev[i] * stdev[j] - * ((loglap2 + loglap1) /(2 * ctilt * ctilt) - 1.0); - sigmadata[j + i * m] = sigmadata[i + j * m]; - } - } - } - } -} - -static copypars(double *x, double *f, double *g, MomIPars *ip, MomIPars *dp, - double *x1, double *f1, double *g1, MomIPars *ip1, MomDPars *dp1) -{ - int i, n, m, nf, ng; - - n = ip->max.n; - m = ip->max.m; - nf = 1 + n + n * n; - ng = m + n * m; - - for (i = 0; i < n; i++) x1[i] = x[i]; - for (i = 0; i < nf; i++) f1[i] = f[i]; - for (i = 0; i < ng; i++) g1[i] = g[i]; - *ip1 = *ip; - *dp1 = *dp; -} - -/************************************************************************/ -/** **/ -/** Laplace Margins **/ -/** **/ -/************************************************************************/ - -lapmar1front(char **ff, char **gf, double *x, double *typx, double *fvals, - MaxIPars *ipars, MaxDPars *dpars, double *xmar, double *ymar, - int *nmar) -{ - char *msg; - int i, n, m, mindex; - double h, loglap0, loglap1, xmode, stdev, sigmadata, ctarget[1]; - double *fvals1, *x1, *cvals, *cvals1, *fvals2, *x2, *cvals2; - MaxDPars dp1, dp2, *dpars1 = &dp1, *dpars2 = &dp2; - MaxIPars ip1, ip2, *ipars1 = &ip1, *ipars2 = &ip2; - - n = ipars->n; - m = 1; - h = dpars->h; - - x1 = (double *) S_alloc(n + 1, sizeof(double)); - fvals1 = (double *) S_alloc(1 + n + n * n, sizeof(double)); - cvals = (double *) S_alloc(m + n * m, sizeof(double)); - cvals1 = (double *) S_alloc(m + n * m, sizeof(double)); - x2 = (double *) S_alloc(n + 1, sizeof(double)); - fvals2 = (double *) S_alloc(1 + n + n * n, sizeof(double)); - cvals2 = (double *) S_alloc(m + n * m, sizeof(double)); - - maxfront(ff, nil, nil, x, typx, fvals, nil, nil, nil, - ipars, dpars, nil, &msg); - cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1); - loglapfront(fvals1, nil, ipars1, dpars1, &loglap0); - cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1); - moms1front(gf, &n, &m, x1, fvals1 + n + 1, &xmode, &stdev, &sigmadata, - &h, typx); - - ipars->k = 1; - ipars->vals_suppl = FALSE; - ctarget[0] = xmode; - maxfront(ff, nil, gf, x, typx, fvals, nil, cvals, ctarget, - ipars, dpars, nil, &msg); - - for (mindex = 0; mindex < *nmar && xmar[mindex] <= xmode; mindex++); - - cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1); - for (i = mindex; i >= 0; i--) { - ctarget[0] = xmar[i]; - maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget, - ipars1, dpars1, nil, &msg); - cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2, - fvals2, cvals2, ipars2, dpars2); - loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1); - ymar[i] = exp(loglap1 - loglap0); - } - cpmarpars(x, fvals, cvals, ipars, dpars, x1, fvals1, cvals1, ipars1, dpars1); - for (i = mindex + 1; i < *nmar; i++) { - ctarget[0] = xmar[i]; - maxfront(ff, nil, gf, x1, typx, fvals1, nil, cvals1, ctarget, - ipars1, dpars1, nil, &msg); - cpmarpars(x1, fvals1, cvals1, ipars1, dpars1, x2, - fvals2, cvals2, ipars2, dpars2); - loglapfront(fvals2, cvals2, ipars2, dpars2, &loglap1); - ymar[i] = exp(loglap1 - loglap0); - } -} - -static cpmarpars(double *x, double *f, double *g, MaxIPars *ip, MaxDPars *dp, - double *x1, double *f1, double *g1, MaxIPars *ip1, MaxDPars *dp1) -{ - int i, n, k, nf, ng; - - n = ip->n; - k = ip->k; - nf = 1 + n + n * n; - ng = k + n * k; - - for (i = 0; i < n; i++) x1[i] = x[i]; - for (i = 0; i < nf; i++) f1[i] = f[i]; - for (i = 0; i < ng; i++) g1[i] = g[i]; - *ip1 = *ip; - *dp1 = *dp; -} -#endif /* SBAYES */ - -/* - TODO - - get hessian from gradiant for analytical gradiants - - avoid repeated derivative calls in mimimize. - - 2d margins - - use pos. definiteness info in margins - -*/ diff --git a/lib/gamln.c b/lib/gamln.c deleted file mode 100644 index 7e5510d..0000000 --- a/lib/gamln.c +++ /dev/null @@ -1,35 +0,0 @@ -#include "xmath.h" - -/* -log gamma function from Numerical Recipes -*/ - -static double cof[6] = { - 76.18009173, - -86.50532033, - 24.01409822, - -1.231739516, - 0.120858003e-2, - -0.536382e-5 -}; - -double gamma(xx) - double xx; -{ - double x, tmp, ser; - int j; - - if (xx < 1.0) return(gamma(1.0 + xx) - log(xx)); - else { - x = xx - 1.0; - tmp = x + 5.5; - tmp -= (x + 0.5) * log(tmp); - ser = 1.0; - for (j = 0; j < 6; j++) { - x += 1.0; - ser += cof[j] / x; - } - return(-tmp + log(2.50662827465 * ser)); - } -} - diff --git a/lib/gammabase.c b/lib/gammabase.c deleted file mode 100644 index 66f833e..0000000 --- a/lib/gammabase.c +++ /dev/null @@ -1,289 +0,0 @@ -#include "xmath.h" - -#define TRUE 1 -#define FALSE 0 - -extern void normbase(double *,double *); -extern double gamma(), ppnd(); -static double gammp(), gser(), gcf(), gnorm(), ppchi2(); - - -void -gammabase(double *x, double *a, double *p) -{ - *p = gammp(*a, *x); -} - -double -ppgamma(double p, double a, int *ifault) -{ - double x, v, g; - - v = 2.0 * a; - g = gamma(a); - x = ppchi2(&p, &v, &g, ifault); - return(x / 2.0); -} - -/* - Static Routines -*/ - -/* - From Numerical Recipes, with normal approximation from Appl. Stat. 239 -*/ - -#define EPSILON 1.0e-14 -#define LARGE_A 10000.0 -#define ITMAX 1000 - -static double gammp(a, x) - double a, x; -{ - double gln, p; - - if (x <= 0.0 || a <= 0.0) p = 0.0; - else if (a > LARGE_A) p = gnorm(a, x); - else { - gln = gamma(a); - if (x < (a + 1.0)) p = gser(a, x, gln); - else p = 1.0 - gcf(a, x, gln); - } - return(p); -} - -/* compute gamma cdf by a normal approximation */ -static double gnorm(a, x) - double a, x; -{ - double p, sx; - - if (x <= 0.0 || a <= 0.0) p = 0.0; - else { - sx = sqrt(a) * 3.0 * (pow(x / a, 1.0 / 3.0) + 1.0 / (a * 9.0) - 1.0); - normbase(&sx, &p); - } - return(p); -} - -/* compute gamma cdf by its series representation */ -static double gser(a, x, gln) - double a, x, gln; -{ - double p, sum, del, ap; - int n, done = FALSE; - - if (x <= 0.0 || a <= 0.0) p = 0.0; - else { - ap = a; - del = sum = 1.0 / a; - for (n = 1; ! done && n < ITMAX; n++) { - ap += 1.0; - del *= x / ap; - sum += del; - if (fabs(del) < EPSILON) done = TRUE; - } - p = sum * exp(- x + a * log(x) - gln); - } - return(p); -} - -/* compute complementary gamma cdf by its continued fraction expansion */ -static double gcf(a, x, gln) - double a, x, gln; -{ - double gold = 0.0, g, fac = 1.0, b1 = 1.0; - double b0 = 0.0, anf, ana, an, a1, a0 = 1.0; - double p; - int done = FALSE; - - a1 = x; - p = 0.0; - for(an = 1.0; ! done && an <= ITMAX; an += 1.0) { - ana = an - a; - a0 = (a1 + a0 * ana) * fac; - b0 = (b1 + b0 * ana) * fac; - anf = an * fac; - a1 = x * a0 + anf * a1; - b1 = x * b0 + anf * b1; - if (a1 != 0.0) { - fac = 1.0 / a1; - g = b1 * fac; - if (fabs((g - gold) / g) < EPSILON) { - p = exp(-x + a * log(x) - gln) * g; - done = TRUE; - } - gold = g; - } - } - return(p); -} - -static double gammad(x, a, iflag) - double *x, *a; - int *iflag; -{ - double cdf; - - gammabase(x, a, &cdf); - return(cdf); -} - -/* - ppchi2.f -- translated by f2c and modified - - Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35 - To evaluate the percentage points of the chi-squared - probability distribution function. - - p must lie in the range 0.000002 to 0.999998, - (but I am using it for 0 < p < 1 - seems to work) - v must be positive, - g must be supplied and should be equal to ln(gamma(v/2.0)) - - Auxiliary routines required: ppnd = AS 111 (or AS 241) and gammad. -*/ - -static double ppchi2(p, v, g, ifault) - double *p, *v, *g; - int *ifault; -{ - /* Initialized data */ - - static double aa = .6931471806; - static double six = 6.; - static double c1 = .01; - static double c2 = .222222; - static double c3 = .32; - static double c4 = .4; - static double c5 = 1.24; - static double c6 = 2.2; - static double c7 = 4.67; - static double c8 = 6.66; - static double c9 = 6.73; - static double e = 5e-7; - static double c10 = 13.32; - static double c11 = 60.; - static double c12 = 70.; - static double c13 = 84.; - static double c14 = 105.; - static double c15 = 120.; - static double c16 = 127.; - static double c17 = 140.; - static double c18 = 1175.; - static double c19 = 210.; - static double c20 = 252.; - static double c21 = 2264.; - static double c22 = 294.; - static double c23 = 346.; - static double c24 = 420.; - static double c25 = 462.; - static double c26 = 606.; - static double c27 = 672.; - static double c28 = 707.; - static double c29 = 735.; - static double c30 = 889.; - static double c31 = 932.; - static double c32 = 966.; - static double c33 = 1141.; - static double c34 = 1182.; - static double c35 = 1278.; - static double c36 = 1740.; - static double c37 = 2520.; - static double c38 = 5040.; - static double zero = 0.; - static double half = .5; - static double one = 1.; - static double two = 2.; - static double three = 3.; - -/* - static double pmin = 2e-6; - static double pmax = .999998; -*/ - static double pmin = 0.0; - static double pmax = 1.0; - - /* System generated locals */ - double ret_val, d_1, d_2; - - /* Local variables */ - static double a, b, c, q, t, x, p1, p2, s1, s2, s3, s4, s5, s6, ch; - static double xx; - static int if1; - - - /* test arguments and initialise */ - ret_val = -one; - *ifault = 1; - if (*p <= pmin || *p >= pmax) return ret_val; - *ifault = 2; - if (*v <= zero) return ret_val; - *ifault = 0; - xx = half * *v; - c = xx - one; - - if (*v < -c5 * log(*p)) { - /* starting approximation for small chi-squared */ - ch = pow(*p * xx * exp(*g + xx * aa), one / xx); - if (ch < e) { - ret_val = ch; - return ret_val; - } - } - else if (*v > c3) { - /* call to algorithm AS 111 - note that p has been tested above. */ - /* AS 241 could be used as an alternative. */ - x = ppnd(*p, &if1); - - /* starting approximation using Wilson and Hilferty estimate */ - p1 = c2 / *v; - /* Computing 3rd power */ - d_1 = x * sqrt(p1) + one - p1; - ch = *v * (d_1 * d_1 * d_1); - - /* starting approximation for p tending to 1 */ - if (ch > c6 * *v + six) - ch = -two * (log(one - *p) - c * log(half * ch) + *g); - } - else{ - /* starting approximation for v less than or equal to 0.32 */ - ch = c4; - a = log(one - *p); - do { - q = ch; - p1 = one + ch * (c7 + ch); - p2 = ch * (c9 + ch * (c8 + ch)); - d_1 = -half + (c7 + two * ch) / p1; - d_2 = (c9 + ch * (c10 + three * ch)) / p2; - t = d_1 - d_2; - ch -= (one - exp(a + *g + half * ch + c * aa) * p2 / p1) / t; - } while (fabs(q / ch - one) > c1); - } - - do { - /* call to gammad and calculation of seven term Taylor series */ - q = ch; - p1 = half * ch; - p2 = *p - gammad(&p1, &xx, &if1); - if (if1 != 0) { - *ifault = 3; - return ret_val; - } - t = p2 * exp(xx * aa + *g + p1 - c * log(ch)); - b = t / ch; - a = half * t - b * c; - s1 = (c19 + a * (c17 + a * (c14 + a * (c13 + a * (c12 + c11 * a))))) / c24; - s2 = (c24 + a * (c29 + a * (c32 + a * (c33 + c35 * a)))) / c37; - s3 = (c19 + a * (c25 + a * (c28 + c31 * a))) / c37; - s4 = (c20 + a * (c27 + c34 * a) + c * (c22 + a * (c30 + c36 * a))) / c38; - s5 = (c13 + c21 * a + c * (c18 + c26 * a)) / c37; - s6 = (c15 + c * (c23 + c16 * c)) / c38; - d_1 = (s3 - b * (s4 - b * (s5 - b * s6))); - d_1 = (s1 - b * (s2 - b * d_1)); - ch += t * (one + half * t * s1 - b * c * d_1); - } while (fabs(q / ch - one) > e); - - ret_val = ch; - return ret_val; -} diff --git a/lib/kernel.c b/lib/kernel.c deleted file mode 100644 index 5504a47..0000000 --- a/lib/kernel.c +++ /dev/null @@ -1,91 +0,0 @@ -#include "xmath.h" - -#ifndef ROOT2PI -#define ROOT2PI 2.50662827463100050241 -#endif /* ROOT2PI*/ - -#ifndef nil -#define nil 0L -#endif /*nil */ - -/* -#ifdef __APPLE__ -#include -#endif -*/ - -#include - -static double -kernel(double x, double y, double w, int type) -{ - double z, k; - - if (w > 0.0) { - z = (x - y) / w; - switch (type) { - case 'B': - w = 2.0 * w; - z = 0.5 * z; - if (-0.5 < z && z < 0.5) { - z = (1.0 - 4 * z * z); - k = 15.0 * z * z / 8.0; - } - else k = 0.0; - break; - case 'G': - w = 0.25 * w; - z = 4.0 * z; - k = exp(- 0.5 * z * z) / ROOT2PI; - break; - case 'U': - w = 1.5 * w; - z = .75 * z; - k = (fabs(z) < 0.5) ? 1.0 : 0.0; - break; - case 'T': - if (-1.0 <= z && z < 0.0) k = 1.0 + z; - else if (0.0 <= z && z < 1.0) k = 1.0 - z; - else k = 0.0; - break; - default: k = 0.0; break; - } - k = k / w; - } - else k = 0.0; - - return(k); -} - -int -kernel_smooth(double *x, double *y, size_t n, - double width, double *wts, double *wds, - double *xs, double *ys, size_t ns, int ktype) -{ - size_t i, j; - double wsum, ysum, lwidth, lwt, xmin, xmax; - - if (n < 1) return(1); - if (width <= 0.0) { - if (n < 2) return(1); - for (xmin = xmax = x[0], i = 1; i < n; i++) { - if (xmin > x[i]) xmin = x[i]; - if (xmax < x[i]) xmax = x[i]; - } - width = (xmax - xmin) / (1 + log((double) n)); - } - - for (i = 0; i < ns; i++) { - for (j = 0, wsum = 0.0, ysum = 0.0; j < n; j++) { - lwidth = (wds != nil) ? width * wds[j] : width; - lwt = kernel(xs[i], x[j], lwidth, ktype); - if (wts != nil) lwt *= wts[j]; - wsum += lwt; - if (y != nil) ysum += lwt * y[j]; - } - if (y != nil) ys[i] = (wsum > 0.0) ? ysum / wsum : 0.0; - else ys[i] = wsum / n; - } - return(0); -} - diff --git a/lib/linalg.h b/lib/linalg.h deleted file mode 100644 index f8e65c9..0000000 --- a/lib/linalg.h +++ /dev/null @@ -1,27 +0,0 @@ -# include "xmath.h" -# include "xlisp.h" -# include "complex.h" - -extern void *calloc(); /*conflict? was (char *) */ -extern double macheps(); - -#define nil 0L - -typedef char **Matrix, *Vector; -typedef int **IMatrix, *IVector; -typedef double **RMatrix, *RVector; -typedef Complex **CMatrix, *CVector; - -#define IN 0 -#define RE 1 -#define CX 2 - -/* external functions */ -extern IVector ivector(); -extern RVector rvector(); -extern CVector cvector(); -extern IMatrix imatrix(); -extern RMatrix rmatrix(); -extern CMatrix cmatrix(); - -extern void free_vector(double *); diff --git a/lib/linalgdata.c b/lib/linalgdata.c deleted file mode 100644 index db1bcc7..0000000 --- a/lib/linalgdata.c +++ /dev/null @@ -1,92 +0,0 @@ -/* linalgdata - allocation support for basic linear algebra routines. */ -/* Copyright (c) 1990, by Luke Tierney */ - -#include "linalg.h" - -#ifdef INTPTR -typedef int PTR; -#else -typedef char *PTR; -#endif - -extern PTR la_allocate(size_t, size_t); -extern void la_free_alloc(PTR); -extern void xlfail(char *); - -/************************************************************************/ -/** **/ -/** Storage Allocation Functions **/ -/** **/ -/************************************************************************/ - -static char -*allocate(size_t n, size_t m) -{ - char *p = (char *) la_allocate(n, m); - if (p == nil) xlfail("allocation failed"); - return(p); -} - -static void -free_alloc(char *p) -{ - if (p != nil) la_free_alloc((PTR) p); -} - -IVector -ivector(size_t n) -{ - return((IVector) allocate(n, sizeof(int))); -} - -double -*rvector(size_t n) -{ - return((double *) allocate(n, sizeof(double))); -} - -CVector cvector(size_t n) -{ - return((CVector) allocate(n, sizeof(Complex))); -} - -void -free_vector(double *v) -{ - free_alloc((char *)v); -} - -IMatrix imatrix(size_t n, size_t m) -{ - size_t i; - IMatrix mat = (IMatrix) allocate(n, sizeof(IVector)); - for (i = 0; i < n; i++) mat[i] = (IVector) allocate(m, sizeof(int)); - return(mat); -} - -RMatrix -rmatrix(size_t n, size_t m) -{ - size_t i; - RMatrix mat = (RMatrix) allocate(n, sizeof(RVector)); - for (i = 0; i < n; i++) mat[i] = (RVector) allocate(m, sizeof(double)); - return(mat); -} - -CMatrix -cmatrix(size_t n, size_t m) -{ - size_t i; - CMatrix mat = (CMatrix) allocate(n, sizeof(CVector)); - for (i = 0; i < n; i++) mat[i] = (CVector) allocate(m, sizeof(Complex)); - return(mat); -} - -void -free_matrix(double **mat, int n) /* Matrix?? Not RMatrix?? */ -{ - size_t i; - - if (mat != nil) for (i = 0; i < n; i++) free_alloc((char *)mat[i]); - free_alloc((char *)mat); -} diff --git a/lib/lowess.c b/lib/lowess.c deleted file mode 100644 index c0780f2..0000000 --- a/lib/lowess.c +++ /dev/null @@ -1,170 +0,0 @@ -/* - Translated from RATFOR lowess code of W. S. Cleveland as obtained from NETLIB -*/ - -#include "xmath.h" -/* -#ifdef __APPLE__ -#include -#endif -*/ -#include - -#define FALSE 0 -#define TRUE 1 - -extern long max(long, long); -extern long min(long, long); - -static double pow2(double x) { return(x * x); } -static double pow3(double x) { return(x * x * x); } -/* static */ double fmax(double x, double y) { return (x > y ? x : y); } - -int -static compar(const void *aa, const void *bb) -{ - const double* a = (double*)aa; - const double* b = (double*)bb; - - if (*a < *b) return(-1); - else if (*a > *b) return(1); - else return(0); -} - -static void -lowest(double *x, double *y, size_t n, double xs, double *ys, long nleft, long nright, - double *w, int userw, double *rw, int *ok) -{ - double range, h, h1, h9, a, b, c, r; - long j, nrt; - - range = x[n - 1] - x[0]; - h = fmax(xs - x[nleft], x[nright] - xs); - h9 = .999 * h; - h1 = .001 * h; - - /* compute weights (pick up all ties on right) */ - a = 0.0; /* sum of weights */ - for(j = nleft; j < n; j++) { - w[j]=0.0; - r = fabs(x[j] - xs); - if (r <= h9) { /* small enough for non-zero weight */ - if (r > h1) w[j] = pow3(1.0-pow3(r/h)); - else w[j] = 1.0; - if (userw) w[j] = rw[j] * w[j]; - a += w[j]; - } - else if (x[j] > xs) break; /* get out at first zero wt on right */ - } - nrt = j - 1; /* rightmost pt (may be greater than nright because of ties) */ - if (a <= 0.0) *ok = FALSE; - else { /* weighted least squares */ - *ok = TRUE; - - /* make sum of w[j] == 1 */ - for (j = nleft; j <= nrt; j++) w[j] = w[j] / a; - - if (h > 0.0) { /* use linear fit */ - - /* find weighted center of x values */ - for (j = nleft, a = 0.0; j <= nrt; j++) a += w[j] * x[j]; - - b = xs - a; - for (j = nleft, c = 0.0; j <= nrt; j++) - c += w[j] * (x[j] - a) * (x[j] - a); - - if(sqrt(c) > .001 * range) { - /* points are spread out enough to compute slope */ - b = b/c; - for (j = nleft; j <= nrt; j++) - w[j] = w[j] * (1.0 + b*(x[j] - a)); - } - } - for (j = nleft, *ys = 0.0; j <= nrt; j++) *ys += w[j] * y[j]; - } -} - -static void -sort(double *x, size_t n) -{ - extern void qsort(); - - qsort(x, n, sizeof(double), compar); -} - - - -int -lowess(double *x, double *y, size_t n, - double f, size_t nsteps, - double delta, double *ys, double *rw, double *res) -{ - int iter, ok; - long i, j, last, m1, m2, nleft, nright, ns; - double d1, d2, denom, alpha, cut, cmad, c9, c1, r; - - if (n < 2) { ys[0] = y[0]; return(1); } - ns = max(min((long) (f * n), n), 2); /* at least two, at most n points */ - for(iter = 1; iter <= nsteps + 1; iter++){ /* robustness iterations */ - nleft = 0; nright = ns - 1; - last = -1; /* index of prev estimated point */ - i = 0; /* index of current point */ - do { - while(nright < n - 1){ - /* move nleft, nright to right if radius decreases */ - d1 = x[i] - x[nleft]; - d2 = x[nright + 1] - x[i]; - /* if d1 <= d2 with x[nright+1] == x[nright], lowest fixes */ - if (d1 <= d2) break; - /* radius will not decrease by move right */ - nleft++; - nright++; - } - lowest(x, y, - n, x[i], - &ys[i], - nleft, nright, - res, (iter > 1), rw, &ok); - /* fitted value at x[i] */ - if (! ok) ys[i] = y[i]; - /* all weights zero - copy over value (all rw==0) */ - if (last < i - 1) { /* skipped points -- interpolate */ - denom = x[i] - x[last]; /* non-zero - proof? */ - for(j = last + 1; j < i; j = j + 1){ - alpha = (x[j] - x[last]) / denom; - ys[j] = alpha * ys[i] + (1.0 - alpha) * ys[last]; - } - } - last = i; /* last point actually estimated */ - cut = x[last] + delta; /* x coord of close points */ - for(i=last + 1; i < n; i++) { /* find close points */ - if (x[i] > cut) break; /* i one beyond last pt within cut */ - if(x[i] == x[last]) { /* exact match in x */ - ys[i] = ys[last]; - last = i; - } - } - i = max(last + 1,i - 1); - /* back 1 point so interpolation within delta, but always go forward */ - } while(last < n - 1); - for (i = 0; i < n; i++) /* residuals */ - res[i] = y[i] - ys[i]; - if (iter > nsteps) break; /* compute robustness weights except last time */ - for (i = 0; i < n; i++) - rw[i] = fabs(res[i]); - sort(rw,n); - m1 = 1 + n / 2; m2 = n - m1 + 1; - cmad = 3.0 * (rw[m1] + rw[m2]); /* 6 median abs resid */ - c9 = .999 * cmad; c1 = .001 * cmad; - for (i = 0; i < n; i++) { - r = fabs(res[i]); - if(r <= c1) rw[i] = 1.0; /* near 0, avoid underflow */ - else if(r > c9) rw[i] = 0.0; /* near 1, avoid underflow */ - else rw[i] = pow2(1.0 - pow2(r / cmad)); - } - } - return(0); -} - - - diff --git a/lib/ludecomp.c b/lib/ludecomp.c deleted file mode 100644 index e659d1c..0000000 --- a/lib/ludecomp.c +++ /dev/null @@ -1,176 +0,0 @@ -/* ludecomp - LU decomposition and backsolving routines. */ -/* Taken from Numerical Recipies. */ -/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */ -/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */ -/* You may give out copies of this software; for conditions see the */ -/* file COPYING included with this distribution. */ - -#include "linalg.h" - -int -crludcmp(Matrix mat, size_t n, IVector indx, int mode, double *d) -{ - int i, imax, j, k; - int singular = FALSE; - double big, temp; - Complex cdum, csum; - double rdum, rsum; - CMatrix cmat = (CMatrix) mat; - RMatrix rmat = (RMatrix) mat; - RVector vv; - - vv = rvector(n); - *d = 1.0; - - /* set up the pivot permutation vector */ - for (i = 0; i < n; i++) indx[i] = i; - - /* get scaling information for implicit pivoting */ - for (i = 0; i < n; i++) { - big = 0.0; - for (j = 0; j < n; j++) { - temp = (mode == RE) ? fabs(rmat[i][j]) : modulus(cmat[i][j]); - if (temp > big) big = temp; - } - if (big == 0.0) { - vv[i] = 1.0; /* no scaling for zero rows */ - singular = TRUE; - } - else vv[i] = 1.0 / big; - } - - /* loop over columns for Crout's method */ - for (j = 0; j < n; j++) { - for (i = 0; i < j; i++) { - if (mode == RE) rsum = rmat[i][j]; - else csum = cmat[i][j]; - - for (k = 0; k < i; k++) - { - if (mode == RE) rsum -= rmat[i][k] * rmat[k][j]; - else csum = csub(csum, cmul(cmat[i][k], cmat[k][j])); - } - - if (mode == RE) rmat[i][j] = rsum; - else cmat[i][j] = csum; - } - big = 0.0; - for (i = j; i < n; i++) { - if (mode == RE) rsum = rmat[i][j]; - else csum = cmat[i][j]; - - for (k = 0; k < j; k++) - { - if (mode == RE) rsum -= rmat[i][k] * rmat[k][j]; - else csum = csub(csum, cmul(cmat[i][k], cmat[k][j])); - } - - if (mode == RE) rmat[i][j] = rsum; - else cmat[i][j] = csum; - - temp = vv[i] * ((mode == RE) ? fabs(rsum) : modulus(csum)); - if (temp >= big) { big = temp; imax = i; } - } - - /* interchange rows if needed */ - if (j != imax) { - for (k = 0; k < n; k++) { - if (mode == RE) { - rdum = rmat[imax][k]; - rmat[imax][k] = rmat[j][k]; - rmat[j][k] = rdum; - } - else { - cdum = cmat[imax][k]; - cmat[imax][k] = cmat[j][k]; - cmat[j][k] = cdum; - } - } - *d = -(*d); - vv[imax] = vv[j]; - } - indx[j] = imax; - - /* divide by the pivot element */ - temp = (mode == RE) ? fabs(rmat[j][j]) : modulus(cmat[j][j]); - if (temp == 0.0) singular = TRUE; - else if (j < n - 1) { - if (mode == RE) { - rdum = 1.0 / rmat[j][j]; - for (i = j + 1; i < n; i++) rmat[i][j] *= rdum; - } - else { - cdum = cdiv(cart2complex(1.0, 0.0), cmat[j][j]); - for (i = j + 1; i < n; i++) - { - cmat[i][j] = cmul(cmat[i][j], cdum); - } - } - } - } - - free_vector(vv); - return(singular); -} - -int -crlubksb(double **a, int n, int *indx, double *b, int mode) -{ - int i, ii, ip, j, singular = FALSE; - CMatrix ca = (CMatrix) a; - CVector cb = (CVector) b; - RMatrix ra = (RMatrix) a; - RVector rb = (RVector) b; - double rsum; - Complex csum; - - /* forward substitute using L part */ - for (i = 0, ii = -1; i < n; i++) { - ip = indx[i]; - if (mode == RE) { - rsum = rb[ip]; - rb[ip] = rb[i]; - } - else { - csum = cb[ip]; - cb[ip] = cb[i]; - } - if (ii >= 0) - for (j = ii; j <= i - 1; j++) - if (mode == RE) rsum -= ra[i][j] * rb[j]; - else csum = csub(csum, cmul(ca[i][j], cb[j])); - else { - if (mode == RE) { - if (rsum != 0.0) ii = i; - } - else if (csum.real != 0.0 || csum.imag != 0.0) ii = i; - } - if (mode == RE) rb[i] = rsum; - else cb[i] = csum; - } - - /* back substitute using the U part */ - for (i = n - 1; i >= 0; i--) { - if (mode == RE) { - rsum = rb[i]; - for (j = i + 1; j < n; j++) rsum -= ra[i][j] * rb[j]; - if (ra[i][i] == 0.0) { - singular = TRUE; - break; - } - else rb[i] = rsum / ra[i][i]; - } - else { - csum = cb[i]; - for (j = i + 1; j < n; j++) csum = csub(csum, cmul(ca[i][j], cb[j])); - if (modulus(ca[i][i]) == 0.0) { - singular = TRUE; - break; - } - else cb[i] = cdiv(csum, ca[i][i]); - } - } - - return(singular); -} - diff --git a/lib/makerotation.c b/lib/makerotation.c deleted file mode 100644 index 6b7decc..0000000 --- a/lib/makerotation.c +++ /dev/null @@ -1,54 +0,0 @@ -/* makerotation - Construct rotation from x to y by alpha. */ -/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */ -/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */ -/* You may give out copies of this software; for conditions see the */ -/* file COPYING included with this distribution. */ - -#include "linalg.h" - -static double inner_product(size_t n, RVector x, RVector y) -{ - double result = 0.0; - - for (; n > 0; n--, x++, y++) result += *x * *y; - return(result); -} - -#define NORM(n, x) (sqrt(inner_product(n, x, x))) - -void -make_rotation(size_t n, double **rot, double *x, double *y, int use_alpha, double alpha) -{ - double nx, ny, xy, c, s; - size_t i, j; - - for (i = 0; i < n; i++) { - for (j = 0; j < n; j++) rot[i][j] = 0.0; - rot[i][i] = 1.0; - } - - nx = NORM(n, x); - ny = NORM(n, y); - if (nx == 0.0 || ny == 0.0) return; - - for (i = 0; i < n; i++) x[i] /= nx; - for (i = 0; i < n; i++) y[i] /= ny; - - xy = inner_product(n, x, y); - c = (use_alpha) ? cos(alpha) : xy; - s = (use_alpha) ? sin(alpha) : sqrt(1 - c * c); - - for (i = 0; i < n; i++) y[i] -= xy * x[i]; - - ny = NORM(n, y); - if (ny == 0.0) return; - for (i = 0; i < n; i++) y[i] /= ny; - - for (i = 0; i < n; i++) { - for (j = 0; j < n; j++) - rot[i][j] = x[j] * ( x[i] * (c - 1.0) + y[i] * s) - + y[j] * (- x[i] * s + y[i] * (c - 1.0)); - rot[i][i] += 1.0; - } -} - diff --git a/lib/minimize.c b/lib/minimize.c deleted file mode 100644 index 05f5028..0000000 --- a/lib/minimize.c +++ /dev/null @@ -1,868 +0,0 @@ -/* minimize - minimization for Bayes routines in XLISP-STAT and S */ -/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */ -/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */ -/* You may give out copies of this software; for conditions see the */ -/* file COPYING included with this distribution. */ - -/* - * Nonlinear optimization modules adapted from Dennis and Schnabel, - * "Numerical Methods for Unconstrained Optimization and Nonlinear - * Equations." - */ - -#include -#include "xmath.h" -extern char buf[200]; -#define PRINTSTR(s) printf(s) - -typedef struct { - size_t n, k; - void (*ffun)(), (*gfun)(); - double f, typf, new_f; - double crit, new_crit; - double *x, *new_x, *sx, *delf, *new_delf, *qnstep, *F; - double **hessf, **H, **L, **new_delg; - double gradtol, steptol, maxstep; - int itncount, itnlimit, maxtaken, consecmax, retcode, termcode; - int use_line_search, verbose, values_supplied, change_sign; - double diagadd; -} Iteration; - - -extern double macheps(); - - -extern double choldecomp(); -static void eval_gradient(Iteration *); -static void eval_next_gradient(Iteration *); - - -/* -This is a little more advanced section on functions, but is very useful. Take this for example: - int applyeqn(int F(int), int max, int min) { - int itmp; - - itmp = F(int) + min; - itmp = itmp - max; - return itmp; - } - - What does this function do if we call it with applyeqn(square(x), - y, z);? What happens is that the int F(int) is a reference to the - function that is passed in as a parameter. Thus inside applyeqn where - there is a call to F, it actually is a call to square! This is very - useful if we have one set function, but wish to vary the input - according to a particular function. So if we had a different function - called cube we could change how we call applyeqn by calling the - function by applyeqn(cube(x), y, z);. */ - - -/************************************************************************/ -/** **/ -/** Definitions and Globals **/ -/** **/ -/************************************************************************/ - -# define nil 0L -# define FALSE 0 -# define TRUE 1 - -# define INIT_GRAD_FRAC .001 -# define CONSEC_MAX_LIMIT 5 -# define ALPHA .0001 -# define MAX_STEP_FACTOR 1000 -# define GRADTOL_POWER 1.0 / 3.0 -# define STEPTOL_POWER 2.0 / 3.0 -# define ITNLIMIT 100 -# define VERBOSE 0 -# define USE_SEARCH TRUE - -typedef double **RMatrix, *RVector; - - -static char *termcodes[] = {"not yet terminated", - "gradient size is less than gradient tolerance", - "step size is less than step tolerance", - "no satisfactory step found in backtracking", - "iteration limit exceeded", - "maximum size step taken 5 iterations in a row"}; - -/************************************************************************/ -/** **/ -/** Utility Functions **/ -/** **/ -/************************************************************************/ - -static double Max(a, b) - double a, b; -{ - return(a > b ? a : b); -} - -static double Min(a, b) - double a, b; -{ - return(a > b ? b : a); -} - -/************************************************************************/ -/** **/ -/** Cholesky Solving Functions **/ -/** **/ -/************************************************************************/ - - -/* solve Ly = b for y */ -static void -lsolve(size_t n, double *b, double **L, double *y) -{ - size_t i, j; - - for (i = 0; i < n; i++) { - y[i] = b[i]; - for (j = 0; j < i; j++) y[i] -= L[i][j] * y[j]; - if (L[i][i] != 0) y[i] /= L[i][i]; - } -} - -/* solve (L^T)x = y for x */ -static void -ltsolve(size_t n, double *y, double **L, double *x) -{ - size_t i, j; - - for (i = n - 1; i >= 0; i--) { - x[i] = y[i]; - for (j = i + 1; j < n; j++) x[i] -= L[j][i] * x[j]; - if (L[i][i] != 0) x[i] /= L[i][i]; - } -} - -/* solve (L L^T) s = -g for s */ -static void -cholsolve(size_t n, - double *g, double **L, - double *s) -{ - size_t i; - - /* solve Ly = g */ - lsolve(n, g, L, s); - - /* solve L^Ts = y */ - ltsolve(n, s, L, s); - - for (i = 0; i < n; i++) s[i] = -s[i]; -} - -static void -modelhess(size_t n, double *sx, double **H, double **L, double *diagadd) -{ - size_t i, j; - double sqrteps, maxdiag, mindiag, maxoff, maxoffl, maxposdiag, mu, - maxadd, maxev, minev, offrow, sdd; - - /* scale H on both sides with sx */ - for (i = 0; i < n; i++) - for (j = 0; j < n; j++) H[i][j] /= sx[i] * sx[j]; - - /* - * find mu to add to diagonal so no diagonal elements are negative - * and largest diagonal dominates largest off diagonal element in H - */ - sqrteps = sqrt(macheps()); - for (maxdiag = H[0][0], i = 1; i < n; i++) - maxdiag = Max(maxdiag, H[i][i]); - for (mindiag = H[0][0], i = 1; i < n; i++) - mindiag = Min(mindiag, H[i][i]); - maxposdiag = Max(0.0, maxdiag); - - if (mindiag <= sqrteps * maxposdiag) { - mu = 2 * (maxposdiag - mindiag) * sqrteps - mindiag; - maxdiag += mu; - } - else mu = 0.0; - - if (n > 1) { - for (maxoff = fabs(H[0][1]), i = 0; i < n; i++) - for (j = i + 1; j < n; j++) - maxoff = Max(maxoff, fabs(H[i][j])); - } - else maxoff = 0.0; - - if (maxoff * (1 + 2 * sqrteps) > maxdiag) { - mu += (maxoff - maxdiag) + 2 * sqrteps * maxoff; - maxdiag = maxoff * (1 + 2 * sqrteps); - } - - if (maxdiag == 0.0) { - mu = 1; - maxdiag = 1; - } - - if (mu > 0) for (i = 0; i < n; i++) H[i][i] += mu; - - maxoffl = sqrt(Max(maxdiag, maxoff / n)); - - /* - * compute the perturbed Cholesky decomposition - */ - for (i = 0; i < n; i++) - for (j = 0; j < n; j++) L[i][j] = H[i][j]; - choldecomp(L, n, maxoffl, &maxadd); - - /* - * add something to diagonal, if needed, to make positive definite - * and recompute factorization - */ - if (maxadd > 0) { - maxev = H[0][0]; - minev = H[0][0]; - for (i = 0; i < n; i++) { - for (offrow = 0.0, j = 0; j < n; j++) - if (i != j) offrow += fabs(H[i][j]); - maxev = Max(maxev, H[i][i] + offrow); - minev = Min(minev, H[i][i] - offrow); - } - sdd = (maxev - minev) * sqrteps - minev; - sdd = Max(sdd, 0.0); - mu = Min(maxadd, sdd); - for (i = 0; i < n; i++) H[i][i] += mu; - for (i = 0; i < n; i++) - for (j = 0; j < n; j++) L[i][j] = H[i][j]; - choldecomp(L, n, maxoffl, &maxadd); - *diagadd = mu; - } - else *diagadd = 0.0; - - /* unscale H and L */ - for (i = 0; i < n; i++) - for (j = 0; j < n; j++) { - H[i][j] *= sx[i] * sx[j]; - L[i][j] *= sx[i]; - } -} - -/************************************************************************/ -/** **/ -/** Stopping Criteria **/ -/** **/ -/************************************************************************/ - -static double -gradsize(Iteration *iter, int new) -{ - size_t n, i; - double size, term, crit, typf; - double *x, *delf, *sx; - - n = iter->n + iter->k; - crit = iter->crit; - typf = iter->typf; - x = iter->x; - sx = iter->sx; - delf = (new) ? iter->new_delf : iter->delf; - - for (i = 0, size = 0.0; i < n; i++) { - term = fabs(delf[i]) * Max(x[i], 1.0 / sx[i]) / Max(fabs(crit), typf); - size = Max(size, term); - } - return(size); -} - -static double -incrsize(Iteration *iter) -{ - size_t n, i; - double size, term; - double *x, *new_x, *sx; - - new_x = iter->new_x; - n = iter->n + iter->k; - x = iter->x; - sx = iter->sx; - - for (i = 0, size = 0.0; i < n; i++) { - term = fabs(x[i] - new_x[i]) / Max(fabs(x[i]), 1.0 / sx[i]); - size = Max(size, term); - } - return(size); -} - -static int /* AJR:TESTME guess! */ -stoptest0(Iteration *iter) -{ - iter->consecmax = 0; - - if (gradsize(iter, FALSE) <= INIT_GRAD_FRAC * iter->gradtol) - iter->termcode = 1; - else iter->termcode = 0; - - return(iter->termcode); -} - -static int -stoptest(Iteration *iter) -{ - int termcode, retcode, itncount, itnlimit, maxtaken, consecmax; - double gradtol, steptol; - - retcode = iter->retcode; - gradtol = iter->gradtol; - steptol = iter->steptol; - itncount = iter->itncount; - itnlimit = iter->itnlimit; - maxtaken = iter->maxtaken; - consecmax = iter->consecmax; - - termcode = 0; - if (retcode == 1) termcode = 3; - else if (gradsize(iter, TRUE) <= gradtol) termcode = 1; - else if (incrsize(iter) <= steptol) termcode = 2; - else if (itncount >= itnlimit) termcode = 4; - else if (maxtaken) { - consecmax++; - if (consecmax >= CONSEC_MAX_LIMIT) termcode = 5; - } - else consecmax = 0; - - iter->consecmax = consecmax; - iter->termcode = termcode; - - return(termcode); -} - -/************************************************************************/ -/** **/ -/** Function and Derivative Evaluation **/ -/** **/ -/************************************************************************/ - -static void -eval_funval(Iteration *iter) -{ - int i; - - (*(iter->ffun))(iter->x, &iter->f, nil, nil); - if (iter->k == 0) iter->crit = iter->f; - else { - eval_gradient(iter); - for (i = 0, iter->crit = 0.0; i < iter->n + iter->k; i++) - iter->crit += 0.5 * pow(fabs(iter->delf[i]), 2.0); - } -} - -static void -eval_next_funval(Iteration *iter) -{ - int i; - - (*(iter->ffun))(iter->new_x, &iter->new_f, nil, nil); - if (iter->k == 0) iter->new_crit = iter->new_f; - else { - eval_next_gradient(iter); - for (i = 0, iter->new_crit = 0.0; i < iter->n + iter->k; i++) - iter->new_crit += 0.5 * pow(fabs(iter->new_delf[i]), 2.0); - } -} - -static void -eval_gradient(Iteration *iter) -{ - size_t i, j, n, k; - - n = iter->n; - k = iter->k; - - (*(iter->ffun))(iter->x, nil, iter->delf, nil); - if (k > 0) { - (*(iter->gfun))(iter->x, iter->delf + n, nil, nil); - (*(iter->gfun))(iter->x, nil, iter->new_delg, nil); - for (i = 0; i < n; i++) { - for (j = 0; j < k; j++) - iter->delf[i] += iter->x[n + j] * iter->new_delg[j][i]; - for (j = 0; j < k; j++) { - iter->hessf[n + j][i] = iter->new_delg[j][i]; - iter->hessf[i][n + j] = iter->new_delg[j][i]; - } - } - } -} - -static void -eval_next_gradient(Iteration *iter) -{ - size_t i, j, n, k; - - n = iter->n; - k = iter->k; - (*(iter->ffun))(iter->new_x, nil, iter->new_delf, nil); - if (k > 0) { - (*(iter->gfun))(iter->new_x, iter->new_delf + n, nil, nil); - (*(iter->gfun))(iter->new_x, nil, iter->new_delg, nil); - for (i = 0; i < n; i++) { - for (j = 0; j < k; j++) - iter->new_delf[i] += iter->new_x[n + j] * iter->new_delg[j][i]; - } - } -} - -static void eval_hessian(Iteration *iter) -{ - size_t i, j, n, k; - - n = iter->n; - k = iter->k; - (*(iter->ffun))(iter->x, nil, nil, iter->hessf); - for (i = n; i < n + k; i++) - for (j = n; j < n + k; j++) iter->hessf[i][j] = 0.0; -} - -/************************************************************************/ -/** **/ -/** Backtracking Line Search **/ -/** **/ -/************************************************************************/ - -static void -linesearch(Iteration *iter) -{ - size_t i, n; - double newtlen, maxstep, initslope, rellength, lambda, minlambda, - lambdatemp, lambdaprev, a, b, disc, critprev, f1, f2, a11, a12, a21, a22, - del; - double *qnstep, *delf, *x, *new_x, *sx; - - n = iter->n + iter->k; - if (! iter->use_line_search) { - iter->maxtaken = FALSE; - for (i = 0; i < n; i++) - iter->new_x[i] = iter->x[i] + iter->qnstep[i]; - eval_next_funval(iter); - iter->retcode = 0; - } - else{ - qnstep = iter->qnstep; - maxstep = iter->maxstep; - delf = (iter->k == 0) ? iter->delf : iter->F; - x = iter->x; - new_x = iter->new_x; - sx = iter->sx; - - iter->maxtaken = FALSE; - iter->retcode = 2; - - for (i = 0, newtlen = 0.0; i < n; i++) - newtlen += pow(sx[i] * qnstep[i], 2.0); - newtlen = sqrt(newtlen); - - if (newtlen > maxstep) { - for (i = 0; i < n; i++) qnstep[i] *= (maxstep / newtlen); - newtlen = maxstep; - } - - for (i = 0, initslope = 0.0; i < n; i++) initslope += delf[i] * qnstep[i]; - - for (i = 0, rellength = 0.0; i < n; i++) - rellength = Max(rellength, fabs(qnstep[i]) / Max(fabs(x[i]), 1.0 / sx[i])); - - minlambda = (rellength == 0.0) ? 2.0 : iter->steptol / rellength; - - lambda = 1.0; - lambdaprev = 1.0; /* to make compilers happy */ - critprev = 1.0; /* to make compilers happy */ - while (iter->retcode > 1) { - for (i = 0; i < n; i++) new_x[i] = x[i] + lambda * qnstep[i]; - eval_next_funval(iter); - if (iter->new_crit <= iter->crit + ALPHA * lambda * initslope) { - iter->retcode = 0; - if (lambda == 1.0 && newtlen > 0.99 * maxstep) iter->maxtaken = TRUE; - } - else if (lambda < minlambda) { - iter->retcode = 1; - iter->new_crit = iter->crit; - for (i = 0; i < n; i++) new_x[i] = x[i]; - } - else { - if (lambda == 1.0) { /* first backtrack, quadratic fit */ - lambdatemp = - initslope - / (2 * (iter->new_crit - iter->crit - initslope)); - } - else { /* all subsequent backtracks, cubic fit */ - del = lambda - lambdaprev; - f1 = iter->new_crit - iter->crit - lambda * initslope; - f2 = critprev - iter->crit - lambdaprev * initslope; - a11 = 1.0 / (lambda * lambda); - a12 = -1.0 / (lambdaprev * lambdaprev); - a21 = -lambdaprev * a11; - a22 = -lambda * a12; - a = (a11 * f1 + a12 * f2) / del; - b = (a21 * f1 + a22 * f2) / del; - disc = b * b - 3 * a * initslope; - if (a == 0) { /* cubic is a quadratic */ - lambdatemp = - initslope / (2 * b); - } - else { /* legitimate cubic */ - lambdatemp = (-b + sqrt(disc)) / (3 * a); - } - lambdatemp = Min(lambdatemp, 0.5 * lambda); - } - lambdaprev = lambda; - critprev = iter->new_crit; - lambda = Max(0.1 * lambda, lambdatemp); - if (iter->verbose > 0) { - sprintf(buf, "Backtracking: lambda = %g\n", lambda); - PRINTSTR(buf); - } - } - } - } -} - -/************************************************************************/ -/** **/ -/** Status Printing Functions **/ -/** **/ -/************************************************************************/ - -static void -print_header(Iteration *iter) -{ - if (iter->verbose > 0) { - sprintf(buf, "Iteration %d.\n", iter->itncount); - PRINTSTR(buf); - } -} - -static void -print_status(Iteration *iter) -{ - size_t i, j; - - if (iter->verbose > 0) { - sprintf(buf, "Criterion value = %g\n", - (iter->change_sign) ? -iter->crit : iter->crit); - PRINTSTR(buf); - if (iter->verbose > 1) { - PRINTSTR("Location = <"); - for (i = 0; i < iter->n + iter->k; i++) { - sprintf(buf, (i < iter->n + iter->k - 1) ? "%g " : "%g>\n", iter->x[i]); - PRINTSTR(buf); - } - } - if (iter->verbose > 2) { - PRINTSTR("Gradient = <"); - for (i = 0; i < iter->n + iter->k; i++) { - sprintf(buf, (i < iter->n + iter->k - 1) ? "%g " : "%g>\n", - (iter->change_sign) ? -iter->delf[i] : iter->delf[i]); - PRINTSTR(buf); - } - } - if (iter->verbose > 3) { - PRINTSTR("Hessian:\n"); - for (i = 0; i < iter->n + iter->k; i++) { - for (j = 0; j < iter->n + iter->k; j++) { - sprintf(buf, (j < iter->n + iter->k - 1) ? "%g " : "%g\n", - (iter->change_sign) ? -iter->hessf[i][j] : iter->hessf[i][j]); - PRINTSTR(buf); - } - } - } - } - if (iter->termcode != 0 && iter->verbose > 0) { - sprintf(buf, "Reason for termination: %s.\n", termcodes[iter->termcode]); - PRINTSTR(buf); - } -} - -/************************************************************************/ -/** **/ -/** Iteration Driver **/ -/** **/ -/************************************************************************/ - -static void -findqnstep(Iteration *iter) -{ - size_t i, j, N, l; - - if (iter->k == 0) { - modelhess(iter->n, iter->sx, iter->hessf, iter->L, &iter->diagadd); - cholsolve(iter->n, iter->delf, iter->L, iter->qnstep); - } - else { - N = iter->n + iter->k; - for (i = 0; i < N; i++) { - for (l = 0, iter->F[i] = 0.0; l < N; l++) - iter->F[i] += iter->hessf[i][l] * iter->delf[l]; - for (j = 0; j < N; j++) - for (l = 0, iter->H[i][j] = 0.0; l < N; l++) - iter->H[i][j] += iter->hessf[i][l] * iter->hessf[j][l]; - } - modelhess(N, iter->sx, iter->H, iter->L, &iter->diagadd); - cholsolve(N, iter->F, iter->L, iter->qnstep); - } -} - -static void iterupdate(Iteration *iter) -{ - size_t i, j, n, k; - - n = iter->n; - k = iter->k; - iter->f = iter->new_f; - iter->crit = iter->new_crit; - for (i = 0; i < n + k; i++) { - iter->x[i] = iter->new_x[i]; - iter->delf[i] = iter->new_delf[i]; - } - for (i = 0; i < k; i++) { - for (j = 0; j < n; j++) { - iter->hessf[n + i][j] = iter->new_delg[i][j]; - iter->hessf[j][n + i] = iter->new_delg[i][j]; - } - } -} - -static void -mindriver(Iteration *iter) -{ - iter->consecmax = 0; - iter->itncount = 0; - iter->termcode = 0; - if (! iter->values_supplied) { - eval_funval(iter); - if (iter->k == 0) eval_gradient(iter); - eval_hessian(iter); - } - - stoptest0(iter); - print_header(iter); - print_status(iter); - while (iter->termcode == 0) { - iter->itncount++; - print_header(iter); - findqnstep(iter); - linesearch(iter); - if (iter->k == 0) eval_next_gradient(iter); - stoptest(iter); - iterupdate(iter); - eval_hessian(iter); - print_status(iter); - } -} - -/************************************************************************/ -/** **/ -/** External Interface Routines **/ -/** **/ -/************************************************************************/ - -static Iteration myiter; - -size_t -minworkspacesize(size_t n, size_t k) -{ - size_t size; - - /* x, new_x, sx, delf, new_delf, qnstep and F */ - size = 7 * sizeof(double) * (n + k); - - /* hessf, H and L */ - size += 3 * (sizeof(double *) * (n + k) - + sizeof(double) * (n + k) * (n + k)); - - /* delg and new_delg */ - size += 2 * (sizeof(double *) * k + sizeof(double) * n * k); - - return(size); -} - -char * -minresultstring(int code) -{ - if (code <= 0) return("bad input data"); - else if (code <= 5) return(termcodes[code]); - else return("unknown return code"); -} - -void -minsetup(size_t n, size_t k, - void (*ffun)(), void (*gfun)(), - double *x, double typf, double *typx, char *work) - /* int (*ffun)(), (*gfun)();*/ -{ - Iteration *iter = &myiter; - size_t i, j; - double nx0, ntypx; - - n = (n > 0) ? n : 0; - k = (k > 0) ? k : 0; - - iter->n = n; - iter->k = k; - iter->ffun = ffun; /* FIXME:AJR */ - iter->gfun = gfun; - - iter->x = (double *) work; work += sizeof(double) * (n + k); - iter->new_x = (double *) work; work += sizeof(double) * (n + k); - iter->sx = (double *) work; work += sizeof(double) * (n + k); - iter->delf = (double *) work; work += sizeof(double) * (n + k); - iter->new_delf = (double *) work; work += sizeof(double) * (n + k); - iter->qnstep = (double *) work; work += sizeof(double) * (n + k); - iter->F = (double *) work; work += sizeof(double) * (n + k); - for (i = 0; i < n; i++) { - iter->x[i] = x[i]; - iter->sx[i] = (typx != nil && typx[i] > 0.0) ? 1.0 / typx[i] : 1.0; - } - for (i = 0; i < k; i++) { - iter->x[n + i] = x[n + i]; - iter->sx[n + i] = 1.0; - } - - iter->hessf = (double **) work; work += sizeof(double *) * (n + k); - for (i = 0; i < n + k; i++) { - iter->hessf[i] = (double *) work; - work += sizeof(double) * (n + k); - } - for (i = 0; i < n + k; i++) - for (j = 0; j < n + k; j++) iter->hessf[i][j] = 0.0; - iter->L = (double **) work; work += sizeof(double *) * (n + k); - for (i = 0; i < n + k; i++) { - iter->L[i] = (double *) work; - work += sizeof(double) * (n + k); - } - iter->H = (double **) work; work += sizeof(double *) * (n + k); - for (i = 0; i < n + k; i++) { - iter->H[i] = (double *) work; - work += sizeof(double) * (n + k); - } - - iter->new_delg = (k > 0) ? (double **) work : nil; - work += sizeof(double *) * k; - for (i = 0; i < k; i++) { - iter->new_delg[i] = (double *) work; - work += sizeof(double) * n; - } - - iter->typf = (typf > 0.0) ? typf : 1.0; - - iter->verbose = VERBOSE; - iter->use_line_search = USE_SEARCH; - iter->itncount = 0; - iter->itnlimit = ITNLIMIT; - iter->gradtol = pow(macheps(), GRADTOL_POWER); - iter->steptol = pow(macheps(), STEPTOL_POWER); - for (i = 0, nx0 = 0.0, ntypx = 0.0; i < iter->n; i++) { - nx0 += fabs(iter->new_x[i] / iter->sx[i]); - ntypx += fabs(1.0 / iter->sx[i]); - } - iter->maxstep = MAX_STEP_FACTOR * Max(nx0, ntypx); - iter->values_supplied = FALSE; -} - -void -minsetoptions(double gradtol, double steptol, double maxstep, - int itnlimit, int verbose, int use_search, int change_sign) -{ - Iteration *iter = &myiter; - - if (gradtol > 0.0) iter->gradtol = gradtol; - if (steptol > 0.0) iter->steptol = steptol; - if (maxstep > 0.0) iter->maxstep = maxstep; - if (itnlimit >= 0) iter->itnlimit = itnlimit; - if (verbose >= 0) iter->verbose = verbose; - iter->use_line_search = use_search; - iter->change_sign = change_sign; -} - -void -minsupplyvalues(double f, double *delf, double **hessf, - double *g, double **delg) -{ - Iteration *iter = &myiter; - size_t i, j, n, k; - - n = iter->n; - k = iter->k; - - iter->f = f; - for (i = 0; i < n; i++) { - iter->delf[i] = delf[i]; - for (j = 0; j < k; j++) - iter->delf[i] += iter->x[n + j] * delg[j][i]; - for (j = 0; j < n; j++) iter->hessf[i][j] = hessf[i][j]; - } - for (i = 0; i < k; i++) { - iter->delf[n + i] = g[i]; - for (j = 0; j < n; j++) { - iter->hessf[n + i][j] = delg[i][j]; - iter->hessf[j][n + i] = delg[i][j]; - } - } - if (k == 0) iter->crit = f; - else { - for (i = 0, iter->crit = 0.0; i < n + k; i++) - iter->crit += 0.5 * pow(fabs(iter->delf[i]), 2.0); - } - iter->values_supplied = TRUE; -} - -void -minimize() -{ - mindriver(&myiter); -} - -void -minresults(double *x, double *pf, double *pcrit, double *delf, double **hessf, - double *g, double **delg, int *pcount, int *ptermcode, double *diagadd) -{ - Iteration *iter = &myiter; - size_t i, j, n, k; - - n = iter->n; - k = iter->k; - - if (pf != nil) *pf = iter->f; - if (pcrit != nil) *pcrit = iter->crit; - for (i = 0; i < n; i++) { - if (x != nil) x[i] = iter->x[i]; - if (delf != nil) delf[i] = iter->delf[i]; - for (j = 0; j < n; j++) if (hessf != nil) hessf[i][j] = iter->hessf[i][j]; - } - for (i = 0; i < k; i++) { - if (x != nil) x[n + i] = iter->x[n + i]; - if (g != nil) g[i] = iter->delf[n + i]; - for (j = 0; j < n; j++) - if (delg != nil) delg[i][j] = iter->hessf[n + i][j]; - } - if (pcount != nil) *pcount = iter->itncount; - if (ptermcode != nil) *ptermcode = iter->termcode; - if (diagadd != nil) *diagadd = iter->diagadd; -} - - -double pdlogdet(int n, double **H) -{ - int i; - double logdet, maxadd; - - choldecomp(H, n, 0.0, &maxadd); - for (i = 0, logdet = 0.0; i < n; i++) logdet += 2.0 * log(H[i][i]); - return(logdet); -} - - -/* - TODO - - return amount added to make pos definite - - scaling for constraints - - alternate global strategies - - callback function for verbose mode -*/ - diff --git a/lib/nor.c b/lib/nor.c deleted file mode 100644 index 3e789cb..0000000 --- a/lib/nor.c +++ /dev/null @@ -1,100 +0,0 @@ -# include -# include "xmath.h" - -#define P10 242.66795523053175 -#define P11 21.979261618294152 -#define P12 6.9963834886191355 -#define P13 -.035609843701815385 -#define Q10 215.05887586986120 -#define Q11 91.164905404514901 -#define Q12 15.082797630407787 -#define Q13 1.0 - -#define P20 300.4592610201616005 -#define P21 451.9189537118729422 -#define P22 339.3208167343436870 -#define P23 152.9892850469404039 -#define P24 43.16222722205673530 -#define P25 7.211758250883093659 -#define P26 .5641955174789739711 -#define P27 -.0000001368648573827167067 -#define Q20 300.4592609569832933 -#define Q21 790.9509253278980272 -#define Q22 931.3540948506096211 -#define Q23 638.9802644656311665 -#define Q24 277.5854447439876434 -#define Q25 77.00015293522947295 -#define Q26 12.78272731962942351 -#define Q27 1.0 - -#define P30 -.00299610707703542174 -#define P31 -.0494730910623250734 -#define P32 -.226956593539686930 -#define P33 -.278661308609647788 -#define P34 -.0223192459734184686 -#define Q30 .0106209230528467918 -#define Q31 .191308926107829841 -#define Q32 1.05167510706793207 -#define Q33 1.98733201817135256 -#define Q34 1.0 - -#define SQRT2 1.414213562373095049 -#define SQRTPI 1.772453850905516027 - -void -normbase(double *x, double *phi) -{ - int sn; - double R1, R2, y, y2, y3, y4, y5, y6, y7, erf, erfc, z, z2, z3, z4; - - y = *x / SQRT2; - if (y < 0) { - y = -y; - sn = -1; - } else { - sn = 1; - } - y2 = y * y; - y4 = y2 * y2; - y6 = y4 * y2; - - if(y < 0.46875) { - R1 = P10 + P11 * y2 + P12 * y4 + P13 * y6; - R2 = Q10 + Q11 * y2 + Q12 * y4 + Q13 * y6; - erf = y * R1 / R2; - if(sn == 1) { - *phi = 0.5 + 0.5*erf; - } else { - *phi = 0.5 - 0.5*erf; - } - } else { - if(y < 4.0) { - y3 = y2 * y; - y5 = y4 * y; - y7 = y6 * y; - R1 = P20 + P21 * y + P22 * y2 + P23 * y3 + P24 * y4 + P25 * y5 + - P26 * y6 + P27 * y7; - R2 = Q20 + Q21 * y + Q22 * y2 + Q23 * y3 + Q24 * y4 + Q25 * y5 + - Q26 * y6 + Q27 * y7; - erfc = exp(-y2) * R1 / R2; - if (sn == 1) { - *phi = 1.0 - (0.5 * erfc); - } else { - *phi = 0.5 * erfc; - } - } else { - z = y4; - z2 = z * z; - z3 = z2 * z; - z4 = z2 * z2; - R1 = P30 + P31 * z + P32 * z2 + P33 * z3 + P34 * z4; - R2 = Q30 + Q31 * z + Q32 * z2 + Q33 * z3 + Q34 * z4; - erfc = (exp(-y2)/y) * (1.0 / SQRTPI + R1 / (R2 * y2)); - if(sn == 1) { - *phi = 1.0 - ( 0.5 * erfc); - } else { - *phi = (0.5 * erfc); - } - } - } -} diff --git a/lib/num_sfun.c b/lib/num_sfun.c deleted file mode 100644 index 6e328a3..0000000 --- a/lib/num_sfun.c +++ /dev/null @@ -1,655 +0,0 @@ -/* -(c) Copyright Taiichi Yuasa and Masami Hagiya, 1984. All rights reserved. -Copying of this file is authorized to users who have executed the true and -proper "License Agreement for Kyoto Common LISP" with SIGLISP. -*/ - -#include "include.h" -#include "num_include.h" - -object imag_unit, minus_imag_unit, imag_two; - -int -fixnum_expt(x, y) -{ - int z; - - z = 1; - while (y > 0) - if (y%2 == 0) { - x *= x; - y /= 2; - } else { - z *= x; - --y; - } - return(z); -} - -object -number_exp(x) -object x; -{ - double exp(); - - switch (type_of(x)) { - - case t_fixnum: - case t_bignum: - case t_ratio: - return(make_longfloat(exp(number_to_double(x)))); - - case t_shortfloat: - return(make_shortfloat((shortfloat)exp((double)(sf(x))))); - - case t_longfloat: - return(make_longfloat(exp(lf(x)))); - - case t_complex: - { - object y, y1; - object number_sin(), number_cos(); - vs_mark; - - y = x->cmp.cmp_imag; - x = x->cmp.cmp_real; - x = number_exp(x); - vs_push(x); - y1 = number_cos(y); - vs_push(y1); - y = number_sin(y); - vs_push(y); - y = make_complex(y1, y); - vs_push(y); - x = number_times(x, y); - vs_reset; - return(x); - } - - default: - FEwrong_type_argument(Snumber, x); - } -} - -object -number_expt(x, y) -object x, y; -{ - enum type tx, ty; - object z, number_nlog(); - vs_mark; - - tx = type_of(x); - ty = type_of(y); - if (ty == t_fixnum && fix(y) == 0) - switch (tx) { - case t_fixnum: case t_bignum: case t_ratio: - return(small_fixnum(1)); - - case t_shortfloat: - return(make_shortfloat((shortfloat)1.0)); - - case t_longfloat: - return(make_longfloat(1.0)); - - case t_complex: - z = number_expt(x->cmp.cmp_real, y); - vs_push(z); - z = make_complex(z, small_fixnum(0)); - vs_reset; - return(z); - - default: - FEwrong_type_argument(Snumber, x); - } - if (number_zerop(x)) { - if (!number_plusp(ty==t_complex?y->cmp.cmp_real:y)) - FEerror("Cannot raise zero to the power ~S.", 1, y); - return(number_times(x, y)); - } - if (ty == t_fixnum || ty == t_bignum) { - if (number_minusp(y)) { - z = number_negate(y); - vs_push(z); - z = number_expt(x, z); - vs_push(z); - z = number_divide(small_fixnum(1), z); - vs_reset; - return(z); - } - z = small_fixnum(1); - vs_push(z); - vs_push(Cnil); - vs_push(Cnil); - while (number_plusp(y)) - if (number_evenp(y)) { - x = number_times(x, x); - vs_top[-1] = x; - y = integer_divide1(y, small_fixnum(2)); - vs_top[-2] = y; - } else { - z = number_times(z, x); - vs_top[-3] = z; - y = number_minus(y, small_fixnum(1)); - vs_top[-2] = y; - } - vs_reset; - return(z); - } - z = number_nlog(x); - vs_push(z); - z = number_times(z, y); - vs_push(z); - z = number_exp(z); - vs_reset; - return(z); -} - -object -number_nlog(x) -object x; -{ - double log(); - object r, i, a, p, number_sqrt(), number_atan2(); - vs_mark; - - if (type_of(x) == t_complex) { - r = x->cmp.cmp_real; - i = x->cmp.cmp_imag; - goto COMPLEX; - } - if (number_zerop(x)) - FEerror("Zero is the logarithmic singularity.", 0); - if (number_minusp(x)) { - r = x; - i = small_fixnum(0); - goto COMPLEX; - } - switch (type_of(x)) { - case t_fixnum: - case t_bignum: - case t_ratio: - return(make_longfloat(log(number_to_double(x)))); - - case t_shortfloat: - return(make_shortfloat((shortfloat)log((double)(sf(x))))); - - case t_longfloat: - return(make_longfloat(log(lf(x)))); - - default: - FEwrong_type_argument(Snumber, x); - } - -COMPLEX: - a = number_times(r, r); - vs_push(a); - p = number_times(i, i); - vs_push(p); - a = number_plus(a, p); - vs_push(a); - a = number_nlog(a); - vs_push(a); - a = number_divide(a, small_fixnum(2)); - vs_push(a); - p = number_atan2(i, r); - vs_push(p); - x = make_complex(a, p); - vs_reset; - return(x); -} - -object -number_log(x, y) -object x, y; -{ - object z; - vs_mark; - - if (number_zerop(y)) - FEerror("Zero is the logarithmic singularity.", 0); - if (number_zerop(x)) - return(number_times(x, y)); - x = number_nlog(x); - vs_push(x); - y = number_nlog(y); - vs_push(y); - z = number_divide(y, x); - vs_reset; - return(z); -} - -object -number_sqrt(x) -object x; -{ - object z; - double sqrt(); - vs_mark; - - if (type_of(x) == t_complex) - goto COMPLEX; - if (number_minusp(x)) - goto COMPLEX; - switch (type_of(x)) { - case t_fixnum: - case t_bignum: - case t_ratio: - return(make_longfloat(sqrt(number_to_double(x)))); - - case t_shortfloat: - return(make_shortfloat((shortfloat)sqrt((double)(sf(x))))); - - case t_longfloat: - return(make_longfloat(sqrt(lf(x)))); - - default: - FEwrong_type_argument(Snumber, x); - } - -COMPLEX: - z = make_ratio(small_fixnum(1), small_fixnum(2)); - vs_push(z); - z = number_expt(x, z); - vs_reset; - return(z); -} - -object -number_atan2(y, x) -object y, x; -{ - object z; - double atan(), dy, dx, dz; - - dy = number_to_double(y); - dx = number_to_double(x); - if (dx > 0.0) - if (dy > 0.0) - dz = atan(dy / dx); - else if (dy == 0.0) - dz = 0.0; - else - dz = -atan(-dy / dx); - else if (dx == 0.0) - if (dy > 0.0) - dz = PI / 2.0; - else if (dy == 0.0) - FEerror("Logarithmic singularity.", 0); - else - dz = -PI / 2.0; - else - if (dy > 0.0) - dz = PI - atan(dy / -dx); - else if (dy == 0.0) - dz = PI; - else - dz = -PI + atan(-dy / -dx); - z = make_longfloat(dz); - return(z); -} - -object -number_atan(y) -object y; -{ - object z, z1; - vs_mark; - - if (type_of(y) == t_complex) { - z = number_times(imag_unit, y); - vs_push(z); - z = one_plus(z); - vs_push(z); - z1 = number_times(y, y); - vs_push(z1); - z1 = one_plus(z1); - vs_push(z1); - z1 = number_sqrt(z1); - vs_push(z1); - z = number_divide(z, z1); - vs_push(z); - z = number_nlog(z); - vs_push(z); - z = number_times(minus_imag_unit, z); - vs_reset; - return(z); - } - return(number_atan2(y, small_fixnum(1))); -} - -object -number_sin(x) -object x; -{ - double sin(); - - switch (type_of(x)) { - - case t_fixnum: - case t_bignum: - case t_ratio: - return(make_longfloat(sin(number_to_double(x)))); - - case t_shortfloat: - return(make_shortfloat((shortfloat)sin((double)(sf(x))))); - - case t_longfloat: - return(make_longfloat(sin(lf(x)))); - - case t_complex: - { - object r; - object x0, x1, x2; - vs_mark; - - x0 = number_times(imag_unit, x); - vs_push(x0); - x0 = number_exp(x0); - vs_push(x0); - x1 = number_times(minus_imag_unit, x); - vs_push(x1); - x1 = number_exp(x1); - vs_push(x1); - x2 = number_minus(x0, x1); - vs_push(x2); - r = number_divide(x2, imag_two); - - vs_reset; - return(r); - } - - default: - FEwrong_type_argument(Snumber, x); - - } -} - -object -number_cos(x) -object x; -{ - double cos(); - - switch (type_of(x)) { - - case t_fixnum: - case t_bignum: - case t_ratio: - return(make_longfloat(cos(number_to_double(x)))); - - case t_shortfloat: - return(make_shortfloat((shortfloat)cos((double)(sf(x))))); - - case t_longfloat: - return(make_longfloat(cos(lf(x)))); - - case t_complex: - { - object r; - object x0, x1, x2; - vs_mark; - - x0 = number_times(imag_unit, x); - vs_push(x0); - x0 = number_exp(x0); - vs_push(x0); - x1 = number_times(minus_imag_unit, x); - vs_push(x1); - x1 = number_exp(x1); - vs_push(x1); - x2 = number_plus(x0, x1); - vs_push(x2); - r = number_divide(x2, small_fixnum(2)); - - vs_reset; - return(r); - } - - default: - FEwrong_type_argument(Snumber, x); - - } -} - -object -number_tan(x) -object x; -{ - object r, s, c; - vs_mark; - - s = number_sin(x); - vs_push(s); - c = number_cos(x); - vs_push(c); - if (number_zerop(c) == TRUE) - FEerror("Cannot compute the tangent of ~S.", 1, x); - r = number_divide(s, c); - vs_reset; - return(r); -} - -object -number_asin(x) - object x; -{ - double asin(); - double dx; - - /* check for a real argument in [-1,1] */ - if (type_of(x) != t_complex) { - switch (type_of(x)) { - case t_fixnum: - case t_bignum: - case t_ratio: - dx = number_to_double(x); - break; - case t_shortfloat: - dx = sf(x); - break; - case t_longfloat: - dx = lf(x); - break; - default: - FEwrong_type_argument(Snumber, x); - } - if (-1.0 <= dx && dx <= 1.0) return(make_longfloat(asin(dx))); - } - - /* treat as complex argument, result */ - { - object r; - object x0, x1; - vs_mark; - - x0 = number_times(x, x); - vs_push(x0); - x0 = number_minus(small_fixnum(1), x0); - vs_push(x0); - x0 = number_sqrt(x0); - vs_push(x0); - x1 = number_times(imag_unit, x); - vs_push(x1); - x0 = number_plus(x0, x1); - vs_push(x0); - x0 = number_nlog(x0); - vs_push(x0); - r = number_times(minus_imag_unit, x0); - vs_reset; - return(r); - } -} - -object -number_acos(x) - object x; -{ - double acos(); - double dx; - - /* check for a real argument in [-1,1] */ - if (type_of(x) != t_complex) { - switch (type_of(x)) { - case t_fixnum: - case t_bignum: - case t_ratio: - dx = number_to_double(x); - break; - case t_shortfloat: - dx = sf(x); - break; - case t_longfloat: - dx = lf(x); - break; - default: - FEwrong_type_argument(Snumber, x); - } - if (-1.0 <= dx && dx <= 1.0) return(make_longfloat(acos(dx))); - } - - /* treat as complex argument, result */ - { - object r; - object x0; - vs_mark; - - x0 = number_times(x, x); - vs_push(x0); - x0 = number_minus(small_fixnum(1), x0); - vs_push(x0); - x0 = number_sqrt(x0); - vs_push(x0); - x0 = number_times(imag_unit, x0); - vs_push(x0); - x0 = number_plus(x0, x); - vs_push(x0); - x0 = number_nlog(x0); - vs_push(x0); - r = number_times(minus_imag_unit, x0); - vs_reset; - return(r); - } -} - -Lexp() -{ - check_arg(1); - check_type_number(&vs_base[0]); - vs_base[0] = number_exp(vs_base[0]); -} - -Lexpt() -{ - check_arg(2); - check_type_number(&vs_base[0]); - check_type_number(&vs_base[1]); - vs_base[0] = number_expt(vs_base[0], vs_base[1]); - vs_pop; -} - -Llog() -{ - int narg; - - narg = vs_top - vs_base; - if (narg < 1) - too_few_arguments(); - else if (narg == 1) { - check_type_number(&vs_base[0]); - vs_base[0] = number_nlog(vs_base[0]); - } else if (narg == 2) { - check_type_number(&vs_base[0]); - check_type_number(&vs_base[1]); - vs_base[0] = number_log(vs_base[1], vs_base[0]); - vs_pop; - } else - too_many_arguments(); -} - -Lsqrt() -{ - check_arg(1); - check_type_number(&vs_base[0]); - vs_base[0] = number_sqrt(vs_base[0]); -} - -Lsin() -{ - check_arg(1); - check_type_number(&vs_base[0]); - vs_base[0] = number_sin(vs_base[0]); -} - -Lcos() -{ - check_arg(1); - check_type_number(&vs_base[0]); - vs_base[0] = number_cos(vs_base[0]); -} - -Ltan() -{ - check_arg(1); - check_type_number(&vs_base[0]); - vs_base[0] = number_tan(vs_base[0]); -} - -Latan() -{ - int narg; - - narg = vs_top - vs_base; - if (narg < 1) - too_few_arguments(); - if (narg == 1) { - check_type_number(&vs_base[0]); - vs_base[0] = number_atan(vs_base[0]); - } else if (narg == 2) { - check_type_or_rational_float(&vs_base[0]); - check_type_or_rational_float(&vs_base[1]); - vs_base[0] = number_atan2(vs_base[0], vs_base[1]); - vs_pop; - } else - too_many_arguments(); -} - -Lasin() -{ - check_arg(1); - check_type_number(&vs_base[0]); - vs_base[0] = number_asin(vs_base[0]); -} - -Lacos() -{ - check_arg(1); - check_type_number(&vs_base[0]); - vs_base[0] = number_acos(vs_base[0]); -} - -init_num_sfun() -{ - imag_unit - = make_complex(make_longfloat(0.0), make_longfloat(1.0)); - enter_mark_origin(&imag_unit); - minus_imag_unit - = make_complex(make_longfloat(0.0), make_longfloat(-1.0)); - enter_mark_origin(&minus_imag_unit); - imag_two - = make_complex(make_longfloat(0.0), make_longfloat(2.0)); - enter_mark_origin(&imag_two); - - make_constant("PI", make_longfloat(PI)); - - make_function("EXP", Lexp); - make_function("EXPT", Lexpt); - make_function("LOG", Llog); - make_function("SQRT", Lsqrt); - make_function("SIN", Lsin); - make_function("COS", Lcos); - make_function("TAN", Ltan); - make_function("ATAN", Latan); - make_function("ASIN", Lasin); - make_function("ACOS", Lacos); -} diff --git a/lib/ppnd.c b/lib/ppnd.c deleted file mode 100644 index f8f218a..0000000 --- a/lib/ppnd.c +++ /dev/null @@ -1,60 +0,0 @@ -#include "xmath.h" - -#define zero 0.0 -#define half 0.5 -#define one 1.0 -#define split 0.42e0 -#define a0 2.50662823884e0 -#define a1 -18.61500062529e0 -#define a2 41.39119773534e0 -#define a3 -25.44106049637e0 -#define b1 -8.47351093090e0 -#define b2 23.08336743743e0 -#define b3 -21.06224101826e0 -#define b4 3.13082909833e0 - -#define c0 -2.78718931138e0 -#define c1 -2.29796479134e0 -#define c2 4.85014127135e0 -#define c3 2.32121276858e0 -#define d1 3.54388924762e0 -#define d2 1.63706781897e0 - -/* -c -c Algorithm as 111 Applied statistics (1977), vol 26 no 1 page 121 -c Produces normal deviate corresponding to lower tail area of p -c the hash sums are the sums of the moduli of the coefficients -c they nave no inherent meanings but are incuded for use in -c checking transcriptions. Functions abs,alog and sqrt are used. -c - -derived from AS111 fortran version -*/ - -double ppnd(p, ifault) - double p; - int *ifault; -{ - double q,r,ppn; - - *ifault = 0; - q = p - half; - if( fabs(q) <= split) { - r = q*q; - ppn = q * (((a3 * r + a2) * r + a1) * r + a0) - / ((((b4 * r + b3) * r + b2) * r + b1) * r + one); - } - else { - r = p; - if(q > zero) r = one - p; - if(r <= zero) { - *ifault = 1; - return(0.0); - } - r = sqrt(-log(r)); - ppn = (((c3*r+c2)*r + c1) * r + c0) / ((d2 * r + d1) * r + one); - if( q < zero) ppn = -ppn; - } - return(ppn); -} diff --git a/lib/qrdecomp.c b/lib/qrdecomp.c deleted file mode 100644 index 8e04db8..0000000 --- a/lib/qrdecomp.c +++ /dev/null @@ -1,186 +0,0 @@ -/* adapted from DQRDC of LINPACK */ - -#include "linalg.h" - -#define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) - -static double -NORM2(int i, int j, int n, double **x) -{ - int k; - double maxx, sum, temp; - - for (k = i, maxx = 0.0; k < n; k++) { - temp = fabs(x[k][j]); - if (maxx < temp) maxx = temp; - } - if (maxx == 0.0) return(0.0); - else { - for (k = i, sum = 0.0; k < n; k++) { - temp = x[k][j] / maxx; - sum += temp * temp; - } - return(maxx * sqrt(sum)); - } -} - -static double -DOT(int i, int j, int k, int n, double **x) -{ - int l; - double sum; - - for (l = i, sum = 0.0; l < n; l++) sum += x[l][j] * x[l][k]; - return(sum); -} - -static void -AXPY(int i, int j, int k, int n, - double a, double **x) -{ - int l; - for (l = i; l < n; l++) { - x[l][k] = a * x[l][j] + x[l][k]; - } -} - -static void -SCALE(int i, int j, int n, double a, double **x) -{ - int k; - for (k = i; k < n; k++) { - x[k][j] *= a; - } -} - -static void -SWAP(int i, int j, int n, double **a) -{ - int k; - double temp; - for (k = 0; k < n; k++) { - temp = a[k][i]; - a[k][i] = a[k][j]; - a[k][j] = temp; - } -} - -void -qrdecomp(double **x, int n, int p, - double **v, int *jpvt, int pivot) -{ - int i,j,k,jp,l,lp1,lup,maxj; - double maxnrm,tt,*qraux,*work; - double nrmxl,t; - - if (n < 0) { - return; - } - work = v[0]; - qraux = rvector(p); - - /* - * compute the norms of the free columns. - */ - if (pivot) - for (j = 0; j < p; j++) { - jpvt[j] = j; - qraux[j] = NORM2(0, j, n, x); - work[j] = qraux[j]; - } - /* - * perform the householder reduction of x. - */ - lup = (n < p) ? n : p; - for (l = 0; l < lup; l++) { - if (pivot) { - /* - * locate the column of largest norm and bring it - * into the pivot position. - */ - maxnrm = 0.0; - maxj = l; - for (j = l; j < p; j++) - if (qraux[j]>maxnrm) { - maxnrm = qraux[j]; - maxj = j; - } - if (maxj!=l) { - SWAP(l, maxj, n, x); - qraux[maxj] = qraux[l]; - work[maxj] = work[l]; - jp = jpvt[maxj]; - jpvt[maxj] = jpvt[l]; - jpvt[l] = jp; - } - } - qraux[l] = 0.0; - if (l != n-1) { - /* - * compute the householder transformation for column l. - */ - nrmxl = NORM2(l, l, n, x); - if (nrmxl != 0.0) { - if (x[l][l] != 0.0) - nrmxl = SIGN(nrmxl, x[l][l]); - SCALE(l, l, n, 1.0 / nrmxl, x); - x[l][l] = 1.0+x[l][l]; - /* - * apply the transformation to the remaining columns, - * updating the norms. - */ - lp1 = l+1; - for (j = lp1; j < p; j++) { - t = -DOT(l, l, j, n, x) / x[l][l]; - AXPY(l, l, j, n, t, x); - if (pivot) - if (qraux[j]!=0.0) { - tt = 1.0-(fabs(x[l][j])/qraux[j])*(fabs(x[l][j])/qraux[j]); - if (tt < 0.0) tt = 0.0; - t = tt; - tt = 1.0+0.05*tt*(qraux[j]/work[j])*(qraux[j]/work[j]); - if (tt!=1.0) - qraux[j] = qraux[j]*sqrt(t); - else { - qraux[j] = NORM2(l+1, j, n, x); - work[j] = qraux[j]; - } - } - } - /* - * save the transformation. - */ - qraux[l] = x[l][l]; - x[l][l] = -nrmxl; - } - } - } - - /* copy over the upper triangle of a */ - for (i = 0; i < p; i++) { - for (j = 0; j < i; j++) v[i][j] = 0.0; - for (j = i; j < p; j++) { - v[i][j] = x[i][j]; - } - } - - /* accumulate the Q transformation -- assumes p <= n */ - for (i = 0; i < p; i++) { - x[i][i] = qraux[i]; - for (k = 0; k < i; k++) x[k][i] = 0.0; - } - for (i = p - 1; i >= 0; i--) { - if (i == n - 1) x[i][i] = 1.0; - else { - for (k = i; k < n; k++) x[k][i] = -x[k][i]; - x[i][i] += 1.0; - } - for (j = i - 1; j >= 0; j--) { - if (x[j][j] != 0.0) { - t = -DOT(j, j, i, n, x) / x[j][j]; - AXPY(j, j, i, n, t, x); - } - } - } - free_vector(qraux); -} diff --git a/lib/rcondest.c b/lib/rcondest.c deleted file mode 100644 index 03c5e17..0000000 --- a/lib/rcondest.c +++ /dev/null @@ -1,102 +0,0 @@ -/* rcondest - Estimates reciprocal of condition number. */ -/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */ -/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */ -/* You may give out copies of this software; for conditions see the */ -/* file COPYING included with this distribution. */ - -#include "linalg.h" - -static -double Max(double a, double b) -{ - return(a > b ? a : b); -} - - -void -backsolve(double **a, double *b, size_t n, int mode) -{ - long i, j; /* Must be signed, since counting down in loops... */ - RMatrix ra = (RMatrix) a; - RVector rb = (RVector) b; - CMatrix ca = (CMatrix) a; - CVector cb = (CVector) b; - - switch (mode) { - case RE: - for (i = n - 1; i >= 0; i--) { - if (ra[i][i] != 0.0) rb[i] = rb[i] / ra[i][i]; - for (j = i + 1; j < n; j++) rb[i] -= ra[i][j] * rb[j]; - } - break; - case CX: - for (i = n - 1; i >= 0; i--) { - if (modulus(ca[i][i]) != 0.0) cb[i] = cdiv(cb[i], ca[i][i]); - for (j = i + 1; j < n; j++) - cb[i] = csub(cb[i], cmul(ca[i][j], cb[j])); - } - break; - } -} - -/** Exported function **/ - -double rcondest(double **a, size_t n) -{ - RVector p, pm, x; - double est, xp, xm, temp, tempm, xnorm; - long i, j; - - for (i = 0; i < n; i++) - if (a[i][i] == 0.0) return(0.0); - - p = rvector(n); - pm = rvector(n); - x = rvector(n); - - /* Set est to reciprocal of L1 matrix norm of A */ - est = fabs(a[0][0]); - - for (j = 1; j < n; j++) { - for (i = 0, temp = fabs(a[j][j]); i < j; i++) - temp += fabs(a[i][j]); - est = Max(est, temp); - } - est = 1.0 / est; - - /* Solve A^Tx = e, selecting e as you proceed */ - x[0] = 1.0 / a[0][0]; - for (i = 1; i < n; i++) p[i] = a[0][i] * x[0]; - for (j = 1; j < n; j++) { - /* select ej and calculate x[j] */ - xp = ( 1.0 - p[j]) / a[j][j]; - xm = (-1.0 - p[j]) / a[j][j]; - temp = fabs(xp); - tempm = fabs(xm); - for (i = j + 1; i < n; i++) { - pm[i] = p[i] + a[j][i] * xm; - tempm += fabs(pm[i] / a[i][i]); - p[i] += a[j][i] * xp; - temp += fabs(p[i] / a[i][i]); - } - if (temp >= tempm) x[j] = xp; - else { - x[j] = xm; - for (i = j + 1; i < n; i++) p[i] = pm[i]; - } - } - - for (j = 0, xnorm = 0.0; j < n; j++) xnorm += fabs(x[j]); - est = est * xnorm; - - backsolve(a, x, n, RE); - for (j = 0, xnorm = 0.0; j < n; j++) xnorm += fabs(x[j]); - if (xnorm > 0) est = est / xnorm; - - free_vector(p); - free_vector(pm); - free_vector(x); - - return(est); -} - diff --git a/lib/splines.c b/lib/splines.c deleted file mode 100644 index 585d2c7..0000000 --- a/lib/splines.c +++ /dev/null @@ -1,112 +0,0 @@ -#include "xmath.h" - -/* -#ifdef __APPLE__ -#include -#endif -*/ - -#include - -/* natural cubic spline interpolation based on Numerical Recipes in C */ - -/* calculate second derivatives; assumes strictly increasing x values */ -static void -find_spline_derivs(double *x, double *y, size_t n, - double *y2, double *u) -{ - long i, k; - double p, sig; - - y2[0] = u[0] = 0.0; /* lower boundary condition for natural spline */ - - /* decomposition loop for the tridiagonal algorithm */ - for (i = 1; i < n - 1; i++) { - y2[i] = u[i] = 0.0; /* set in case a zero test is failed */ - if (x[i - 1] < x[i] && x[i] < x[i + 1]) { - sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]); - p = sig * y2[i - 1] + 2.0; - if (p != 0.0) { - y2[i] = (sig - 1.0) / p; - u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - - (y[i] - y[i - 1]) / (x[i] - x[i - 1]); - u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p; - } - } - } - - /* upper boundary condition for natural spline */ - y2[n - 1] = 0.0; - - /* backsubstitution loop of the tridiagonal algorithm */ - for (k = n - 2; k >= 0; k--) - y2[k] = y2[k] * y2[k + 1] + u[k]; -} - -/* interpolate or extrapolate value at x using results of find_spline_derivs */ -static void -spline_interp(double *xa, double *ya, double *y2a, - size_t n, double x, double *y) -{ - long klo, khi, k; - double h, b, a; - - if (x <= xa[0]) { - h = xa[1] - xa[0]; - b = (h > 0.0) ? (ya[1] - ya[0]) / h - h * y2a[1] / 6.0 : 0.0; - *y = ya[0] + b * (x - xa[0]); - } - else if (x >= xa[n - 1]) { - h = xa[n - 1] - xa[n - 2]; - b = (h > 0.0) ? (ya[n - 1] - ya[n - 2]) / h + h * y2a[n - 2] / 6.0 : 0.0; - *y = ya[n - 1] + b * (x - xa[n - 1]); - } - else { - /* try a linear interpolation for equally spaced x values */ - k = (n - 1) * (x - xa[0]) / (xa[n - 1] - xa[0]); - - /* make sure the range is right */ - if (k < 0) k = 0; - if (k > n - 2) k = n - 2; - - /* bisect if necessary until the bracketing interval is found */ - klo = (x >= xa[k]) ? k : 0; - khi = (x < xa[k + 1]) ? k + 1 : n - 1; - while (khi - klo > 1) { - k = (khi + klo) / 2; - if (xa[k] > x) khi = k; - else klo = k; - } - - /* interpolate */ - h = xa[khi] - xa[klo]; - if (h > 0.0) { - a = (xa[khi] - x) / h; - b = (x - xa[klo]) / h; - *y = a * ya[klo] + b * ya[khi] - + ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0; - } - else *y = (ya[klo] + ya[khi]) / 2.0; /* should not be needed */ - } -} - -int -fit_spline(size_t n, double *x, double *y, - size_t ns, double *xs, double *ys, double *work) -{ - size_t i; - double *y2, *u; - - y2 = work; u = work + n; - - if (n < 2 || ns < 1) return (1); /* signal an error */ - for (i = 1; i < n; i++) - if (x[i - 1] >= x[i]) return(1); /* signal an error */ - - find_spline_derivs(x, y, n, y2, u); - - for (i = 0; i < ns; i++) - spline_interp(x, y, y2, n, xs[i], &ys[i]); - - return(0); -} diff --git a/lib/studentbase.c b/lib/studentbase.c deleted file mode 100644 index 4a6d97b..0000000 --- a/lib/studentbase.c +++ /dev/null @@ -1,125 +0,0 @@ -#include "xmath.h" - -#define TWOVRPI 0.636619772367581343 -#define HALF_PI 1.5707963268 -#define EPSILON .000001 - -extern void betabase(double *, double *, double *,long *, long *,double *); -extern void normbase(double *,double *); - -extern double ppnd(), ppbeta(); - -/* CACM Algorithm 395, by G. W. Hill */ - -void -studentbase(double *x, double *df, double *cdf) -{ - double t, y, b, a, z, j, n; - - n = *df; - z = 1.0; - t = (*x) * (*x); - y = t / n; - b = 1.0 + y; - - if (n > floor(n) || (n >= 20 && t < n) || (n > 20)) { - if (n < 2.0 && n != 1.0) { - /* beta integral aproximation for small df */ - double da = 0.5, db = 0.5 * n, dx, dp; - long ia = 0, ib = floor(db); - - dx = db / (db + da * t); - betabase(&dx, &db, &da, &ib, &ia, &dp); - *cdf = (*x >= 0) ? 1.0 - .5 * dp : .5 * dp; - } else { - /* asymptotic series for large or non-integer df */ - if(y > EPSILON) { - y = log(b); - } - a = n - 0.5; - b = 48.0 * a * a; - y = a * y; - y = (((((-0.4 * y - 3.3) * y - 24.0) * y - 85.5) - / (0.8 * y * y + 100.0 + b) + y + 3.0) / b + 1.0) * sqrt(y); - - y = -1.0 * y; - normbase(&y, cdf); - if (*x > 0.0) *cdf = 1.0 - *cdf; - } - } else { - /* nested summation of cosine series */ - if (n < 20 && t < 4.0) { - a = y = sqrt(y); - if(n == 1.0) a = 0.0; - } else { - a = sqrt(b); - y = a * n; - for(j = 2; fabs(a - z) > EPSILON; j += 2.0) { - z = a; - y = (y * (j - 1)) / (b * j); - a = a + y / (n + j); - } - n += 2.0; - z = y = 0.0; - a = -a; - } - for(n = n - 2.0; n > 1.0; n -= 2.0) a = ((n - 1.0) / (b * n)) * a + y; - a = (fabs(n) < EPSILON) ? a/sqrt(b) : TWOVRPI * (atan(y) + a / b); - *cdf = z - a; - if(*x > 0.0) { - *cdf = 1.0 - 0.5 * *cdf; - } else { - *cdf = 0.5 * *cdf; - } - } -} - -/* CACM Algorithm 396, by G. W. Hill */ - -double ppstudent(pp, n, ifault) - double pp, n; - int *ifault; -{ - double sq, p, a, b, c, d, x, y; - - /* convert to double upper tailed probability */ - p = (pp < 0.5) ? 2.0 * pp : 2.0 * (1.0 - pp); - - if (n <= 3.0) { - if (n == 1) sq = tan(HALF_PI * (1.0 - p)); - else if (n == 2.0) sq = sqrt(2.0 / (p * (2.0 - p)) - 2.0); - else { - sq = ppbeta(p, 0.5 * n, 0.5, ifault); - if (sq != 0.0) sq = sqrt((n / sq) - n); - } - } - else { - a = 1.0 / (n - 0.5); - b = 48.0 / (a * a); - c = ((20700.0 * a / b - 98.0) * a - 16) * a + 96.36; - d = ((94.5 / (b + c) - 3.0) / b + 1.0) * sqrt(a * HALF_PI) * n; - x = d * p; - y = pow(x, 2.0 / n); - if (y > 0.05 + a) { - /* asymptotic inverse expansion about normal */ - x = ppnd(0.5 * p, ifault); - y = x * x; - if (n < 5) c = c + 0.3 * (n - 4.5) * (x + 0.6); - c = (((0.05 * d * x - 5.0) * x - 7.0) * x - 2.0) * x + b + c; - y = (((((0.4 * y + 6.3) * y + 36.0) * y + 94.5) / c - y - 3.0) / b - + 1.0) * x; - y = a * y * y; - y = (y > .002) ? exp(y) - 1.0 : 0.5 * y * y + y; - } - else { - y = ((1.0 / (((n + 6.0) / (n * y) - 0.089 * d - 0.822) - * (n + 2.0) * 3.0) + 0.5 / (n + 4.0)) * y - 1.0) - * (n + 1.0) / (n + 2.0) + 1.0 / y; - } - sq = sqrt(n * y); - } - - /* decode based on tail */ - if (pp < 0.5) sq = -sq; - return(sq); -} diff --git a/lib/svdecomp.c b/lib/svdecomp.c deleted file mode 100644 index 5be0073..0000000 --- a/lib/svdecomp.c +++ /dev/null @@ -1,263 +0,0 @@ -/* svdecomp - SVD decomposition routines. */ -/* Taken from Numerical Recipies. */ -/* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney */ -/* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz */ -/* You may give out copies of this software; for conditions see the */ -/* file COPYING included with this distribution. */ - -#include "linalg.h" - -#include - -static double PYTHAG(double a, double b) -{ - double at = fabs(a), bt = fabs(b), ct, result; - - if (at > bt) { ct = bt / at; result = at * sqrt(1.0 + ct * ct); } - else if (bt > 0.0) { ct = at / bt; result = bt * sqrt(1.0 + ct * ct); } - else result = 0.0; - return(result); -} - -#define SWAPD(a, b) (temp = (a), (a) = (b), (b) = temp) - -static void -sort_sv(size_t m, size_t n, size_t k, - double **a, double *w, double **v) -{ - size_t i, j; - double temp; - - for (i = k; (i < n - 1) && (w[i] < w[i+1]); i++) { - SWAPD(w[i], w[i+1]); - for (j = 0; j < m; j++) SWAPD(a[j][i], a[j][i+1]); - for (j = 0; j < n; j++) SWAPD(v[j][i], v[j][i+1]); - } -} - -static double maxarg1, maxarg2; -#define Max(a, b) (maxarg1 = (a), maxarg2 = (b), (maxarg1) > (maxarg2) ? (maxarg1) : (maxarg2)) -#define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) - -int -svdcmp(double **a, size_t m, size_t n, double *w, double **v) -{ - int flag, its; - long i, k, l; /* These must be signed */ - size_t j, jj, nm; - double c, f, h, s, x, y, z; - double anorm = 0.0, g = 0.0, scale = 0.0; - RVector rv1; - - if (m < n) return(FALSE); /* flag an error if m < n */ - - rv1 = rvector(n); - - /* Householder reduction to bidiagonal form */ - for (i = 0; i < n; i++) { - - /* left-hand reduction */ - l = i + 1; - rv1[i] = scale * g; - g = s = scale = 0.0; - if (i < m) { - for (k = i; k < m; k++) scale += fabs(a[k][i]); - if (scale) { - for (k = i; k < m; k++) { - a[k][i] /= scale; - s += a[k][i] * a[k][i]; - } - f = a[i][i]; - g = -SIGN(sqrt(s), f); - h = f * g - s; - a[i][i] = f - g; - if (i != n - 1) { - for (j = l; j < n; j++) { - for (s = 0.0, k = i; k < m; k++) s += a[k][i] * a[k][j]; - f = s / h; - for (k = i; k < m; k++) a[k][j] += f * a[k][i]; - } - } - for (k = i; k < m; k++) a[k][i] *= scale; - } - } - w[i] = scale * g; - - /* right-hand reduction */ - g = s = scale = 0.0; - if (i < m && i != n - 1) { - for (k = l; k < n; k++) scale += fabs(a[i][k]); - if (scale) { - for (k = l; k < n; k++) { - a[i][k] /= scale; - s += a[i][k] * a[i][k]; - } - f = a[i][l]; - g = -SIGN(sqrt(s), f); - h = f * g - s; - a[i][l] = f - g; - for (k = l; k < n; k++) rv1[k] = a[i][k] / h; - if (i != m - 1) { - for (j = l; j < m; j++) { - for (s = 0.0, k = l; k < n; k++) s += a[j][k] * a[i][k]; - for (k = l; k < n; k++) a[j][k] += s * rv1[k]; - } - } - for (k = l; k < n; k++) a[i][k] *= scale; - } - } - anorm = Max(anorm, (fabs(w[i]) + fabs(rv1[i]))); - } - - /* accumulate the right-hand transformation */ - for (i = n - 1; i >= 0; i--) { - if (i < n - 1) { - if (g) { - for (j = l; j < n; j++) - v[j][i] = (a[i][j] / a[i][l]) / g; - - for (j = l; j < n; j++) { - for (s = 0.0, k = l; k < n; k++) s += a[i][k] * v[k][j]; - for (k = l; k < n; k++) v[k][j] += s * v[k][i]; - } - } - for (j = l; j < n; j++) v[i][j] = v[j][i] = 0.0; - } - v[i][i] = 1.0; - g = rv1[i]; - l = i; - } - - /* accumulate the left-hand transformation */ - for (i = n - 1; i >= 0; i--) { - l = i + 1; - g = w[i]; - if (i < n - 1) - for (j = l; j < n; j++) a[i][j] = 0.0; - if (g) { - g = 1.0 / g; - if (i != n - 1) { - for (j = l; j < n; j++) { - for (s = 0.0, k = l; k < m; k++) s += a[k][i] * a[k][j]; - f = (s / a[i][i]) * g; - for (k = i; k < m; k++) a[k][j] += f * a[k][i]; - } - } - for (j = i; j < m; j++) a[j][i] *= g; - } - else { - for (j = i; j < m; j++) a[j][i] = 0.0; - } - ++a[i][i]; - } - - /* diagonalize the bidiagonal form */ - for (k = n - 1; k >= 0; k--) { /* loop over singular values */ - for (its = 0; its < 30; its++) { /* loop over allowed iterations */ - flag = 1; - for (l = k; l >= 0; l--) { /* test for splitting */ - nm = l - 1; - if (fabs(rv1[l]) + anorm == anorm) { - flag = 0; - break; - } - if (fabs(w[nm]) + anorm == anorm) break; - } - if (flag) { - c = 0.0; - s = 1.0; - for (i = l; i <= k; i++) { - f = s * rv1[i]; - if (fabs(f) + anorm != anorm) { - g = w[i]; - h = PYTHAG(f, g); - w[i] = h; - if (h == 0.0) { - char s[100]; - sprintf(s, "h = %f, f = %f, g = %f\n", h, f, g); - fprintf(stdout,"%s",s); - } - h = 1.0 / h; - c = g * h; - s = (- f * h); - for (j = 0; j < m; j++) { - y = a[j][nm]; - z = a[j][i]; - a[j][nm] = y * c + z * s; - a[j][i] = z * c - y * s; - } - } - } - } - z = w[k]; - if (l == k) { /* convergence */ - if (z < 0.0) { /* make singular value nonnegative */ - w[k] = -z; - for (j = 0; j < n; j++) v[j][k] = (-v[j][k]); - } - sort_sv(m, n, k, a, w, v); - break; - } - if (its >= 30) { - free_vector(rv1); - return(FALSE); /* return an error flag */ - } - - /* shift from bottom 2 x 2 minor */ - x = w[l]; - nm = k - 1; - y = w[nm]; - g = rv1[nm]; - h = rv1[k]; - f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y); - g = PYTHAG(f, 1.0); - f = ((x - z) * (x + z) + h * ((y / (f + SIGN(g, f))) - h)) / x; - - /* next QR transformation */ - c = s = 1.0; - for (j = l; j <= nm; j++) { - i = j + 1; - g = rv1[i]; - y = w[i]; - h = s * g; - g = c * g; - z = PYTHAG(f, h); - rv1[j] = z; - c = f / z; - s = h / z; - f = x * c + g * s; - g = g * c - x * s; - h = y * s; - y = y * c; - for (jj = 0; jj < n; jj++) { - x = v[jj][j]; - z = v[jj][i]; - v[jj][j] = x * c + z * s; - v[jj][i] = z * c - x * s; - } - z = PYTHAG(f, h); - w[j] = z; - if (z) { - z = 1.0 / z; - c = f * z; - s = h * z; - } - f = (c * g) + (s * y); - x = (c * y) - (s * g); - for (jj = 0; jj < m; jj++) { - y = a[jj][j]; - z = a[jj][i]; - a[jj][j] = y * c + z * s; - a[jj][i] = z * c - y * s; - } - } - rv1[l] = 0.0; - rv1[k] = f; - w[k] = x; - } - } - - free_vector(rv1); - return(TRUE); -} - diff --git a/lib/xlisp.h b/lib/xlisp.h deleted file mode 100644 index 073648b..0000000 --- a/lib/xlisp.h +++ /dev/null @@ -1,9 +0,0 @@ -#include - -#define FALSE 0 -#define TRUE 1 - -#ifndef IN_KCL_GLUE -typedef void *object; -typedef object LVAL; -#endif diff --git a/lib/xmath.h b/lib/xmath.h deleted file mode 100644 index 716e607..0000000 --- a/lib/xmath.h +++ /dev/null @@ -1,5 +0,0 @@ -#include - -#ifndef HUGE -#define HUGE 1e38 -#endif -- 2.11.4.GIT