/* ccmath.h CCMATH mathematics library source code.
*
* Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
* This code may be redistributed under the terms of the GNU library
* public license (LGPL). ( See the lgpl.license file for details.)
* ------------------------------------------------------------------------
*/
/*
CCM
Numerical Analysis Toolkit Header File
ELF Shared Library Version
*/
/* Required for Shared Library */
#define XMATH 1
/* Define File Pointers and Standard Library */
#include <stdio.h>
#include <stdlib.h>
/* Definitions of Types */
#ifndef NULL
#define NULL ((void *)0
#endif
/* Complex Types */
#ifndef CPX
struct complex {double re,im;};
typedef struct complex Cpx;
#define CPX 1
#endif
/* Orthogonal Polynomial Type */
#ifndef OPOL
struct orpol {double cf,hs,df;};
typedef struct orpol Opol;
#define OPOL 1
#endif
/* Tree Types */
#ifdef BAL
struct tnode {char *key,*rec; int bal; struct tnode *pr,*pl;};
#else
struct tnode {char *key,*rec; struct tnode *pr,*pl;};
#endif
/* Time Series Types */
struct mcof {double cf; int lag;};
struct fmod {int fac; double val;};
/* List Definition */
struct llst {char *pls; struct llst *pt;};
/* Hash Table Definition */
struct tabl {char *key,*val; struct tabl *pt;};
/* Extended Precision Types */
/* XMATH must be defined to use extended precision functions */
#ifdef XMATH
#ifndef XPR
#define XDIM 7
struct xpr {unsigned short nmm[XDIM+1];};
extern unsigned short m_sgn,m_exp;
extern short bias;
extern int itt_div,k_tanh;
extern int ms_exp,ms_trg,ms_hyp;
extern short max_p,k_lin;
extern short d_bias,d_max,d_lex;
extern struct xpr zero,one,two,ten;
extern struct xpr x_huge;
/* Variables used in extended precision arithmetic */
unsigned short m_sgn=0x8000,m_exp=0x7fff;
short bias=16383;
int itt_div=2,k_tanh=5;
int ms_exp=21,ms_hyp=25,ms_trg=31;
short max_p=16*XDIM,k_lin= -8*XDIM;
short d_bias=15360,d_max=2047,d_lex=12;
struct xpr zero={{0x0,0x0}};
struct xpr one={{0x3fff,0x8000}};
struct xpr two={{0x4000,0x8000}};
struct xpr ten={{0x4002,0xa000}};
struct xpr x_huge={{0x7fff,0x0}};
/* Variables used in the extended precision math functions */
struct xpr pi4={{0x3FFE,0xC90F,0xDAA2,0x2168,0xC234,0xC4C6,0x628B,0x80DC}};
struct xpr pi2={{0x3FFF,0xC90F,0xDAA2,0x2168,0xC234,0xC4C6,0x628B,0x80DC}};
struct xpr pi={{0x4000,0xC90F,0xDAA2,0x2168,0xC234,0xC4C6,0x628B,0x80DC}};
struct xpr ee={{0x4000,0xADF8,0x5458,0xA2BB,0x4A9A,0xAFDC,0x5620,0x273D}};
struct xpr ln2={{0x3FFE,0xB172,0x17F7,0xD1CF,0x79AB,0xC9E3,0xB398,0x3F3}};
struct xpr srt2={{0x3FFF,0xB504,0xF333,0xF9DE,0x6484,0x597D,0x89B3,0x754B}};
#define XPR 1
#endif
#endif
/* FUNCTION DECLARATIONS */
/* Linear Algebra */
/* Real Linear Systems */
int minv(double *a,int n) ;
int psinv(double *v,int n) ;
int ruinv(double *a,int n) ;
int solv(double *a,double *b,int n) ;
int solvps(double *s,double *x,int n) ;
int solvru(double *a,double *b,int n) ;
void solvtd(double *a,double *b,double *c,double *x,int m) ;
void eigen(double *a,double *eval,int n) ;
void eigval(double *a,double *eval,int n) ;
double evmax(double *a,double *u,int n) ;
int svdval(double *d,double *a,int m,int n) ;
int sv2val(double *d,double *a,int m,int n) ;
int svduv(double *d,double *a,double *u,int m,double *v,int n) ;
int sv2uv(double *d,double *a,double *u,int m,double *v,int n) ;
int svdu1v(double *d,double *a,int m,double *v,int n) ;
int sv2u1v(double *d,double *a,int m,double *v,int n) ;
void mmul(double *mat,double *a,double *b,int n) ;
void rmmult(double *mat,double *a,double *b,int m,int k,int n) ;
void vmul(double *vp,double *mat,double *v,int n) ;
double vnrm(double *u,double *v,int n) ;
void matprt(double *a,int n,int m,char *fmt) ;
void fmatprt(FILE *fp,double *a,int n,int m,char *fmt) ;
void trnm(double *a,int n) ;
void mattr(double *a,double *b,int m,int n) ;
void otrma(double *at,double *u,double *a,int n) ;
void otrsm(double *st,double *u,double *s0,int n) ;
void mcopy(double *a,double *b,int m) ;
void ortho(double *evc,int n) ;
void smgen(double *a,double *eval,double *evec,int n) ;
/* utility routines for real symmertic eigensystems */
void house(double *a,double *d,double *ud,int n) ;
void housev(double *a,double *d,double *ud,int n) ;
int qreval(double *eval,double *ud,int n) ;
int qrevec(double *eval,double *evec,double *dp,int n) ;
/* utility routines for singular value decomposition */
int qrbdi(double *d, double *e,int n) ;
int qrbdv(double *d, double *e,double *u,int m,double *v,int n) ;
int qrbdu1(double *d, double *e,double *u,int m,double *v,int n) ;
void ldumat(double *a,double *u,int m,int n) ;
void ldvmat(double *a,double *v,int n) ;
void atou1(double *a,int m,int n) ;
void atovm(double *v,int n) ;
/* Complex Matrix Algebra */
int cminv(Cpx *a,int n) ;
int csolv(Cpx *a,Cpx *b,int n) ;
void heigvec(Cpx *a,double *eval,int n) ;
void heigval(Cpx *a,double *eval,int n) ;
double hevmax(Cpx *a,Cpx *u,int n) ;
void cmmul(Cpx *c,Cpx *a,Cpx *b,int n) ;
void cmmult(Cpx *c,Cpx *a,Cpx *b,int m,int k,int n) ;
void cvmul(Cpx *vp,Cpx *mat,Cpx *v,int n) ;
Cpx cvnrm(Cpx *u,Cpx *v,int n) ;
void cmprt(Cpx *a,int n,int m,char *fmt) ;
void trncm(Cpx *a,int n) ;
void hconj(Cpx *u,int n) ;
void cmattr(Cpx *a,Cpx *b,int m,int n) ;
void utrncm(Cpx *at,Cpx *u,Cpx *a,int n) ;
void utrnhm(Cpx *ht,Cpx *u,Cpx *h0,int n) ;
void cmcpy(Cpx *a,Cpx *b,int n) ;
void unitary(Cpx *u,int n) ;
void hmgen(Cpx *h,double *eval,Cpx *u,int n) ;
/* utility routines for hermitian eigen problems */
void chouse(Cpx *a,double *d,double *ud,int n) ;
void chousv(Cpx *a,double *d,double *ud,int n) ;
void qrecvc(double *eval,Cpx *evec,double *ud,int n) ;
/* Geometry */
void crossp(double *h,double *u,double *v) ;
double dotp(double *u,double *v,int m) ;
double metpr(double *u,double *a,double *v,int n) ;
void scalv(double *r,double s,int n) ;
void trvec(double *c,double *a,double *b,int n) ;
double leng(double *a,double *b,int n) ;
void rotax(double *v,double az,double pa,double ang,int k) ;
void euler(double *pv,int m,double a,double b,double c) ;
/* plane trigonometry */
void trgsas(double a,double g,double b,double *ans);
int trgasa(double a,double ss,double b,double *asn);
double trgarea(double a,double b,double c);
int trgsss(double a,double b,double c,double *ang);
int trgssa(double a,double b,double ba,double *an);
/* spherical trigonometry */
void stgsas(double a,double g,double b,double *ang);
int stgasa(double a,double c,double b,double *ang);
int stgsss(double a,double b,double c,double *ang);
int stgaaa(double a,double b,double c,double *ang);
double stgarea(double a,double b,double c);
/* hyperbolic trigonometry */
void htgsas(double a,double g,double b,double *an);
int htgasa(double a,double cc,double b,double *ans);
int htgsss(double a,double b,double c,double *ang);
int htgaaa(double a,double b,double c,double *as);
double htgarea(double a,double b,double c);
/* Numerical Integration */
double fintg(double a,double b,int n,double te,double (*func)()) ;
/* functional form: double (*func)(double) */
double chintg(double *a,int m,double (*func)()) ;
/* functional form: double (*func)(double) */
double fchb(double x,double *a,int m) ;
int deqsy(double *y,int n,double a,double b,int nd,double te,
int (*fsys)()) ;
/* functional form: int (*fsys)(double x,double *y,double *dp) */
/* Optimization and Roots */
int optmiz(double *x,int n,double (*func)(),double de,
double test,int max) ;
/* functional form: double (*func)(double *x) */
double optsch(double (*func)(),double a,double b,double test) ;
/* functional form: double (*func)(double) */
int plrt(double *cof,int n,struct complex *root,double ra,double rb) ;
struct complex polyc(struct complex z,double *cof,int n) ;
double secrt(double (*func)(),double x,double dx,double test) ;
/* functional form: double (*func)(double) */
int solnl(double *x,double *f,double (*fvec[])(),int n,double test) ;
/* functional form: double (*fvec[])(double *x) */
int solnx(double *x,double *f,double (*fvec[])(),double *jm,
int n,double test) ;
/* functional form: double (*fvec[])(double *x) */
/* Curve Fitting and Least Squares */
void chcof(double *c,int m,double (*func)()) ;
/* functional form: double (*func)(double) */
void chpade(double *c,double *a,int m,double *b,int n) ;
double ftch(double x,double *a,int m,double *b,int n) ;
void cspl(double *x,double *y,double *z,int m,double tn) ;
void csplp(double *x,double *y,double *z,int m,double tn) ;
double csfit(double w,double *x,double *y,double *z,int m) ;
double tnsfit(double w,double *x,double *y,double *z,
int m,double tn) ;
double dcspl(double x,double *u,double *v,double *z,int m) ;
/* polynominal least squares functions use the Opol structure. */
void plsq(double *x,double *y,int n,Opol *c,double *ssq,int m) ;
double pplsq(double *x,double *y,int n,double *b,int m) ;
double evpsq(double x,Opol *c,int m) ;
double evpsqv(double x, Opol *c,int m,double *sig,double sqv) ;
void psqcf(double *pc,Opol *c,int m) ;
void psqvar(double *var,double s,Opol *c,int m) ;
/* QR transformation for linear least squares. */
double qrlsq(double *a,double *b,int m,int n,int *f) ;
double qrvar(double *v,int m,int n,double ssq) ;
/* singular value decomposition least squares. */
double lsqsv(double *x,int *pr,double *var,double *d,double *b,
double *v,int m,int n,double th) ;
int svdlsq(double *d,double *a,double *b,int m,double *v,int n) ;
int sv2lsq(double *d,double *a,double *b,int m,double *v,int n) ;
/* utility called by svdlsq and sv2lsq. */
int qrbdbv(double *d,double *e,double *b,double *v,int n) ;
/* nonlinear least squares */
double seqlsq(double *x,double *y,int n,double *par,double *var,
int m,double de,double (*func)(),int kf) ;
/* functional form: double (*func)(double x,double *par) */
double gnlsq(double *x,double *y,int n,double *par,
double *var,int m,double de,double (*func)()) ;
/* functional form: double (*func)(double x,double *par) */
double fitval(double x,double *s,double *par,double (*fun)(),
double *v,int n) ;
/* functional form: double (*func)(double x,double *par) */
void setfval(int i,int n) ;
/* Fourier Analysis */
void fft2(struct complex *ft,int m,int inv) ;
void fft2_d(struct complex *a,int m,int n,int f) ;
void fftgc(struct complex **pc,struct complex *ft,int n,
int *kk,int inv) ;
void fftgr(double *x,struct complex *ft,int n,int *kk,int inv) ;
void ftuns(struct complex **pt,int n) ;
int pfac(int n,int *kk,int fe) ;
void pshuf(Cpx **pa,Cpx **pb,int *kk,int n) ;
int pwspec(double *x,int n,int m) ;
void smoo(double *x,int n,int m) ;
/* Simulation Support */
double *autcor(double *x,int n,int lag) ;
int *hist(double *x,int n,double xmin,double xmax,
int kbin,double *bin) ;
unsigned int lran1() ;
void setlran1(unsigned int seed) ;
unsigned int lrand() ;
void setlrand(unsigned int seed) ;
int bran(int n) ;
void setbran(unsigned int seed) ;
int bran2(int n) ;
void setbran2(unsigned int seed) ;
double unfl() ;
void setunfl(unsigned int seed) ;
double unfl2() ;
void setunfl2(unsigned int seed) ;
double nrml() ;
void setnrml(unsigned int seed) ;
void norm(double *err) ;
void setnorm(unsigned int seed) ;
void norm2(double *err) ;
void setnorm2(unsigned int seed) ;
void sampl(void **s,int n,void **d,int m) ;
void shuffl(void **s,int n) ;
/* utility routines used for 2^31 - 1 modular arithmetic */
unsigned int lrana(unsigned int s) ;
unsigned int lranb(unsigned int s) ;
/* Sorts and Searches */
int batdel(char *kin,struct tnode *hd) ;
struct tnode *batins(char *kin,struct tnode *hd) ;
struct tnode *btsearch(char *kin,struct tnode *hd) ;
void btsort(struct tnode *hd,struct tnode **ar) ;
void prbtree(struct tnode *hd,int m) ;
int btdel(char *kin,struct tnode *hd) ;
struct tnode *btins(char *kin,struct tnode *hd) ;
struct tnode *tsearch(char *kin,struct tnode *hd) ;
void tsort(struct tnode *hd,struct tnode **ar) ;
void prtree(struct tnode *hd,int m) ;
int hashdel(char *kin,struct tabl *harr[],int mh) ;
struct tabl *hashins(char *kin,struct tabl *harr[],int mh) ;
struct tabl *hfind(char *kin,struct tabl *harr[],int mh) ;
int hval(char *key,int mh) ;
struct llst *msort(struct llst *st,int dim,int (*comp)()) ;
void qsrt(void *v,int i,int j,int (*comp)()) ;
void hsort(void *v,int n,int (*comp)()) ;
void ssort(void *v,int n,int (*comp)()) ;
/* comparison functions for sort routines. */
/* define the functional form of int (*comp)() */
int dubcmp(double *x,double *y) ;
int intcmp(int *x,int *y) ;
int unicmp(unsigned *x,unsigned *y) ;
/* the standard library function strcmp will also work
with these sorts */
/* Statistical Distributions */
double qnorm(double x) ;
double pctn(double pc) ;
double qgama(double x,double a) ;
double pctg(double pc,double a) ;
double qbeta(double x,double a,double b) ;
double pctb(double pc,double a,double b) ;
double qgnc(double x,double a,double d) ;
double pctgn(double pc,double a,double d) ;
double qbnc(double x,double a,double b,double d) ;
double pctbn(double pc,double a,double b,double d) ;
/* Special Functions */
/* elliptic integrals and functions */
double nome(double k,double *pk,double *pkp) ;
double amelp(double u,double k) ;
double theta(double u,int n) ;
void stheta(double k) ;
double felp(double an,double k,double *pk,double *pz,double *ph) ;
double gelp(double an,double k,double as,double bs,
double ds,double *pg,double *pf,double *pk) ;
double g2elp(double an,double bn,double k,double as,
double bs,double ds) ;
/* bessel functions */
double jbes(double v,double x) ;
double ibes(double v,double x) ;
double kbes(double v,double x) ;
double nbes(double v,double x) ;
double drbes(double x,double v,int f,double *p) ;
double rcbes() ;
void setrcb(double u,double y,int fl,int dr,double *pf,
double *ph) ;
/* spherical bessel functions */
double jspbes(int n,double x) ;
double kspbes(int n,double x) ;
double yspbes(int n,double x) ;
double drspbes(double x,int n,int f,double *p) ;
double rcspbs() ;
void setrcsb(int n,double y,int fl,int dr,double *pf,double *ph) ;
/* airy functions */
double airy(double x,int df) ;
double biry(double x,int df) ;
/* gamma and related functions */
double gaml(double x) ;
double psi(int m) ;
double psih(double v) ;
/* support routines for evaluation of elliptic integrals */
double gsng(double *pa,double *pb,double *pc,double b,double an) ;
double gsng2(double *pa,double *pb,double *pc,double b,
double an,double bn) ;
/* Complex Arithmetic */
struct complex cmul(struct complex s,struct complex t) ;
struct complex cdiv(struct complex s,struct complex t) ;
struct complex cadd(struct complex s,struct complex t) ;
struct complex csub(struct complex s,struct complex t) ;
struct complex crmu(double a,struct complex z) ;
struct complex cimu(double b,struct complex z) ;
struct complex ccng(struct complex z) ;
struct complex cdef(double r,double i) ;
double cabs(struct complex c) ;
double cnrm(struct complex z) ;
struct complex cexp(struct complex z) ;
struct complex clog(struct complex z) ;
struct complex csinh(struct complex z) ;
struct complex ccosh(struct complex z) ;
struct complex ctanh(struct complex z) ;
struct complex casinh(struct complex z) ;
struct complex cacosh(struct complex z) ;
struct complex catanh(struct complex z) ;
struct complex casin(struct complex z) ;
struct complex cacos(struct complex z) ;
struct complex catan(struct complex z) ;
struct complex csqrt(struct complex z) ;
struct complex csin(struct complex z) ;
struct complex ccos(struct complex z) ;
struct complex ctan(struct complex z) ;
/* Time Series */
double sarma(double er) ;
void setsim(int k) ;
double parma(double *x,double *e) ;
double evfmod(struct fmod y) ;
void setevf(int k) ;
double drfmod(struct fmod y,double *dr) ;
void setdrf(int k) ;
double seqtsf(struct fmod *x,int n,double *var,int kf) ;
double fixtsf(struct fmod *x,int n,double *var,double *cr) ;
double evmod(double y) ;
void setev(int k) ;
double drmod(double y,double *dr) ;
void setdr(int k) ;
double seqts(double *x,int n,double *var,int kf) ;
double fixts(double *x,int n,double *var,double *cr) ;
int resid(double *x,int n,int lag,double **pau,int nbin,
double xa,double xb,int **phs,int *cks) ;
int sany(double *x,int n,double *pm,double *cd,double *ci,
int nd,int ms,int lag) ;
double sdiff(double y,int nd,int k) ;
double sintg(double y,int nd,int k) ;
double xmean(double *x,int n) ;
/* Extended Precision Arithmetic */
/* XMATH must be defined to use these functions */
#ifdef XMATH
struct xpr xadd(struct xpr s,struct xpr t,int f) ;
struct xpr xmul(struct xpr s,struct xpr t) ;
struct xpr xdiv(struct xpr s,struct xpr t) ;
double xtodub(struct xpr s) ;
struct xpr dubtox(double y) ;
struct xpr inttox(int n) ;
int xprcmp(struct xpr *pa,struct xpr *pb) ;
struct xpr xneg(struct xpr s) ;
struct xpr xabs(struct xpr s) ;
int xex(struct xpr *ps) ;
int neg(struct xpr *ps) ;
struct xpr xfrex(struct xpr s,int *p) ;
struct xpr xfmod(struct xpr s,struct xpr t,int *p) ;
struct xpr xsqrt(struct xpr z) ;
struct xpr xexp(struct xpr z) ;
struct xpr xlog(struct xpr z) ;
struct xpr xpwr(struct xpr s,int n) ;
struct xpr xpr2(struct xpr s,int m) ;
struct xpr xtan(struct xpr z) ;
struct xpr xcos(struct xpr z) ;
struct xpr xsin(struct xpr z) ;
struct xpr xatan(struct xpr z) ;
struct xpr xasin(struct xpr z) ;
struct xpr xacos(struct xpr z) ;
struct xpr xtanh(struct xpr z) ;
struct xpr xsinh(struct xpr z) ;
struct xpr xcosh(struct xpr z) ;
struct xpr atox(char *s) ;
void prxpr(struct xpr u,int lim) ;
void xprint(struct xpr x) ;
/* special applications */
void xchcof(struct xpr *cf,int m,struct xpr (*xfunc)()) ;
/* functional form: xpr (*xfunc)(xpr *cf) */
struct xpr xevtch(struct xpr z,struct xpr *a,int m) ;
/* utility operations on extended precision numbers */
struct xpr sfmod(struct xpr s,int *p) ;
void lshift(int n,unsigned short *pm,int m) ;
void rshift(int n,unsigned short *pm,int m) ;
#endif
/* Utility Operations (on Bits) */
unsigned short bset(unsigned short x,unsigned short n) ;
int bget(unsigned short x,unsigned short n) ;
int bcnt(unsigned short x) ;
unsigned int lbset(unsigned int x,int n) ;
int lbget(unsigned int x,int n) ;
int lbcnt(unsigned int x) ;
void bitpc(unsigned char x) ;
void bitps(unsigned short x) ;
void bitpl(unsigned int x) ;
void bitpf(float x);
void bitpd(double x) ;
#ifdef XMATH
void bpatx(struct xpr x) ;
#endif
double pwr(double y,int n) ;
/*
special declarations required for shared library
*/
int np,nma,nar,nfc,ndif;
struct mcof *par,*pma,*pfc;
syntax highlighted by Code2HTML, v. 0.9.1