#include "header_cpp.h" //#include "mathfunc.h" #include #include "mathnr.h" #include "util.h" void nrerror(char error_text[]) { fprintf(stderr, "%s\n", error_text); exit(0); } double erfcc(double x) { double t,z,ans; z=fabs(x); t=1.0/(1.0+0.5*z); ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+ t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+ t*(-0.82215223+t*0.17087277))))))))); return x >= 0.0 ? ans : 2.0-ans; } /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ double gammln(double xx) { double x,y,tmp,ser; static double cof[6]={76.18009172947146,-86.50532032941677, 24.01409824083091,-1.231739572450155, 0.1208650973866179e-2,-0.5395239384953e-5}; int j; y=x=xx; tmp=x+5.5; tmp -= (x+0.5)*log(tmp); ser=1.000000000190015; for (j=0;j<=5;j++) ser += cof[j]/++y; return -tmp+log(2.5066282746310005*ser/x); } /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #include #define MAXIT 100 #define EPS 3.0e-7 #define FPMIN 1.0e-30 double betacf(double a, double b, double x) { void nrerror(char error_text[]); int m,m2; double aa,c,d,del,h,qab,qam,qap; qab=a+b; qap=a+1.0; qam=a-1.0; c=1.0; d=1.0-qab*x/qap; if (fabs(d) < FPMIN) d=FPMIN; d=1.0/d; h=d; for (m=1;m<=MAXIT;m++) { m2=2*m; aa=m*(b-m)*x/((qam+m2)*(a+m2)); d=1.0+aa*d; if (fabs(d) < FPMIN) d=FPMIN; c=1.0+aa/c; if (fabs(c) < FPMIN) c=FPMIN; d=1.0/d; h *= d*c; aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2)); d=1.0+aa*d; if (fabs(d) < FPMIN) d=FPMIN; c=1.0+aa/c; if (fabs(c) < FPMIN) c=FPMIN; d=1.0/d; del=d*c; h *= del; if (fabs(del-1.0) < EPS) break; } if (m > MAXIT) nrerror("a or b too big, or MAXIT too small in betacf"); return h; } #undef MAXIT #undef EPS #undef FPMIN /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ double betai(double a, double b, double x) { double betacf(double a, double b, double x); double gammln(double xx); void nrerror(char error_text[]); double bt; if (x < 0.0 || x > 1.0) nrerror("Bad x in routine betai"); if (x == 0.0 || x == 1.0) bt=0.0; else bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x)); if (x < (a+1.0)/(a+b+2.0)) return bt*betacf(a,b,x)/a; else return 1.0-bt*betacf(b,a,1.0-x)/b; } /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #define TINY 1.0e-20 void pearsn(double x[], double y[], unsigned long n, double *r, double *prob, double *z) { double betai(double a, double b, double x); double erfcc(double x); unsigned long j; double yt,xt,t,df; double syy=0.0,sxy=0.0,sxx=0.0,ay=0.0,ax=0.0; for (j=1;j<=n;j++) { ax += x[j]; ay += y[j]; } ax /= n; ay /= n; for (j=1;j<=n;j++) { xt=x[j]-ax; yt=y[j]-ay; sxx += xt*xt; syy += yt*yt; sxy += xt*yt; } *r=sxy/sqrt(sxx*syy); *z=0.5*log((1.0+(*r)+TINY)/(1.0-(*r)+TINY)); df=n-2; t=(*r)*sqrt(df/((1.0-(*r)+TINY)*(1.0+(*r)+TINY))); *prob=betai(0.5*df,0.5,df/(df+t*t)); } #undef TINY #define MBIG 1000000000 #define MSEED 161803398 #define MZ 0 #define FAC (1.0/MBIG) double ran3(long *idum) { static int inext,inextp; static long ma[56]; static int iff=0; long mj,mk; int i,ii,k; if (*idum < 0 || iff == 0) { iff=1; mj=MSEED-(*idum < 0 ? -*idum : *idum); mj %= MBIG; ma[55]=mj; mk=1; for (i=1;i<=54;i++) { ii=(21*i) % 55; ma[ii]=mk; mk=mj-mk; if (mk < MZ) mk += MBIG; mj=ma[ii]; } for (k=1;k<=4;k++) for (i=1;i<=55;i++) { ma[i] -= ma[1+(i+30) % 55]; if (ma[i] < MZ) ma[i] += MBIG; } inext=0; inextp=31; *idum=1; } if (++inext == 56) inext=1; if (++inextp == 56) inextp=1; mj=ma[inext]-ma[inextp]; if (mj < MZ) mj += MBIG; ma[inext]=mj; return mj*FAC; } #undef MBIG #undef MSEED #undef MZ #undef FAC /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ void moment(double data[], int n, double *ave, double *adev, double *sdev, double *var, double *skew, double *curt) { void nrerror(char error_text[]); int j; double ep=0.0,s,p; if (n <= 1) nrerror("n must be at least 2 in moment"); s=0.0; for (j=1;j<=n;j++) s += data[j]; *ave=s/n; *adev=(*var)=(*skew)=(*curt)=0.0; for (j=1;j<=n;j++) { *adev += fabs(s=data[j]-(*ave)); *var += (p=s*s); *skew += (p *= s); *curt += (p *= s); } *adev /= n; *var=(*var-ep*ep/n)/(n-1); *sdev=sqrt(*var); if (*var) { *skew /= (n*(*var)*(*sdev)); *curt=(*curt)/(n*(*var)*(*var))-3.0; } else nrerror("No skew/kurtosis when variance = 0 (in moment)"); } /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ void moment2(double data[], int n, double *ave, double *sdev) { void nrerror(char error_text[]); int j; double ep=0.0,s,p; double *var = new double [1]; *var = 0; if (n <= 1) nrerror("n must be at least 2 in moment"); s=0.0; for (j=1;j<=n;j++) s += data[j]; *ave=s/n; for (j=1;j<=n;j++) { s=data[j]-(*ave); *var += (p=s*s); } *var=(*var-ep*ep/n)/(n-1); *sdev=sqrt(*var); delete [] var; } /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #define NRANSI double pythag(double a, double b) { double absa,absb; absa=fabs(a); absb=fabs(b); if (absa > absb) return absa*sqrt(1.0+DSQR(absb/absa)); else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+DSQR(absa/absb))); } #undef NRANSI /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #define RADIX 2.0 void balanc(double **a, int n) { int last,j,i; double s,r,g,f,c,sqrdx; sqrdx=RADIX*RADIX; last=0; while (last == 0) { last=1; for (i=1;i<=n;i++) { r=c=0.0; for (j=1;j<=n;j++) if (j != i) { c += fabs(a[j][i]); r += fabs(a[i][j]); } if (c && r) { g=r/RADIX; f=1.0; s=c+r; while (cg) { f /= RADIX; c /= sqrdx; } if ((c+r)/f < 0.95*s) { last=0; g=1.0/f; for (j=1;j<=n;j++) a[i][j] *= g; for (j=1;j<=n;j++) a[j][i] *= f; } } } } } #undef RADIX /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #include #define SWAP(g,h) {y=(g);(g)=(h);(h)=y;} void elmhes(double **a, int n) { int m,j,i; double y,x; for (m=2;m fabs(x)) { x=a[j][m-1]; i=j; } } if (i != m) { for (j=m-1;j<=n;j++) SWAP(a[i][j],a[m][j]) for (j=1;j<=n;j++) SWAP(a[j][i],a[j][m]) } if (x) { for (i=m+1;i<=n;i++) { if ((y=a[i][m-1]) != 0.0) { y /= x; a[i][m-1]=y; for (j=m;j<=n;j++) a[i][j] -= y*a[m][j]; for (j=1;j<=n;j++) a[j][m] += y*a[j][i]; } } } } } #undef SWAP /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #include #define NRANSI void hqr(double **a, int n, double wr[], double wi[]) { int nn,m,l,k,j,its,i,mmin; double z,y,x,w,v,u,t,s,r,q,p,anorm; anorm=fabs(a[1][1]); for (i=2;i<=n;i++) for (j=(i-1);j<=n;j++) anorm += fabs(a[i][j]); nn=n; t=0.0; while (nn >= 1) { its=0; do { for (l=nn;l>=2;l--) { s=fabs(a[l-1][l-1])+fabs(a[l][l]); if (s == 0.0) s=anorm; if ((double)(fabs(a[l][l-1]) + s) == s) break; } x=a[nn][nn]; if (l == nn) { wr[nn]=x+t; wi[nn--]=0.0; } else { y=a[nn-1][nn-1]; w=a[nn][nn-1]*a[nn-1][nn]; if (l == (nn-1)) { p=0.5*(y-x); q=p*p+w; z=sqrt(fabs(q)); x += t; if (q >= 0.0) { z=p+SIGN(z,p); wr[nn-1]=wr[nn]=x+z; if (z) wr[nn]=x-w/z; wi[nn-1]=wi[nn]=0.0; } else { wr[nn-1]=wr[nn]=x+p; wi[nn-1]= -(wi[nn]=z); } nn -= 2; } else { if (its == 30) nrerror("Too many iterations in hqr"); if (its == 10 || its == 20) { t += x; for (i=1;i<=nn;i++) a[i][i] -= x; s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]); y=x=0.75*s; w = -0.4375*s*s; } ++its; for (m=(nn-2);m>=l;m--) { z=a[m][m]; r=x-z; s=y-z; p=(r*s-w)/a[m+1][m]+a[m][m+1]; q=a[m+1][m+1]-z-r-s; r=a[m+2][m+1]; s=fabs(p)+fabs(q)+fabs(r); p /= s; q /= s; r /= s; if (m == l) break; u=fabs(a[m][m-1])*(fabs(q)+fabs(r)); v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1])); if ((double)(u+v) == v) break; } for (i=m+2;i<=nn;i++) { a[i][i-2]=0.0; if (i != (m+2)) a[i][i-3]=0.0; } for (k=m;k<=nn-1;k++) { if (k != m) { p=a[k][k-1]; q=a[k+1][k-1]; r=0.0; if (k != (nn-1)) r=a[k+2][k-1]; if ((x=fabs(p)+fabs(q)+fabs(r)) != 0.0) { p /= x; q /= x; r /= x; } } if ((s=SIGN(sqrt(p*p+q*q+r*r),p)) != 0.0) { if (k == m) { if (l != m) a[k][k-1] = -a[k][k-1]; } else a[k][k-1] = -s*x; p += s; x=p/s; y=q/s; z=r/s; q /= p; r /= p; for (j=k;j<=nn;j++) { p=a[k][j]+q*a[k+1][j]; if (k != (nn-1)) { p += r*a[k+2][j]; a[k+2][j] -= p*z; } a[k+1][j] -= p*y; a[k][j] -= p*x; } mmin = nn #define NRANSI void svdcmp(double **a, int m, int n, double w[], double **v) { double pythag(double a, double b); int flag,i,its,j,jj,k,l,nm; double anorm,c,f,g,h,s,scale,x,y,z,*rv1; rv1=dvector(n); g=scale=anorm=0.0; for (i=1;i<=n;i++) { 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; 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; g=s=scale=0.0; if (i <= m && i != n) { 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; 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=DMAX(anorm,(fabs(w[i])+fabs(rv1[i]))); } for (i=n;i>=1;i--) { if (i < n) { 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; } for (i=IMIN(m,n);i>=1;i--) { l=i+1; g=w[i]; for (j=l;j<=n;j++) a[i][j]=0.0; if (g) { g=1.0/g; 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]; } for (k=n;k>=1;k--) { for (its=1;its<=30;its++) { flag=1; for (l=k;l>=1;l--) { nm=l-1; if ((double)(fabs(rv1[l])+anorm) == anorm) { flag=0; break; } if ((double)(fabs(w[nm])+anorm) == anorm) break; } if (flag) { c=0.0; s=1.0; for (i=l;i<=k;i++) { f=s*rv1[i]; rv1[i]=c*rv1[i]; if ((double)(fabs(f)+anorm) == anorm) break; g=w[i]; h=pythag(f,g); w[i]=h; h=1.0/h; c=g*h; s = -f*h; for (j=1;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) { if (z < 0.0) { w[k] = -z; for (j=1;j<=n;j++) v[j][k] = -v[j][k]; } break; } if (its == 30) nrerror("no convergence in 30 svdcmp iterations"); 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; 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 *= c; for (jj=1;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=1;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; } } delete [] rv1; //free_vector(rv1,1,n); } #undef NRANSI /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #include #define NRANSI #define TINY 1.0e-20; void ludcmp(double **a, int n, int *indx, double *d) { int i,imax,j,k; double big,dum,sum,temp; double *vv; //vv=vector(1,n); vv = dvector(n); *d=1.0; for (i=1;i<=n;i++) { big=0.0; for (j=1;j<=n;j++) if ((temp=fabs(a[i][j])) > big) big=temp; if (big == 0.0) nrerror("Singular matrix in routine ludcmp"); vv[i]=1.0/big; } for (j=1;j<=n;j++) { for (i=1;i= big) { big=dum; imax=i; } } if (j != imax) { for (k=1;k<=n;k++) { dum=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=dum; } *d = -(*d); vv[imax]=vv[j]; } indx[j]=imax; if (a[j][j] == 0.0) a[j][j]=TINY; if (j != n) { dum=1.0/(a[j][j]); for (i=j+1;i<=n;i++) a[i][j] *= dum; } } //free_vector(vv,1,n); delete [] vv; } #undef TINY #undef NRANSI /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ void lubksb(double **a, int n, int *indx, double b[]) { int i,ii=0,ip,j; double sum; for (i=1;i<=n;i++) { ip=indx[i]; sum=b[ip]; b[ip]=b[i]; if (ii) for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j]; else if (sum) ii=i; b[i]=sum; } for (i=n;i>=1;i--) { sum=b[i]; for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j]; b[i]=sum/a[i][i]; } } /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #define NRANSI //#include "nrutil.h" #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp; #define M 7 #define NSTACK 50 void sort(unsigned long n, double arr[]) { unsigned long i,ir=n,j,k,l=1; int jstack=0,*istack; double a,temp; //istack=ivector(1,NSTACK); istack = ivector(NSTACK); for (;;) { if (ir-l < M) { for (j=l+1;j<=ir;j++) { a=arr[j]; for (i=j-1;i>=1;i--) { if (arr[i] <= a) break; arr[i+1]=arr[i]; } arr[i+1]=a; } if (jstack == 0) break; ir=istack[jstack--]; l=istack[jstack--]; } else { k=(l+ir) >> 1; SWAP(arr[k],arr[l+1]) if (arr[l+1] > arr[ir]) { SWAP(arr[l+1],arr[ir]) } if (arr[l] > arr[ir]) { SWAP(arr[l],arr[ir]) } if (arr[l+1] > arr[l]) { SWAP(arr[l+1],arr[l]) } i=l+1; j=ir; a=arr[l]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAP(arr[i],arr[j]); } arr[l]=arr[j]; arr[j]=a; jstack += 2; if (jstack > NSTACK) nrerror("NSTACK too small in sort."); if (ir-i+1 >= j-l) { istack[jstack]=ir; istack[jstack-1]=i; ir=j-1; } else { istack[jstack]=j-1; istack[jstack-1]=l; l=i; } } } //free_ivector(istack,1,NSTACK); delete [] istack; } #undef M #undef NSTACK #undef SWAP #undef NRANSI /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */ #define NRANSI //#include "nrutil.h" #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp; #define M 7 #define NSTACK 50 /* template void sort(unsigned long n, kind arr[]) { unsigned long i,ir=n,j,k,l=1; int jstack=0,*istack; kind a,temp; //istack=ivector(1,NSTACK); istack = ivector(NSTACK); for (;;) { if (ir-l < M) { for (j=l+1;j<=ir;j++) { a=arr[j]; for (i=j-1;i>=1;i--) { if (arr[i] <= a) break; arr[i+1]=arr[i]; } arr[i+1]=a; } if (jstack == 0) break; ir=istack[jstack--]; l=istack[jstack--]; } else { k=(l+ir) >> 1; SWAP(arr[k],arr[l+1]) if (arr[l+1] > arr[ir]) { SWAP(arr[l+1],arr[ir]) } if (arr[l] > arr[ir]) { SWAP(arr[l],arr[ir]) } if (arr[l+1] > arr[l]) { SWAP(arr[l+1],arr[l]) } i=l+1; j=ir; a=arr[l]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAP(arr[i],arr[j]); } arr[l]=arr[j]; arr[j]=a; jstack += 2; if (jstack > NSTACK) nrerror("NSTACK too small in sort."); if (ir-i+1 >= j-l) { istack[jstack]=ir; istack[jstack-1]=i; ir=j-1; } else { istack[jstack]=j-1; istack[jstack-1]=l; l=i; } } } //free_ivector(istack,1,NSTACK); delete [] istack; } #undef M #undef NSTACK #undef SWAP #undef NRANSI */ /* (C) Copr. 1986-92 Numerical Recipes Software 21].,t45316. */