From 3bdd28c4ae3de28dce32d8b9158c1f8d1b2e3924 Mon Sep 17 00:00:00 2001 From: AJ Rossini Date: Sun, 21 Sep 2008 20:12:34 +0200 Subject: [PATCH] Brent Fulgham's changes to source commited after verifying working C code with SBCL/i32 combination; BF suggested it works with Clozure on 64 bit. --- lib/Makefile.macosx | 6 ++- lib/betabase.c | 6 +-- lib/cbayes.c | 29 ++++-------- lib/cdists.c | 83 ++++++++++++++++------------------ lib/cffi-glue.c | 128 +++++++++++++++++++++++++++------------------------- lib/cfft.c | 51 ++++++++++----------- lib/cholesky.c | 4 +- lib/clinalg.c | 46 +++++++++---------- lib/derivatives.c | 16 +++++-- lib/functions.c | 93 ++++++++++++++++---------------------- lib/kernel.c | 14 ++++-- lib/linalgdata.c | 26 +++++------ lib/lowess.c | 36 +++++++++------ lib/ludecomp.c | 17 +++++-- lib/makerotation.c | 8 ++-- lib/minimize.c | 50 ++++++++++---------- lib/rcondest.c | 14 +++--- lib/splines.c | 22 ++++++--- lib/studentbase.c | 4 +- lib/svdecomp.c | 21 +++++---- 20 files changed, 348 insertions(+), 326 deletions(-) diff --git a/lib/Makefile.macosx b/lib/Makefile.macosx index e56e8b7..771a67b 100644 --- a/lib/Makefile.macosx +++ b/lib/Makefile.macosx @@ -1,4 +1,6 @@ ## -*- mode: Makefile -*- +ARCHS = -arch i386 -arch x86_64 +SPEC_WARNINGS = -Wconversion -Wformat -Wshorten-64-to-32 LISPSTAT-SO = liblispstat.dylib @@ -14,7 +16,7 @@ CBAYESOBJS = cbayes.o functions.o derivatives.o minimize.o CGLUE = cffi-glue.o ## CFLAGS = -g -DINTPTR -fPIC -Wall -CFLAGS = -g -Wall +CFLAGS = -g -Wall $(ARCHS) $(SPEC_WARNINGS) cffi : $(LISPSTAT-SO) @@ -29,7 +31,7 @@ clib.a: ${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) ranlib clib.a $(LISPSTAT-SO) :${CDISTOBJS} ${CLINALGOBJS} ${CBAYESOBJS} $(CGLUE) - gcc -dynamiclib -single_module -undefined dynamic_lookup -o $(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) *~ diff --git a/lib/betabase.c b/lib/betabase.c index f7eed68..7e5e6d4 100644 --- a/lib/betabase.c +++ b/lib/betabase.c @@ -10,7 +10,7 @@ extern double macheps(), gamma(); static double logbeta(), betai(), xinbta(); void -betabase(double *x, double *a, double *b, int *gia, int *gib, double *cdf) +betabase(double *x, double *a, double *b, long *gia, long *gib, double *cdf) { *cdf = betai(*x, *a, *b); } @@ -52,7 +52,7 @@ static double betai(double x, double pin, double qin) */ 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; - int i, n, ib; + long i, n, ib; /* I'm not sure these tolerances are optimal */ if (eps == 0.0) { @@ -171,7 +171,7 @@ static double xinbta(double *p, double *q, double *beta, double *alpha, int *ifa static int indx; static double prev, a, g, h, r, s, t, w, y, yprev, pp, qq; static double sq, tx, adj, acu; - static int iex; + static long iex; static double fpu, xin; /* Define accuracy and initialise. */ diff --git a/lib/cbayes.c b/lib/cbayes.c index f6402ff..8a842a4 100644 --- a/lib/cbayes.c +++ b/lib/cbayes.c @@ -10,8 +10,8 @@ typedef int PTR; typedef char *PTR; #endif -extern void maximize_callback(int, PTR, PTR, PTR, PTR, PTR); -extern void evalfront(char **, int *, double *, double *, double *, +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(); @@ -67,7 +67,7 @@ meminit() } static int -makespace(void **pptr, int size) /* why are we using **char? */ +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. */ @@ -91,9 +91,9 @@ makespace(void **pptr, int size) /* why are we using **char? */ /************************************************************************/ static void -callLminfun(int n, +callLminfun(size_t n, double *x, double *fval, double *grad, double *hess, - int *derivs) + long *derivs) { maximize_callback(n, (PTR) x, (PTR) fval, (PTR) grad, (PTR) hess, (PTR) derivs); @@ -103,7 +103,7 @@ void call_S(char *fun, long narg, char **args, char **mode, long *length,char **names, long nvals, char **values) { - int n = length[0], derivs; + long n = length[0], derivs; static double *fval = nil, *grad = nil, *hess = nil; makespace((void **)&fval, 1 * sizeof(double)); /* probably should test the @@ -135,10 +135,7 @@ Recover(s, w) /************************************************************************/ void -numgrad_front(n, px, pgrad, h, pscale) - int n; - PTR px, pgrad, pscale; - double h; +numgrad_front(size_t n, PTR px, PTR pgrad, double h, PTR pscale) { LVAL f = nil; double fval; @@ -148,10 +145,7 @@ numgrad_front(n, px, pgrad, h, pscale) } void -numhess_front(n, px, pf, pgrad, phess, h, pscale) - int n; - PTR px, pf, pgrad, phess, pscale; - double h; +numhess_front(size_t n, PTR px, PTR pf, PTR pgrad, PTR phess, double h, PTR pscale) { LVAL f = nil; @@ -181,8 +175,7 @@ numhess_front(n, px, pf, pgrad, phess, h, pscale) #define TSCAL_POS 10 #define MULT_POS 11 -static MaxIPars getMaxIPars(ipars) - int *ipars; +static MaxIPars getMaxIPars(int *ipars) { MaxIPars ip; @@ -219,9 +212,7 @@ static MaxDPars getMaxDPars(dpars) } void -minfo_maximize(px, pfvals, pscale, pip, pdp, verbose) - PTR px, pfvals, pscale, pip, pdp; - int verbose; +minfo_maximize(PTR px, PTR pfvals, PTR pscale, PTR pip, PTR pdp, int verbose) { LVAL f = nil; MaxIPars ip; diff --git a/lib/cdists.c b/lib/cdists.c index 7be4379..e8cad37 100644 --- a/lib/cdists.c +++ b/lib/cdists.c @@ -3,9 +3,9 @@ extern void normbase(double *,double *); extern void gammabase(double *,double *, double *); extern void studentbase(double *,double *, double *); -extern void betabase(double *, double *,double *,int *, int *,double *); -extern int max(int, int); -extern int min(int, int); +extern void betabase(double *, double *,double *,long *, long *,double *); +extern long max(long, long); +extern long min(long, long); extern void xlfail(char *); @@ -18,8 +18,13 @@ 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 @@ -30,13 +35,12 @@ extern double tdens(); #ifdef mips #ifdef UX31 -double tan(x) - double x; +double tan(double x) { return(sin(x) / cos(x)); } -double floor(x) - double x; + +double floor(double x) { long ix = x; double dx = ix; @@ -68,7 +72,7 @@ checkdf(double df) } static void -checkprob(double p, int zerostrict, int onestrict) +checkprob(double p, long zerostrict, long onestrict) { if (zerostrict) { if (p <= 0.0) xlfail("non-positive probability argument"); @@ -99,7 +103,7 @@ checkpoisson(double L) } static void -checkbinomial(int n, double p) +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"); @@ -132,8 +136,7 @@ double unirand() /*****************************************************************************/ /* standard normal cdf */ -double normalcdf(x) - double x; +double normalcdf(double x) { double p; normbase(&x, &p); @@ -141,8 +144,7 @@ double normalcdf(x) } /* standard normal quantile function */ -double normalquant(p) - double p; +double normalquant(double p) { int flag; double x; @@ -154,8 +156,7 @@ double normalquant(p) } /* standard normal density */ -double normaldens(x) - double x; +double normaldens(double x) { return(exp(- 0.5 * x * x) / sqrt(2.0 * PI)); } @@ -180,8 +181,7 @@ double normalrand() } /* bivariate normal cdf */ -double bnormcdf(x, y, r) - double x, y, r; +double bnormcdf(double x, double y, double r) { checkrho(r); return(bivnor(-x, -y, r)); @@ -280,7 +280,7 @@ gammarand(double a) { double x, u0, u1, u2, v, w, c, c1, c2, c3, c4, c5; static double e = -1.0; - int done; + long done; checkexp(a); if (e < 0.0) e = exp(1.0); @@ -385,7 +385,7 @@ double betacdf(double x, double a, double b) { double p; - int ia, ib; + long ia, ib; checkexp(a); checkexp(b); ia = a; ib = b; @@ -579,7 +579,7 @@ double frand(double ndf, double ddf) /*****************************************************************************/ static double -poisson_cdf(int k, double L) +poisson_cdf(long k, double L) { double dp, dx; @@ -601,20 +601,20 @@ double poissoncdf(double k, double L) { checkpoisson(L); - return(poisson_cdf((int) floor(k), L)); + return(poisson_cdf((long) floor(k), L)); } -int +long poissonquant(double x, double L) { - int k, k1, k2, del, ia; + 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, (int) (0.2 * s)); + del = max(1, (long) (0.2 * s)); if (x == 0.0) { k = 0.0; @@ -651,7 +651,7 @@ poissonquant(double x, double L) } double -poissonpmf(int k, double L) +poissonpmf(long k, double L) { double dx, dp; @@ -670,11 +670,11 @@ poissonpmf(int k, double L) } /* poisson random generator from Numerical Recipes */ -int poissonrand(double xm) +long poissonrand(double xm) { static double sqrt2xm, logxm, expxm, g, oldxm = -1.0; double t, y; - int k; + long k; checkpoisson(xm); if (xm < 12.0) { @@ -710,10 +710,10 @@ int poissonrand(double xm) /*****************************************************************************/ static double -binomial_cdf(int k, int n, double p) +binomial_cdf(long k, long n, double p) { double da, db, dp; - int ia, ib; + long ia, ib; if (k < 0) { dp = 0.0; @@ -741,17 +741,16 @@ binomial_cdf(int k, int n, double p) } double -binomialcdf(double k, int n, double p) +binomialcdf(double k, long n, double p) { checkbinomial(n, p); - return(binomial_cdf((int) floor(k), n, p)); - + return(binomial_cdf((long) floor(k), n, p)); } -int -binomialquant(double x, int n, double p) +long +binomialquant(double x, long n, double p) { - int k, k1, k2, del, ia; + long k, k1, k2, del, ia; double m, s, p1, p2, pk; checkbinomial(n, p); @@ -759,7 +758,7 @@ binomialquant(double x, int n, double p) m = n * p; s = sqrt(n * p * (1 - p)); - del = max(1, (int) (0.2 * s)); + del = max(1, (long) (0.2 * s)); if (x == 0.0) k = 0.0; else if (x == 1.0) k = n; @@ -787,9 +786,7 @@ binomialquant(double x, int n, double p) return(k2); } -double binomialpmf(k, n, p) - int k, n; - double p; +double binomialpmf(long k, long n, double p) { double dx, dp; @@ -804,12 +801,10 @@ double binomialpmf(k, n, p) } /* binomial random generator from Numerical Recipes */ -int binomialrand(n, pp) - int n; - double pp; +long binomialrand(long n, double pp) { - int j, k; - static int nold = -1; + 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; diff --git a/lib/cffi-glue.c b/lib/cffi-glue.c index c4232dd..12f2065 100644 --- a/lib/cffi-glue.c +++ b/lib/cffi-glue.c @@ -1,8 +1,14 @@ +/* #ifdef __APPLE__ #include #else #include #endif +*/ +#include +#include + + #include @@ -20,9 +26,10 @@ extern double betacdf(), betaquant(), betadens(), betarand(); extern double tcdf(), tquant(), tdens(), trand(); extern double fcdf(), fquant(), fdens(), frand(); extern double poissoncdf(), poissonpmf(); -extern int poissonquant(), poissonrand(); +extern long poissonquant(), poissonrand(); extern double binomialcdf(), binomialpmf(); -extern int binomialquant(), binomialrand(); +extern long binomialquant(), binomialrand(); +extern double rcondest_front(); void xlfail(char *); @@ -32,7 +39,6 @@ extern int lu_solve_front(); extern int lu_inverse_front(); extern int sv_decomp_front(); extern int qr_decomp_front(); -extern int rcondest_front(); extern int make_rotation_front(); extern int eigen_front(); extern int range_to_rseq(); @@ -53,10 +59,10 @@ extern int minfo_maximize(); /** Callback Support Functions **/ /***********************************************************************/ -static int ccl_integer_value; +static long ccl_integer_value; static double ccl_double_value; static PTR ccl_ptr_value; -void ccl_store_integer(int x) { ccl_integer_value = x; } +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; } @@ -71,10 +77,9 @@ static void (*free_ptr)(); register_new_ptr(f) void (*f)(); { new_ptr = f; } register_free_ptr(f) void (*f)(); { free_ptr = f; } -char *calloc(n, m) - int n, m; +char* calloc(size_t n, size_t m) { - int i, N = n * m; + size_t i, N = n * m; char *p; (*new_ptr)(N); @@ -86,8 +91,7 @@ char *calloc(n, m) malloc() { xlfail("malloc not available yet"); } realloc() { xlfail("realloc not available yet"); } -void free(p) - char *p; +void free(char *p) { (*free_ptr)(p); } @@ -100,25 +104,25 @@ void free(p) /***************************************************************************/ PTR -la_base_allocate(unsigned n, unsigned m) +la_base_allocate(size_t n, size_t m) { char *p = calloc(n, m); if (p == nil) xlfail("allocation failed"); return((PTR) p); } -int +long la_base_free_alloc(PTR p) { if (p) free((char *) p); return(0); } -int +size_t la_mode_size(int mode) { switch (mode) { - case IN: return(sizeof(int)); + case IN: return(sizeof(long)); case RE: return(sizeof(double)); case CX: return(sizeof(Complex)); } @@ -137,7 +141,7 @@ 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(int n, int m) +la_allocate(size_t n, size_t m) { (*ccl_la_allocate)(n, m); return(ccl_ptr_value); @@ -155,34 +159,34 @@ la_free_alloc(PTR p) /** **/ /***************************************************************************/ -int -la_get_integer(PTR p, int i) +long +la_get_integer(PTR p, size_t i) { - return(*(((int *) p) + i)); + return(*(((long *) p) + i)); } double -la_get_double(PTR p, int i) +la_get_double(PTR p, size_t i) { return(*(((double *) p) + i)); } double -la_get_complex_real(PTR p, int i) +la_get_complex_real(PTR p, size_t i) { Complex *c = ((Complex *) p) + i; return(c->real); } double -la_get_complex_imag(PTR p, int i) +la_get_complex_imag(PTR p, size_t i) { Complex *c = ((Complex *) p) + i; return(c->imag); } PTR -la_get_pointer(PTR p, int i) +la_get_pointer(PTR p, size_t i) { return(*(((PTR *) p) + i)); } @@ -194,20 +198,20 @@ la_get_pointer(PTR p, int i) /***************************************************************************/ int -la_put_integer(PTR p, int i, int x) +la_put_integer(PTR p, size_t i, long x) { - *(((int *) p) + i) = x; + *(((long *) p) + i) = x; return(0); } -int la_put_double(PTR p, int i, double x) +int la_put_double(PTR p, size_t i, double x) { - *(((double *) p) + i) = x; + *(((double *) p) + i) = x; return(0); } int -la_put_complex(PTR p, int i, double x, double y) +la_put_complex(PTR p, size_t i, double x, double y) { Complex *c = ((Complex *) p) + i; c->real = x; @@ -216,7 +220,7 @@ la_put_complex(PTR p, int i, double x, double y) } int -la_put_pointer(PTR p, int i, PTR x) +la_put_pointer(PTR p, size_t i, PTR x) { *(((PTR *) p) + i) = x; return(0); @@ -249,7 +253,7 @@ resetbuf() static void prbuf(char *s) { - int i, n; + size_t i, n; n = strlen(s); for (i = 0; i y) ? x : y); } @@ -61,26 +61,26 @@ double macheps() /************************************************************************/ void -chol_decomp_front(PTR mat, int n, PTR dpars) +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, int n, PTR iv, int mode, PTR dp) +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, int n, PTR indx, PTR b, int mode) +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, int n, PTR piv, PTR pv, int mode, PTR pinv) +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; @@ -118,31 +118,31 @@ lu_inverse_front(PTR pmat, int n, PTR piv, PTR pv, int mode, PTR pinv) } int -sv_decomp_front(PTR mat, int m, int n, PTR w, PTR v) +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, int m, int n, PTR v, PTR jpvt, int pivot) +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, int n) +rcondest_front(PTR mat, size_t n) { return(rcondest((char **) mat, n)); } void -make_rotation_front(int n, PTR rot, PTR x, PTR y, int use_alpha, double alpha) +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, int n, PTR w, PTR z, PTR fv1) +eigen_front(PTR a, size_t n, PTR w, PTR z, PTR fv1) { int ierr; @@ -151,23 +151,23 @@ eigen_front(PTR a, int n, PTR w, PTR z, PTR fv1) } void -fft_front(int n, PTR x, PTR work, int isign) +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, int n, double f, - int nsteps, double delta, PTR ys, PTR rw, PTR res) +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(int n, PTR px, int ns, PTR pxs) +range_to_rseq(size_t n, PTR px, long ns, PTR pxs) { - int i; + size_t i; double xmin, xmax, *x, *xs; x = (double *) px; @@ -181,7 +181,7 @@ range_to_rseq(int n, PTR px, int ns, PTR pxs) } int -spline_front(int n, PTR x, PTR y, int ns, +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, @@ -189,7 +189,7 @@ spline_front(int n, PTR x, PTR y, int ns, } int -kernel_dens_front(PTR x, int n, double width, PTR xs, PTR ys, int ns, int code) +kernel_dens_front(PTR x, size_t n, double width, PTR xs, PTR ys, size_t ns, int code) { int ktype; @@ -203,8 +203,8 @@ kernel_dens_front(PTR x, int n, double width, PTR xs, PTR ys, int ns, int code) } int -kernel_smooth_front(PTR x, PTR y, int n, double width, - PTR xs, PTR ys, int ns, int code) +kernel_smooth_front(PTR x, PTR y, size_t n, double width, + PTR xs, PTR ys, size_t ns, int code) { int ktype; diff --git a/lib/derivatives.c b/lib/derivatives.c index 6a70383..f7ff8c4 100644 --- a/lib/derivatives.c +++ b/lib/derivatives.c @@ -4,6 +4,14 @@ /* 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(); @@ -12,12 +20,12 @@ extern double macheps(); # define TRUE 1 void -numergrad(int n, double *x, double *grad, double *fsum, +numergrad(size_t n, double *x, double *grad, double *fsum, int ffun(double *, double *, double *, double **), double h, double *typx) /* int (*ffun)();*/ { - int i; + size_t i; double old_xi, f1, f2, hi; for (i = 0; i < n; i++) { @@ -34,11 +42,11 @@ numergrad(int n, double *x, double *grad, double *fsum, } void -numerhess(int n, double *x, double **hess, double f, double *fsum, +numerhess(size_t n, double *x, double **hess, double f, double *fsum, int *ffun(double *, double *, double *, double *), double h, double *typx) { - int i, j; + size_t i, j; double old_xi, old_xj, f1, f2, hi, hj; for (i = 0; i < n; i++) { diff --git a/lib/functions.c b/lib/functions.c index a6b42ec..6d23f8d 100644 --- a/lib/functions.c +++ b/lib/functions.c @@ -20,7 +20,7 @@ extern void minsupplyvalues(double , double *, double **, extern void minimize(); extern void minresults(double *, double *, double *, double *, double **, double *, double **, int *, int *, double *); -extern void minsetup(int n, int k, +extern void minsetup(size_t n, size_t k, int *ffun(), int *gfun(), double *x, double typf, double *typx, char *work); /* int (*ffun)(), (*gfun)();*/ @@ -37,22 +37,22 @@ extern void call_S(char *fun, long narg, char **args, char **mode, long *length, /* next 2 from derivatives.c */ /* -extern void numergrad(int n, double *x, double *grad, double *fsum, +extern void numergrad(size_t n, double *x, double *grad, double *fsum, int *ffun(), double h, double *typx); -extern void numerhess(int n, double *x, double **hess, double f, +extern void numerhess(size_t n, double *x, double **hess, double f, double *fsum, void ffun(), double h, double *typx); */ -extern void numergrad(int , double *, double *, double *, +extern void numergrad(size_t , double *, double *, double *, int ffun(double *, double *, double *, double **), double , double *); -extern void numerhess(int , double *, double **, double , +extern void numerhess(size_t , double *, double **, double , double *, int ffun(double *, double *, double *, double **), double , double *); /* next from minimize */ -extern int minworkspacesize(int, int); +extern size_t minworkspacesize(size_t, size_t); /************************************************************************/ /** **/ @@ -76,7 +76,7 @@ extern int minworkspacesize(int, int); typedef struct{ char *f, **sf, **g; - int n, k; + size_t n, k; int change_sign, fderivs; int *gderivs; double typf, h, dflt; @@ -99,9 +99,7 @@ static Fundata func, gfuncs, cfuncs; /* from S. It attempts to avoid the danger of dangling callocs. */ void static -makespace(pptr, size) - char **pptr; - int size; +makespace(char **pptr, size_t size) { if (size <= 0) return; if (*pptr == nil) *pptr = calloc(size, 1); @@ -122,7 +120,7 @@ makespace(pptr, size) /* install log posterior function */ static void -install_func(char *f, char **sf, int n, int change_sign, /*?? */ +install_func(char *f, char **sf, size_t n, int change_sign, /*?? */ double typf, double h, double *typx, double dflt) { int i; @@ -150,9 +148,9 @@ install_func(char *f, char **sf, int n, int change_sign, /*?? */ /* install tilt functions */ static void -install_gfuncs(char **g, int n, int k, int change_sign, double h, double *typx) +install_gfuncs(char **g, size_t n, size_t k, int change_sign, double h, double *typx) { - int i; + size_t i; static int inited = FALSE; static double *gfsumdata = nil; @@ -180,9 +178,9 @@ install_gfuncs(char **g, int n, int k, int change_sign, double h, double *typx) /* install constraint functions */ static void -install_cfuncs(char **g, int n, int k, double *ctarget, double h, double *typx) +install_cfuncs(char **g, size_t n, size_t k, double *ctarget, double h, double *typx) { - int i; + size_t i; static int inited = FALSE; if (! inited) { @@ -206,7 +204,7 @@ install_cfuncs(char **g, int n, int k, double *ctarget, double h, double *typx) } static int -in_support(char **ff, int n, double *x) +in_support(char **ff, size_t n, double *x) { char *args[1], *values[1]; int *result; @@ -287,7 +285,7 @@ evalfunc(double *x, double *pval, double *grad, double **hess) /* callback for tilt function evaluation */ -static int which_gfunc; +static size_t which_gfunc; static int evalgfunc(double *x, double *pval, double *grad, double **hess) @@ -399,10 +397,10 @@ evalcfunc(double *x, double *pval, double *grad, double **hess) /* S front end for logposterior evaluation */ void -evalfront(char **ff, int *n, double *x, double *val, double *grad, +evalfront(char **ff, size_t *n, double *x, double *val, double *grad, double *phess, double *h, double *typx) { - int i; + size_t i; static double **hess = nil; install_func(ff[0], nil, *n, FALSE, 1.0, *h, typx, 0.0); @@ -416,10 +414,10 @@ evalfront(char **ff, int *n, double *x, double *val, double *grad, /* S front end for tilt function evaluation */ void -gevalfront(char **gg, int *n, int *m, double *x, double *h, +gevalfront(char **gg, size_t *n, size_t *m, double *x, double *h, double *typx, double *val, double *grad) { - int i; + size_t i; install_gfuncs(gg, *n, *m, FALSE, *h, typx); for (i = 0; i < *m; i++, val++) { @@ -483,9 +481,9 @@ check_derivs(double *x, double drvtol) /* joint density of normal-cauchy mixture */ static double -dncmix(double *x, int n, double p) +dncmix(double *x, size_t n, double p) { - int i; + size_t i; double dens; for (i = 0, dens = 1.0; i < n; i++) { @@ -501,7 +499,7 @@ dncmix(double *x, int n, double p) */ void samplefront(char **ff, char **sf, char **rf, - double *p, int *n, + double *p, size_t *n, double *x, double *ch, int *N, double *y, double *w) { double val; @@ -564,7 +562,7 @@ static void add_tilt(double *x, double *pval, double *grad, double **hess, double tilt, int exptilt) { - int i, j, k, n = func.n, m = gfuncs.k; + size_t i, j, k, n = func.n, m = gfuncs.k; double *gval, *ggrad, **ghess, etilt; if (m == 0) return; @@ -607,11 +605,11 @@ add_tilt(double *x, double *pval, double *grad, double **hess, } static void -set_tilt_info(int n, int m, +set_tilt_info(size_t n, size_t m, double tilt, int exptilt, double *tscale) { static double *hessdata = nil, *graddata = nil; - int i; + size_t i; static int inited = FALSE; if (! inited) { @@ -637,7 +635,7 @@ set_tilt_info(int n, int m, static void minfunc(double *x, double *pval, double *grad, double **hess) { - int k = gfuncs.k; + size_t k = gfuncs.k; if (evalfunc(x, pval, grad, hess) && (k > 0)) add_tilt(x, pval, grad, hess, tiltinfo.tilt, tiltinfo.exptilt); @@ -666,7 +664,7 @@ maxfront(char **ff, char **gf, char **cf, static char *work = nil; static double **H = nil, **cJ = nil; double *pf, *grad, *c; - int i, n, m, k; + size_t i, n, m, k; void (*cfun)(); if (ipars->verbose > 0) PRINTSTR("maximizing...\n"); @@ -797,10 +795,8 @@ loglapfront(fvals, cvals, ipars, dpars, val) /** **/ /************************************************************************/ -moms1front(gf, n, m, x, hessdata, mean, stdev, sigmadata, h, typx) - char *gf; - int *n, *m; - double *x, *hessdata, *mean, *stdev, *sigmadata, *h, *typx; +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; @@ -855,13 +851,9 @@ typedef struct { double mgfdel; } MomDPars; -moms2front(ff, gf, gnum, x, typx, fvals, gvals, ipars, dpars, - mean, stdev, sigmadata) - char **ff, **gf; - int *gnum; - double *x, *typx, *fvals, *gvals, *mean, *stdev, *sigmadata; - MomIPars *ipars; - MomDPars *dpars; +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; @@ -1007,10 +999,8 @@ moms2front(ff, gf, gnum, x, typx, fvals, gvals, ipars, dpars, } } -static copypars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1) - double *x, *f, *g, *x1, *f1, *g1; - MomIPars *ip, *ip1; - MomDPars *dp, *dp1; +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; @@ -1032,12 +1022,9 @@ static copypars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1) /** **/ /************************************************************************/ -lapmar1front(ff, gf, x, typx, fvals, ipars, dpars, xmar, ymar, nmar) - char **ff, **gf; - double *x, *typx, *fvals, *xmar, *ymar; - MaxIPars *ipars; - MaxDPars *dpars; - int *nmar; +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; @@ -1096,10 +1083,8 @@ lapmar1front(ff, gf, x, typx, fvals, ipars, dpars, xmar, ymar, nmar) } } -static cpmarpars(x, f, g, ip, dp, x1, f1, g1, ip1, dp1) - double *x, *f, *g, *x1, *f1, *g1; - MaxIPars *ip, *ip1; - MaxDPars *dp, *dp1; +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; diff --git a/lib/kernel.c b/lib/kernel.c index e00550f..5504a47 100644 --- a/lib/kernel.c +++ b/lib/kernel.c @@ -8,6 +8,14 @@ #define nil 0L #endif /*nil */ +/* +#ifdef __APPLE__ +#include +#endif +*/ + +#include + static double kernel(double x, double y, double w, int type) { @@ -50,11 +58,11 @@ kernel(double x, double y, double w, int type) } int -kernel_smooth(double *x, double *y, int n, +kernel_smooth(double *x, double *y, size_t n, double width, double *wts, double *wds, - double *xs, double *ys, int ns, int ktype) + double *xs, double *ys, size_t ns, int ktype) { - int i, j; + size_t i, j; double wsum, ysum, lwidth, lwt, xmin, xmax; if (n < 1) return(1); diff --git a/lib/linalgdata.c b/lib/linalgdata.c index 85a9e68..db1bcc7 100644 --- a/lib/linalgdata.c +++ b/lib/linalgdata.c @@ -9,11 +9,10 @@ typedef int PTR; typedef char *PTR; #endif -extern PTR la_allocate(int, int); +extern PTR la_allocate(size_t, size_t); extern void la_free_alloc(PTR); extern void xlfail(char *); - /************************************************************************/ /** **/ /** Storage Allocation Functions **/ @@ -21,7 +20,7 @@ extern void xlfail(char *); /************************************************************************/ static char -*allocate(unsigned n, unsigned m) +*allocate(size_t n, size_t m) { char *p = (char *) la_allocate(n, m); if (p == nil) xlfail("allocation failed"); @@ -35,19 +34,18 @@ free_alloc(char *p) } IVector -ivector(unsigned n) +ivector(size_t n) { return((IVector) allocate(n, sizeof(int))); } double -*rvector(unsigned n) +*rvector(size_t n) { return((double *) allocate(n, sizeof(double))); } -CVector cvector(n) - unsigned n; +CVector cvector(size_t n) { return((CVector) allocate(n, sizeof(Complex))); } @@ -58,27 +56,27 @@ free_vector(double *v) free_alloc((char *)v); } -IMatrix imatrix(unsigned n, unsigned m) +IMatrix imatrix(size_t n, size_t m) { - int i; + 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(unsigned n, unsigned m) +rmatrix(size_t n, size_t m) { - int i; + 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(unsigned n, unsigned m) +cmatrix(size_t n, size_t m) { - int i; + 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); @@ -87,7 +85,7 @@ cmatrix(unsigned n, unsigned m) void free_matrix(double **mat, int n) /* Matrix?? Not RMatrix?? */ { - int i; + 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 index 01ea67b..c0780f2 100644 --- a/lib/lowess.c +++ b/lib/lowess.c @@ -3,31 +3,40 @@ */ #include "xmath.h" +/* +#ifdef __APPLE__ +#include +#endif +*/ +#include #define FALSE 0 #define TRUE 1 -extern int max(int, int); -extern int min(int, int); +extern long max(long, long); +extern long min(long, long); -static double pow2(x) double x; { return(x * x); } -static double pow3(x) double x; { return(x * x * x); } -/* static */ double fmax(x,y) double x, y; { return (x > y ? x : y); } +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(double *a, double *b) +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, int n, double xs, double *ys, int nleft, int nright, +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; - int j, nrt; + long j, nrt; range = x[n - 1] - x[0]; h = fmax(xs - x[nleft], x[nright] - xs); @@ -76,7 +85,7 @@ lowest(double *x, double *y, int n, double xs, double *ys, int nleft, int nright } static void -sort(double *x, int n) +sort(double *x, size_t n) { extern void qsort(); @@ -86,15 +95,16 @@ sort(double *x, int n) int -lowess(double *x, double *y, int n, - double f, int nsteps, +lowess(double *x, double *y, size_t n, + double f, size_t nsteps, double delta, double *ys, double *rw, double *res) { - int iter, ns, ok, nleft, nright, i, j, last, m1, m2; + 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((int) (f * n), n), 2); /* at least two, at most n points */ + 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 */ diff --git a/lib/ludecomp.c b/lib/ludecomp.c index 19fa1d0..e659d1c 100644 --- a/lib/ludecomp.c +++ b/lib/ludecomp.c @@ -8,9 +8,10 @@ #include "linalg.h" int -crludcmp(Matrix mat, int n, IVector indx, int mode, double *d) +crludcmp(Matrix mat, size_t n, IVector indx, int mode, double *d) { - int i, imax, j, k, singular = FALSE; + int i, imax, j, k; + int singular = FALSE; double big, temp; Complex cdum, csum; double rdum, rsum; @@ -44,9 +45,11 @@ crludcmp(Matrix mat, int n, IVector indx, int mode, double *d) if (mode == RE) rsum = rmat[i][j]; else csum = cmat[i][j]; - for (k = 0; k < i; k++) + 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; @@ -57,8 +60,10 @@ crludcmp(Matrix mat, int n, IVector indx, int mode, double *d) 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; @@ -96,10 +101,14 @@ crludcmp(Matrix mat, int n, IVector indx, int mode, double *d) } 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); + for (i = j + 1; i < n; i++) + { + cmat[i][j] = cmul(cmat[i][j], cdum); + } } } } + free_vector(vv); return(singular); } diff --git a/lib/makerotation.c b/lib/makerotation.c index 006d53b..6b7decc 100644 --- a/lib/makerotation.c +++ b/lib/makerotation.c @@ -6,9 +6,7 @@ #include "linalg.h" -static double inner_product(n, x, y) - int n; - RVector x, y; +static double inner_product(size_t n, RVector x, RVector y) { double result = 0.0; @@ -19,10 +17,10 @@ static double inner_product(n, x, y) #define NORM(n, x) (sqrt(inner_product(n, x, x))) void -make_rotation(int n, double **rot, double *x, double *y, int use_alpha, double alpha) +make_rotation(size_t n, double **rot, double *x, double *y, int use_alpha, double alpha) { double nx, ny, xy, c, s; - int i, j; + size_t i, j; for (i = 0; i < n; i++) { for (j = 0; j < n; j++) rot[i][j] = 0.0; diff --git a/lib/minimize.c b/lib/minimize.c index e6440c8..05f5028 100644 --- a/lib/minimize.c +++ b/lib/minimize.c @@ -16,7 +16,7 @@ extern char buf[200]; #define PRINTSTR(s) printf(s) typedef struct { - int n, k; + size_t n, k; void (*ffun)(), (*gfun)(); double f, typf, new_f; double crit, new_crit; @@ -114,9 +114,9 @@ static double Min(a, b) /* solve Ly = b for y */ static void -lsolve(int n, double *b, double **L, double *y) +lsolve(size_t n, double *b, double **L, double *y) { - int i, j; + size_t i, j; for (i = 0; i < n; i++) { y[i] = b[i]; @@ -127,9 +127,9 @@ lsolve(int n, double *b, double **L, double *y) /* solve (L^T)x = y for x */ static void -ltsolve(int n, double *y, double **L, double *x) +ltsolve(size_t n, double *y, double **L, double *x) { - int i, j; + size_t i, j; for (i = n - 1; i >= 0; i--) { x[i] = y[i]; @@ -140,11 +140,11 @@ ltsolve(int n, double *y, double **L, double *x) /* solve (L L^T) s = -g for s */ static void -cholsolve(int n, +cholsolve(size_t n, double *g, double **L, double *s) { - int i; + size_t i; /* solve Ly = g */ lsolve(n, g, L, s); @@ -156,9 +156,9 @@ cholsolve(int n, } static void -modelhess(int n, double *sx, double **H, double **L, double *diagadd) +modelhess(size_t n, double *sx, double **H, double **L, double *diagadd) { - int i, j; + size_t i, j; double sqrteps, maxdiag, mindiag, maxoff, maxoffl, maxposdiag, mu, maxadd, maxev, minev, offrow, sdd; @@ -252,7 +252,7 @@ modelhess(int n, double *sx, double **H, double **L, double *diagadd) static double gradsize(Iteration *iter, int new) { - int n, i; + size_t n, i; double size, term, crit, typf; double *x, *delf, *sx; @@ -273,7 +273,7 @@ gradsize(Iteration *iter, int new) static double incrsize(Iteration *iter) { - int n, i; + size_t n, i; double size, term; double *x, *new_x, *sx; @@ -369,7 +369,7 @@ eval_next_funval(Iteration *iter) static void eval_gradient(Iteration *iter) { - int i, j, n, k; + size_t i, j, n, k; n = iter->n; k = iter->k; @@ -392,7 +392,7 @@ eval_gradient(Iteration *iter) static void eval_next_gradient(Iteration *iter) { - int i, j, n, k; + size_t i, j, n, k; n = iter->n; k = iter->k; @@ -409,7 +409,7 @@ eval_next_gradient(Iteration *iter) static void eval_hessian(Iteration *iter) { - int i, j, n, k; + size_t i, j, n, k; n = iter->n; k = iter->k; @@ -427,7 +427,7 @@ static void eval_hessian(Iteration *iter) static void linesearch(Iteration *iter) { - int i, n; + size_t i, n; double newtlen, maxstep, initslope, rellength, lambda, minlambda, lambdatemp, lambdaprev, a, b, disc, critprev, f1, f2, a11, a12, a21, a22, del; @@ -537,7 +537,7 @@ print_header(Iteration *iter) static void print_status(Iteration *iter) { - int i, j; + size_t i, j; if (iter->verbose > 0) { sprintf(buf, "Criterion value = %g\n", @@ -584,7 +584,7 @@ print_status(Iteration *iter) static void findqnstep(Iteration *iter) { - int i, j, N, l; + size_t i, j, N, l; if (iter->k == 0) { modelhess(iter->n, iter->sx, iter->hessf, iter->L, &iter->diagadd); @@ -606,7 +606,7 @@ findqnstep(Iteration *iter) static void iterupdate(Iteration *iter) { - int i, j, n, k; + size_t i, j, n, k; n = iter->n; k = iter->k; @@ -660,10 +660,10 @@ mindriver(Iteration *iter) static Iteration myiter; -int -minworkspacesize(int n, int k) +size_t +minworkspacesize(size_t n, size_t k) { - int size; + size_t size; /* x, new_x, sx, delf, new_delf, qnstep and F */ size = 7 * sizeof(double) * (n + k); @@ -687,13 +687,13 @@ minresultstring(int code) } void -minsetup(int n, int k, +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; - int i, j; + size_t i, j; double nx0, ntypx; n = (n > 0) ? n : 0; @@ -781,7 +781,7 @@ minsupplyvalues(double f, double *delf, double **hessf, double *g, double **delg) { Iteration *iter = &myiter; - int i, j, n, k; + size_t i, j, n, k; n = iter->n; k = iter->k; @@ -819,7 +819,7 @@ minresults(double *x, double *pf, double *pcrit, double *delf, double **hessf, double *g, double **delg, int *pcount, int *ptermcode, double *diagadd) { Iteration *iter = &myiter; - int i, j, n, k; + size_t i, j, n, k; n = iter->n; k = iter->k; diff --git a/lib/rcondest.c b/lib/rcondest.c index ccb0e31..03c5e17 100644 --- a/lib/rcondest.c +++ b/lib/rcondest.c @@ -14,9 +14,9 @@ double Max(double a, double b) void -backsolve(double **a, double *b, int n, int mode) +backsolve(double **a, double *b, size_t n, int mode) { - int i, j; + long i, j; /* Must be signed, since counting down in loops... */ RMatrix ra = (RMatrix) a; RVector rb = (RVector) b; CMatrix ca = (CMatrix) a; @@ -41,12 +41,12 @@ backsolve(double **a, double *b, int n, int mode) /** Exported function **/ -double rcondest(double **a, int n) +double rcondest(double **a, size_t n) { RVector p, pm, x; double est, xp, xm, temp, tempm, xnorm; - int i, j; - + long i, j; + for (i = 0; i < n; i++) if (a[i][i] == 0.0) return(0.0); @@ -56,6 +56,7 @@ double rcondest(double **a, int 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]); @@ -87,6 +88,7 @@ double rcondest(double **a, int n) 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; @@ -94,7 +96,7 @@ double rcondest(double **a, int n) free_vector(p); free_vector(pm); free_vector(x); - + return(est); } diff --git a/lib/splines.c b/lib/splines.c index b9c95b4..585d2c7 100644 --- a/lib/splines.c +++ b/lib/splines.c @@ -1,13 +1,21 @@ #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, int n, +find_spline_derivs(double *x, double *y, size_t n, double *y2, double *u) { - int i, k; + long i, k; double p, sig; y2[0] = u[0] = 0.0; /* lower boundary condition for natural spline */ @@ -38,9 +46,9 @@ find_spline_derivs(double *x, double *y, int n, /* interpolate or extrapolate value at x using results of find_spline_derivs */ static void spline_interp(double *xa, double *ya, double *y2a, - int n, double x, double *y) + size_t n, double x, double *y) { - int klo, khi, k; + long klo, khi, k; double h, b, a; if (x <= xa[0]) { @@ -83,10 +91,10 @@ spline_interp(double *xa, double *ya, double *y2a, } int -fit_spline(int n, double *x, double *y, - int ns, double *xs, double *ys, double *work) +fit_spline(size_t n, double *x, double *y, + size_t ns, double *xs, double *ys, double *work) { - int i; + size_t i; double *y2, *u; y2 = work; u = work + n; diff --git a/lib/studentbase.c b/lib/studentbase.c index 907f6c6..4a6d97b 100644 --- a/lib/studentbase.c +++ b/lib/studentbase.c @@ -4,7 +4,7 @@ #define HALF_PI 1.5707963268 #define EPSILON .000001 -extern void betabase(double *, double *, double *,int *, int *,double *); +extern void betabase(double *, double *, double *,long *, long *,double *); extern void normbase(double *,double *); extern double ppnd(), ppbeta(); @@ -26,7 +26,7 @@ studentbase(double *x, double *df, double *cdf) if (n < 2.0 && n != 1.0) { /* beta integral aproximation for small df */ double da = 0.5, db = 0.5 * n, dx, dp; - int ia = 0, ib = floor(db); + long ia = 0, ib = floor(db); dx = db / (db + da * t); betabase(&dx, &db, &da, &ib, &ia, &dp); diff --git a/lib/svdecomp.c b/lib/svdecomp.c index ed8aad2..5be0073 100644 --- a/lib/svdecomp.c +++ b/lib/svdecomp.c @@ -7,8 +7,9 @@ #include "linalg.h" -static double PYTHAG(a, b) - double a, b; +#include + +static double PYTHAG(double a, double b) { double at = fabs(a), bt = fabs(b), ct, result; @@ -21,10 +22,10 @@ static double PYTHAG(a, b) #define SWAPD(a, b) (temp = (a), (a) = (b), (b) = temp) static void -sort_sv(int m, int n, int k, +sort_sv(size_t m, size_t n, size_t k, double **a, double *w, double **v) { - int i, j; + size_t i, j; double temp; for (i = k; (i < n - 1) && (w[i] < w[i+1]); i++) { @@ -39,15 +40,17 @@ static double maxarg1, maxarg2; #define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) int -svdcmp(double **a, int m, int n, double *w, double **v) +svdcmp(double **a, size_t m, size_t n, double *w, double **v) { - int flag, i, its, j, jj, k, l, nm; + 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 */ @@ -112,6 +115,7 @@ svdcmp(double **a, int m, int n, double *w, double **v) 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]; @@ -252,6 +256,7 @@ svdcmp(double **a, int m, int n, double *w, double **v) w[k] = x; } } + free_vector(rv1); return(TRUE); } -- 2.11.4.GIT