开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 生物技术 > 查看源码
Recipes.cpp
资源名称：ClusterTreeViewCode.zip　[点击查看]
上传用户：szb0815
上传日期：2007-06-13
资源大小：338k
文件大小：35k
源码类别：

生物技术

开发平台：

C++ Builder

Recipes.cpp：源码内容
  1. //---------------------------------------------------------------------------
  2. #pragma hdrstop
  4. #include "Recipes.h"
  6. float *circx,*circy,*circw;
  7. int circn;
  9. //---------------------------------------------------------------------------
  10. #pragma package(smart_init)
  11. #include <math.h>
  12. #define NRANSI #define EPS 1.0e-7 void medfit(float x[], float y[], int ndata, float *a, float *b, float *abdev) { float rofunc(float b); int j; float bb,b1,b2,del,f,f1,f2,sigb,temp; float sx=0.0,sy=0.0,sxy=0.0,sxx=0.0,chisq=0.0; ndatat=ndata; xt=x; yt=y; for (j=1;j<=ndata;j++) { sx += x[j]; sy += y[j]; sxy += x[j]*y[j]; sxx += x[j]*x[j]; } del=ndata*sxx-sx*sx; aa=(sxx*sy-sx*sxy)/del; bb=(ndata*sxy-sx*sy)/del; for (j=1;j<=ndata;j++) chisq += (temp=y[j]-(aa+bb*x[j]),temp*temp); sigb=sqrt(chisq/del); b1=bb; f1=rofunc(b1); b2=bb+SIGN(3.0*sigb,f1); f2=rofunc(b2); if (b2 == b1) { *a=aa; *b=bb; *abdev=abdevt/ndata; return; } while (f1*f2 > 0.0) { bb=b2+1.6*(b2-b1); b1=b2; f1=f2; b2=bb; f2=rofunc(b2); } sigb=0.01*sigb; while (fabs(b2-b1) > sigb) { bb=b1+0.5*(b2-b1); if (bb == b1 || bb == b2) break; f=rofunc(bb); if (f*f1 >= 0.0) { f1=f; b1=bb; } else { f2=f; b2=bb; } } *a=aa; *b=bb; *abdev=abdevt/ndata; }
  13. float rofunc(float b) { //float select(unsigned long k, unsigned long n, float arr[]); int j; float *arr,d,sum=0.0; arr = new float[ndatat+1]; for (j=1;j<=ndatat;j++)         {             arr[j]=yt[j]-b*xt[j];         } if (ndatat & 1) { aa=select((ndatat+1)>>1,ndatat,arr); } else { j=ndatat >> 1; aa=0.5*(select(j,ndatat,arr)+select(j+1,ndatat,arr)); } abdevt=0.0; for (j=1;j<=ndatat;j++) { d=yt[j]-(b*xt[j]+aa); abdevt += fabs(d); if (yt[j] != 0.0) d /= fabs(yt[j]); if (fabs(d) > EPS) sum += (d >= 0.0 ? xt[j] : -xt[j]); }         delete arr; return sum; } float select(unsigned long k, unsigned long n, float arr[]) { unsigned long i,ir,j,l,mid; float a,temp; l=1; ir=n; for (;;) { if (ir <= l+1) { if (ir == l+1 && arr[ir] < arr[l]) { SWAP(arr[l],arr[ir]) } return arr[k]; } else { mid=(l+ir) >> 1; SWAP(arr[mid],arr[l+1]) if (arr[l] > arr[ir]) { SWAP(arr[l],arr[ir]) } if (arr[l+1] > arr[ir]) { SWAP(arr[l+1],arr[ir]) } if (arr[l] > arr[l+1]) { SWAP(arr[l],arr[l+1]) } i=l+1; j=ir; a=arr[l+1]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAP(arr[i],arr[j]) } arr[l+1]=arr[j]; arr[j]=a; if (j >= k) ir=j-1; if (j <= k) l=i; } } } #undef SWAP #undef EPS #undef NRANSI
  16. float fcircfit(float p[])
  17. {
  18.         int i;
  19.         float val = 0;
  20.         float val1,val2;
  22.         for (i=1;i<=circn;i++)
  23.         {
  24.                 val1 = sqrt(pow(circx[i]-p[1],2)+pow(circy[i]-p[2],2));
  25.                 val2 = val1 - p[3];
  26.                 val += circw[i] * pow(val2,2);
  27.         }
  29.         return val;
  30. }
  32. void fcircderiv(float p[], float xi[])
  33. {
  34.         int i;
  35.         float val1,val2;
  37.         xi[1] = 0;
  38.         xi[2] = 0;
  39.         xi[3] = 0;
  41.         for (i=1;i<=circn;i++)
  42.         {
  43.                 val1 = sqrt(pow(circx[i]-p[1],2)+pow(circy[i]-p[2],2));
  44.                 val2 = val1 - p[3];
  46.                 if (val1 > 0)
  47.                 {
  49.                 xi[1] += circw[i] * 2*p[1]*val2*(circx[i] - p[1])/val1;
  50.                 xi[2] += circw[i] * 2*p[2]*val2*(circy[i] - p[2])/val1;
  51.                 xi[3] += circw[i] * 2 * val2;
  52.                 }
  53.         }
  54. }
  56. #define ITMAX 200
  57. #define EPS 1.0e-10 #define FREEALL free_vector(xi,1,n);free_vector(h,1,n);free_vector(g,1,n); void frprmn(float p[], int n, float ftol, int *iter, float *fret, float (*func)(float []), void (*dfunc)(float [], float [])) { void linmin(float p[], float xi[], int n, float *fret, float (*func)(float [])); int j,its; float gg,gam,fp,dgg; float *g,*h,*xi; g=vector(1,n); h=vector(1,n); xi=vector(1,n); fp=(*func)(p); (*dfunc)(p,xi); for (j=1;j<=n;j++) { g[j] = -xi[j]; xi[j]=h[j]=g[j]; } for (its=1;its<=ITMAX;its++) { *iter=its; linmin(p,xi,n,fret,func); if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) { FREEALL return; } fp=(*func)(p); (*dfunc)(p,xi); dgg=gg=0.0; for (j=1;j<=n;j++) { gg += g[j]*g[j]; dgg += (xi[j]+g[j])*xi[j]; } if (gg == 0.0) { FREEALL return; } gam=dgg/gg; for (j=1;j<=n;j++) { g[j] = -xi[j]; xi[j]=h[j]=g[j]+gam*h[j]; } } nrerror("Too many iterations in frprmn"); } #undef ITMAX #undef EPS #undef FREEALL #define TOL 2.0e-4 void linmin(float p[], float xi[], int n, float *fret, float (*func)(float [])) { float brent(float ax, float bx, float cx, float (*f)(float), float tol, float *xmin); float f1dim(float x); void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb, float *fc, float (*func)(float)); int j; float xx,xmin,fx,fb,fa,bx,ax; ncom=n; pcom=vector(1,n); xicom=vector(1,n); nrfunc=func; for (j=1;j<=n;j++) { pcom[j]=p[j]; xicom[j]=xi[j]; } ax=0.0; xx=1.0; mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim); *fret=brent(ax,xx,bx,f1dim,TOL,&xmin); for (j=1;j<=n;j++) { xi[j] *= xmin; p[j] += xi[j]; } free_vector(xicom,1,n); free_vector(pcom,1,n); }
  58. #undef TOL
  59. #define ITMAX 100
  60. #define CGOLD 0.3819660 #define ZEPS 1.0e-10 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); float brent(float ax, float bx, float cx, float (*f)(float), float tol, float *xmin) { int iter; float a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; float e=0.0; a=(ax < cx ? ax : cx); b=(ax > cx ? ax : cx); x=w=v=bx; fw=fv=fx=(*f)(x); for (iter=1;iter<=ITMAX;iter++) { xm=0.5*(a+b); tol2=2.0*(tol1=tol*fabs(x)+ZEPS); if (fabs(x-xm) <= (tol2-0.5*(b-a))) { *xmin=x; return fx; } if (fabs(e) > tol1) { r=(x-w)*(fx-fv); q=(x-v)*(fx-fw); p=(x-v)*q-(x-w)*r; q=2.0*(q-r); if (q > 0.0) p = -p; q=fabs(q); etemp=e; e=d; if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) d=CGOLD*(e=(x >= xm ? a-x : b-x)); else { d=p/q; u=x+d; if (u-a < tol2 || b-u < tol2) d=SIGN(tol1,xm-x); } } else { d=CGOLD*(e=(x >= xm ? a-x : b-x)); } u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d)); fu=(*f)(u); if (fu <= fx) { if (u >= x) a=x; else b=x; SHFT(v,w,x,u) SHFT(fv,fw,fx,fu) } else { if (u < x) a=u; else b=u; if (fu <= fw || w == x) { v=w; w=u; fv=fw; fw=fu; } else if (fu <= fv || v == x || v == w) { v=u; fv=fu; } } } nrerror("Too many iterations in brent"); }
  62. #undef ITMAX
  63. #undef CGOLD
  64. #undef ZEPS
  65. #undef SHFT
  67. #define GOLD 1.618034
  68. #define GLIMIT 100.0 #define TINY 1.0e-20 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); void mnbrak(float *ax, float *bx, float *cx, float *fa, float *fb, float *fc, float (*func)(float)) { float ulim,u,r,q,fu,dum; *fa=(*func)(*ax); *fb=(*func)(*bx); if (*fb > *fa) { SHFT(dum,*ax,*bx,dum) SHFT(dum,*fb,*fa,dum) } *cx=(*bx)+GOLD*(*bx-*ax); *fc=(*func)(*cx); while (*fb > *fc) { r=(*bx-*ax)*(*fb-*fc); q=(*bx-*cx)*(*fb-*fa); u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); ulim=(*bx)+GLIMIT*(*cx-*bx); if ((*bx-u)*(u-*cx) > 0.0) { fu=(*func)(u); if (fu < *fc) { *ax=(*bx); *bx=u; *fa=(*fb); *fb=fu; return; } else if (fu > *fb) { *cx=u; *fc=fu; return; } u=(*cx)+GOLD*(*cx-*bx); fu=(*func)(u); } else if ((*cx-u)*(u-ulim) > 0.0) { fu=(*func)(u); if (fu < *fc) { SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) SHFT(*fb,*fc,fu,(*func)(u)) } } else if ((u-ulim)*(ulim-*cx) >= 0.0) { u=ulim; fu=(*func)(u); } else { u=(*cx)+GOLD*(*cx-*bx); fu=(*func)(u); } SHFT(*ax,*bx,*cx,u) SHFT(*fa,*fb,*fc,fu) } } #undef GOLD #undef GLIMIT #undef TINY #undef SHFT
  70. float f1dim(float x)
  71. { int j; float f,*xt; xt=vector(1,ncom); for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j]; f=(*nrfunc)(xt); free_vector(xt,1,ncom); return f; }
  73. void kstwo(float data1[], unsigned long n1, float data2[], unsigned long n2,
  74. float *d, float *prob) { float probks(float alam); void sort(unsigned long n, float arr[]); unsigned long j1=1,j2=1; float d1,d2,dt,en1,en2,en,fn1=0.0,fn2=0.0; sort(n1,data1); sort(n2,data2); en1=n1; en2=n2; *d=0.0; while (j1 <= n1 && j2 <= n2) { if ((d1=data1[j1]) <= (d2=data2[j2])) fn1=j1++/en1; if (d2 <= d1) fn2=j2++/en2; if ((dt=fabs(fn2-fn1)) > *d) *d=dt; } en=sqrt(en1*en2/(en1+en2)); *prob=probks((en+0.12+0.11/en)*(*d)); }
  76. #define EPS1 0.001
  77. #define EPS2 1.0e-8 float probks(float alam) { int j; float a2,fac=2.0,sum=0.0,term,termbf=0.0; a2 = -2.0*alam*alam; for (j=1;j<=100;j++) { term=fac*exp(a2*j*j); sum += term; if (fabs(term) <= EPS1*termbf || fabs(term) <= EPS2*sum) return sum; fac = -fac; termbf=fabs(term); } return 1.0; } #undef EPS1 #undef EPS2
  79. #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
  80. #define M 7 #define NSTACK 50 void sort(unsigned long n, float arr[]) { unsigned long i,ir=n,j,k,l=1; int jstack=0,*istack; float a,temp; istack=ivector(1,NSTACK); for (;;) { if (ir-l < M) { for (j=l+1;j<=ir;j++) { a=arr[j]; for (i=j-1;i>=l;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] > arr[ir]) { SWAP(arr[l],arr[ir]) } if (arr[l+1] > arr[ir]) { SWAP(arr[l+1],arr[ir]) } if (arr[l] > arr[l+1]) { SWAP(arr[l],arr[l+1]) } i=l+1; j=ir; a=arr[l+1]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAP(arr[i],arr[j]); } arr[l+1]=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); } #undef M #undef NSTACK #undef SWAP #undef NRANSI
  82. #define NRANSI
  83. #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp; #define M 7 #define NSTACK 50
  85. void sort2fi(unsigned long n, float arr[], int brr[])
  86. { unsigned long i,ir=n,j,k,l=1; int *istack,jstack=0; float a,temp;         int b; istack=ivector(1,NSTACK); for (;;) { if (ir-l < M) { for (j=l+1;j<=ir;j++) { a=arr[j]; b=brr[j]; for (i=j-1;i>=l;i--) { if (arr[i] <= a) break; arr[i+1]=arr[i]; brr[i+1]=brr[i]; } arr[i+1]=a; brr[i+1]=b; } if (!jstack) { free_ivector(istack,1,NSTACK); return; } ir=istack[jstack]; l=istack[jstack-1]; jstack -= 2; } else { k=(l+ir) >> 1; SWAP(arr[k],arr[l+1]) SWAP(brr[k],brr[l+1]) if (arr[l] > arr[ir]) { SWAP(arr[l],arr[ir]) SWAP(brr[l],brr[ir]) } if (arr[l+1] > arr[ir]) { SWAP(arr[l+1],arr[ir]) SWAP(brr[l+1],brr[ir]) } if (arr[l] > arr[l+1]) { SWAP(arr[l],arr[l+1]) SWAP(brr[l],brr[l+1]) } i=l+1; j=ir; a=arr[l+1]; b=brr[l+1]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAP(arr[i],arr[j]) SWAP(brr[i],brr[j]) } arr[l+1]=arr[j]; arr[j]=a; brr[l+1]=brr[j]; brr[j]=b; jstack += 2; if (jstack > NSTACK) nrerror("NSTACK too small in sort2."); 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; } } } } #undef M #undef NSTACK #undef SWAP #undef NRANSI
  88. #define NRANSI
  89. #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp; #define M 7 #define NSTACK 50
  91. void sort3fff(unsigned long n, float arr[], float brr[], float crr[])
  92. { unsigned long i,ir=n,j,k,l=1; int *istack,jstack=0; float a,temp;         float b,c; istack=ivector(1,NSTACK); for (;;) { if (ir-l < M) { for (j=l+1;j<=ir;j++) { a=arr[j]; b=brr[j];                                 c=crr[j]; for (i=j-1;i>=l;i--) { if (arr[i] <= a) break; arr[i+1]=arr[i]; brr[i+1]=brr[i];                                         crr[i+1]=crr[i]; } arr[i+1]=a; brr[i+1]=b;                                 crr[i+1]=c; } if (!jstack) { free_ivector(istack,1,NSTACK); return; } ir=istack[jstack]; l=istack[jstack-1]; jstack -= 2; } else { k=(l+ir) >> 1; SWAP(arr[k],arr[l+1]) SWAP(brr[k],brr[l+1])                         SWAP(crr[k],crr[l+1]) if (arr[l] > arr[ir]) { SWAP(arr[l],arr[ir]) SWAP(brr[l],brr[ir])                                 SWAP(crr[l],crr[ir]) } if (arr[l+1] > arr[ir]) { SWAP(arr[l+1],arr[ir]) SWAP(brr[l+1],brr[ir]) SWAP(crr[l+1],crr[ir]) } if (arr[l] > arr[l+1]) { SWAP(arr[l],arr[l+1]) SWAP(brr[l],brr[l+1]) SWAP(crr[l],crr[l+1]) } i=l+1; j=ir; a=arr[l+1]; b=brr[l+1];                         c=crr[l+1]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAP(arr[i],arr[j]) SWAP(brr[i],brr[j]) SWAP(crr[i],crr[j]) } arr[l+1]=arr[j]; arr[j]=a; brr[l+1]=brr[j]; brr[j]=b; crr[l+1]=crr[j]; crr[j]=c; jstack += 2; if (jstack > NSTACK) nrerror("NSTACK too small in sort2."); 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; } } } } #undef M #undef NSTACK #undef SWAP #undef NRANSI
  94. #include <math.h>
  95. #define TINY 1.0e-20 void pearsn(float x[], float y[], unsigned long n, float *r, float *prob, float *z) { float betai(float a, float b, float x); float erfcc(float x); unsigned long j; float yt,xt,t,df; float 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)+TINY); *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
  97. float betai(float a, float b, float x)
  98. { float betacf(float a, float b, float x); float gammln(float xx); void nrerror(char error_text[]); float 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; }
  100. float gammln(float 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); }
  102. #define MAXIT 100
  103. #define EPS 3.0e-7 #define FPMIN 1.0e-30 float betacf(float a, float b, float x) { void nrerror(char error_text[]); int m,m2; float 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
  104. //---------------------------------------------------------------------------
  105. #include <math.h>
  106. #define NRANSI
  107. #include "nrutil.h"
  109. void spear(float data1[], float data2[], unsigned long n, float *d, float *zd,
  110. float *probd, float *rs, float *probrs)
  111. {
  112. float betai(float a, float b, float x);
  113. void crank(unsigned long n, float w[], float *s);
  114. float erfcc(float x);
  115. void sort2(unsigned long n, float arr[], float brr[]);
  116. unsigned long j;
  117. float vard,t,sg,sf,fac,en3n,en,df,aved,*wksp1,*wksp2;
  119. wksp1=vector(1,n);
  120. wksp2=vector(1,n);
  121. for (j=1;j<=n;j++) {
  122. wksp1[j]=data1[j];
  123. wksp2[j]=data2[j];
  124. }
  125. sort2(n,wksp1,wksp2);
  126. crank(n,wksp1,&sf);
  127. sort2(n,wksp2,wksp1);
  128. crank(n,wksp2,&sg);
  129. *d=0.0;
  130. for (j=1;j<=n;j++)
  131. *d += SQR(wksp1[j]-wksp2[j]);
  132. en=n;
  133. en3n=en*en*en-en;
  134. aved=en3n/6.0-(sf+sg)/12.0;
  135. fac=(1.0-sf/en3n)*(1.0-sg/en3n);
  136. vard=((en-1.0)*en*en*SQR(en+1.0)/36.0)*fac;
  137. *zd=(*d-aved)/sqrt(vard);
  138. *probd=erfcc(fabs(*zd)/1.4142136);
  139. *rs=(1.0-(6.0/en3n)*(*d+(sf+sg)/12.0))/sqrt(fac);
  140. fac=(*rs+1.0)*(1.0-(*rs));
  141. if (fac > 0.0) {
  142. t=(*rs)*sqrt((en-2.0)/fac);
  143. df=en-2.0;
  144. *probrs=betai(0.5*df,0.5,df/(df+t*t));
  145. } else
  146. *probrs=0.0;
  147. free_vector(wksp2,1,n);
  148. free_vector(wksp1,1,n);
  149. }
  150. #undef NRANSI
  151. #include <math.h>
  153. float gammln(float xx)
  154. {
  155. double x,y,tmp,ser;
  156. static double cof[6]={76.18009172947146,-86.50532032941677,
  157. 24.01409824083091,-1.231739572450155,
  158. 0.1208650973866179e-2,-0.5395239384953e-5};
  159. int j;
  161. y=x=xx;
  162. tmp=x+5.5;
  163. tmp -= (x+0.5)*log(tmp);
  164. ser=1.000000000190015;
  165. for (j=0;j<=5;j++) ser += cof[j]/++y;
  166. return -tmp+log(2.5066282746310005*ser/x);
  167. }
  169. #define NR_END 1
  170. #define FREE_ARG char*
  172. void nrerror(char error_text[])
  173. /* Numerical Recipes standard error handler */
  174. {
  175. fprintf(stderr,"Numerical Recipes run-time error...n");
  176. fprintf(stderr,"%sn",error_text);
  177. fprintf(stderr,"...now exiting to system...n");
  178. exit(1);
  179. }
  181. float *vector(long nl, long nh)
  182. /* allocate a float vector with subscript range v[nl..nh] */
  183. {
  184. float *v;
  186. v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
  187. if (!v) nrerror("allocation failure in vector()");
  188. return v-nl+NR_END;
  189. }
  191. int *ivector(long nl, long nh)
  192. /* allocate an int vector with subscript range v[nl..nh] */
  193. {
  194. int *v;
  196. v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
  197. if (!v) nrerror("allocation failure in ivector()");
  198. return v-nl+NR_END;
  199. }
  201. unsigned char *cvector(long nl, long nh)
  202. /* allocate an unsigned char vector with subscript range v[nl..nh] */
  203. {
  204. unsigned char *v;
  206. v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
  207. if (!v) nrerror("allocation failure in cvector()");
  208. return v-nl+NR_END;
  209. }
  211. unsigned long *lvector(long nl, long nh)
  212. /* allocate an unsigned long vector with subscript range v[nl..nh] */
  213. {
  214. unsigned long *v;
  216. v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
  217. if (!v) nrerror("allocation failure in lvector()");
  218. return v-nl+NR_END;
  219. }
  221. double *dvector(long nl, long nh)
  222. /* allocate a double vector with subscript range v[nl..nh] */
  223. {
  224. double *v;
  226. v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
  227. if (!v) nrerror("allocation failure in dvector()");
  228. return v-nl+NR_END;
  229. }
  231. float **matrix(long nrl, long nrh, long ncl, long nch)
  232. /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
  233. {
  234. long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  235. float **m;
  237. /* allocate pointers to rows */
  238. m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
  239. if (!m) nrerror("allocation failure 1 in matrix()");
  240. m += NR_END;
  241. m -= nrl;
  243. /* allocate rows and set pointers to them */
  244. m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
  245. if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
  246. m[nrl] += NR_END;
  247. m[nrl] -= ncl;
  249. for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  251. /* return pointer to array of pointers to rows */
  252. return m;
  253. }
  255. double **dmatrix(long nrl, long nrh, long ncl, long nch)
  256. /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
  257. {
  258. long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  259. double **m;
  261. /* allocate pointers to rows */
  262. m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
  263. if (!m) nrerror("allocation failure 1 in matrix()");
  264. m += NR_END;
  265. m -= nrl;
  267. /* allocate rows and set pointers to them */
  268. m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
  269. if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
  270. m[nrl] += NR_END;
  271. m[nrl] -= ncl;
  273. for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  275. /* return pointer to array of pointers to rows */
  276. return m;
  277. }
  279. int **imatrix(long nrl, long nrh, long ncl, long nch)
  280. /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
  281. {
  282. long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
  283. int **m;
  285. /* allocate pointers to rows */
  286. m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
  287. if (!m) nrerror("allocation failure 1 in matrix()");
  288. m += NR_END;
  289. m -= nrl;
  292. /* allocate rows and set pointers to them */
  293. m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
  294. if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
  295. m[nrl] += NR_END;
  296. m[nrl] -= ncl;
  298. for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
  300. /* return pointer to array of pointers to rows */
  301. return m;
  302. }
  304. float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
  305. long newrl, long newcl)
  306. /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
  307. {
  308. long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
  309. float **m;
  311. /* allocate array of pointers to rows */
  312. m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
  313. if (!m) nrerror("allocation failure in submatrix()");
  314. m += NR_END;
  315. m -= newrl;
  317. /* set pointers to rows */
  318. for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
  320. /* return pointer to array of pointers to rows */
  321. return m;
  322. }
  324. float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
  325. /* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
  326. declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
  327. and ncol=nch-ncl+1. The routine should be called with the address
  328. &a[0][0] as the first argument. */
  329. {
  330. long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
  331. float **m;
  333. /* allocate pointers to rows */
  334. m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
  335. if (!m) nrerror("allocation failure in convert_matrix()");
  336. m += NR_END;
  337. m -= nrl;
  339. /* set pointers to rows */
  340. m[nrl]=a-ncl;
  341. for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
  342. /* return pointer to array of pointers to rows */
  343. return m;
  344. }
  346. float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
  347. /* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
  348. {
  349. long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
  350. float ***t;
  352. /* allocate pointers to pointers to rows */
  353. t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
  354. if (!t) nrerror("allocation failure 1 in f3tensor()");
  355. t += NR_END;
  356. t -= nrl;
  358. /* allocate pointers to rows and set pointers to them */
  359. t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
  360. if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
  361. t[nrl] += NR_END;
  362. t[nrl] -= ncl;
  364. /* allocate rows and set pointers to them */
  365. t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
  366. if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
  367. t[nrl][ncl] += NR_END;
  368. t[nrl][ncl] -= ndl;
  370. for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
  371. for(i=nrl+1;i<=nrh;i++) {
  372. t[i]=t[i-1]+ncol;
  373. t[i][ncl]=t[i-1][ncl]+ncol*ndep;
  374. for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
  375. }
  377. /* return pointer to array of pointers to rows */
  378. return t;
  379. }
  381. void free_vector(float *v, long nl, long nh)
  382. /* free a float vector allocated with vector() */
  383. {
  384. free((FREE_ARG) (v+nl-NR_END));
  385. }
  387. void free_ivector(int *v, long nl, long nh)
  388. /* free an int vector allocated with ivector() */
  389. {
  390. free((FREE_ARG) (v+nl-NR_END));
  391. }
  393. void free_cvector(unsigned char *v, long nl, long nh)
  394. /* free an unsigned char vector allocated with cvector() */
  395. {
  396. free((FREE_ARG) (v+nl-NR_END));
  397. }
  399. void free_lvector(unsigned long *v, long nl, long nh)
  400. /* free an unsigned long vector allocated with lvector() */
  401. {
  402. free((FREE_ARG) (v+nl-NR_END));
  403. }
  405. void free_dvector(double *v, long nl, long nh)
  406. /* free a double vector allocated with dvector() */
  407. {
  408. free((FREE_ARG) (v+nl-NR_END));
  409. }
  411. void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
  412. /* free a float matrix allocated by matrix() */
  413. {
  414. free((FREE_ARG) (m[nrl]+ncl-NR_END));
  415. free((FREE_ARG) (m+nrl-NR_END));
  416. }
  418. void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
  419. /* free a double matrix allocated by dmatrix() */
  420. {
  421. free((FREE_ARG) (m[nrl]+ncl-NR_END));
  422. free((FREE_ARG) (m+nrl-NR_END));
  423. }
  425. void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
  426. /* free an int matrix allocated by imatrix() */
  427. {
  428. free((FREE_ARG) (m[nrl]+ncl-NR_END));
  429. free((FREE_ARG) (m+nrl-NR_END));
  430. }
  432. void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
  433. /* free a submatrix allocated by submatrix() */
  434. {
  435. free((FREE_ARG) (b+nrl-NR_END));
  436. }
  438. void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
  439. /* free a matrix allocated by convert_matrix() */
  440. {
  441. free((FREE_ARG) (b+nrl-NR_END));
  442. }
  444. void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
  445. long ndl, long ndh)
  446. /* free a float f3tensor allocated by f3tensor() */
  447. {
  448. free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
  449. free((FREE_ARG) (t[nrl]+ncl-NR_END));
  450. free((FREE_ARG) (t+nrl-NR_END));
  451. }
  452. #define NRANSI
  453. #include "nrutil.h"
  454. #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
  455. #define M 7
  456. #define NSTACK 50
  458. void sort2(unsigned long n, float arr[], float brr[])
  459. {
  460. unsigned long i,ir=n,j,k,l=1;
  461. int *istack,jstack=0;
  462. float a,b,temp;
  464. istack=ivector(1,NSTACK);
  465. for (;;) {
  466. if (ir-l < M) {
  467. for (j=l+1;j<=ir;j++) {
  468. a=arr[j];
  469. b=brr[j];
  470. for (i=j-1;i>=l;i--) {
  471. if (arr[i] <= a) break;
  472. arr[i+1]=arr[i];
  473. brr[i+1]=brr[i];
  474. }
  475. arr[i+1]=a;
  476. brr[i+1]=b;
  477. }
  478. if (!jstack) {
  479. free_ivector(istack,1,NSTACK);
  480. return;
  481. }
  482. ir=istack[jstack];
  483. l=istack[jstack-1];
  484. jstack -= 2;
  485. } else {
  486. k=(l+ir) >> 1;
  487. SWAP(arr[k],arr[l+1])
  488. SWAP(brr[k],brr[l+1])
  489. if (arr[l] > arr[ir]) {
  490. SWAP(arr[l],arr[ir])
  491. SWAP(brr[l],brr[ir])
  492. }
  493. if (arr[l+1] > arr[ir]) {
  494. SWAP(arr[l+1],arr[ir])
  495. SWAP(brr[l+1],brr[ir])
  496. }
  497. if (arr[l] > arr[l+1]) {
  498. SWAP(arr[l],arr[l+1])
  499. SWAP(brr[l],brr[l+1])
  500. }
  501. i=l+1;
  502. j=ir;
  503. a=arr[l+1];
  504. b=brr[l+1];
  505. for (;;) {
  506. do i++; while (arr[i] < a);
  507. do j--; while (arr[j] > a);
  508. if (j < i) break;
  509. SWAP(arr[i],arr[j])
  510. SWAP(brr[i],brr[j])
  511. }
  512. arr[l+1]=arr[j];
  513. arr[j]=a;
  514. brr[l+1]=brr[j];
  515. brr[j]=b;
  516. jstack += 2;
  517. if (jstack > NSTACK) nrerror("NSTACK too small in sort2.");
  518. if (ir-i+1 >= j-l) {
  519. istack[jstack]=ir;
  520. istack[jstack-1]=i;
  521. ir=j-1;
  522. } else {
  523. istack[jstack]=j-1;
  524. istack[jstack-1]=l;
  525. l=i;
  526. }
  527. }
  528. }
  529. }
  530. #undef M
  531. #undef NSTACK
  532. #undef SWAP
  533. #undef NRANSI
  535. void crank(unsigned long n, float w[], float *s)
  536. {
  537. unsigned long j=1,ji,jt;
  538. float t,rank;
  540. *s=0.0;
  541. while (j < n) {
  542. if (w[j+1] != w[j]) {
  543. w[j]=j;
  544. ++j;
  545. } else {
  546. for (jt=j+1;jt<=n && w[jt]==w[j];jt++);
  547. rank=0.5*(j+jt-1);
  548. for (ji=j;ji<=(jt-1);ji++) w[ji]=rank;
  549. t=jt-j;
  550. *s += t*t*t-t;
  551. j=jt;
  552. }
  553. }
  554. if (j == n) w[n]=n;
  555. }
  557. float erfcc(float x)
  558. {
  559. float t,z,ans;
  561. z=fabs(x);
  562. t=1.0/(1.0+0.5*z);
  563. ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
  564. t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
  565. t*(-0.82215223+t*0.17087277)))))))));
  566. return x >= 0.0 ? ans : 2.0-ans;
  567. }
  569. float betai(float a, float b, float x)
  570. {
  571. float betacf(float a, float b, float x);
  572. float gammln(float xx);
  573. void nrerror(char error_text[]);
  574. float bt;
  576. if (x < 0.0 || x > 1.0) nrerror("Bad x in routine betai");
  577. if (x == 0.0 || x == 1.0) bt=0.0;
  578. else
  579. bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
  580. if (x < (a+1.0)/(a+b+2.0))
  581. return bt*betacf(a,b,x)/a;
  582. else
  583. return 1.0-bt*betacf(b,a,1.0-x)/b;
  584. }
  587. #define MAXIT 100
  588. #define EPS 3.0e-7
  589. #define FPMIN 1.0e-30
  591. float betacf(float a, float b, float x)
  592. {
  593. void nrerror(char error_text[]);
  594. int m,m2;
  595. float aa,c,d,del,h,qab,qam,qap;
  597. qab=a+b;
  598. qap=a+1.0;
  599. qam=a-1.0;
  600. c=1.0;
  601. d=1.0-qab*x/qap;
  602. if (fabs(d) < FPMIN) d=FPMIN;
  603. d=1.0/d;
  604. h=d;
  605. for (m=1;m<=MAXIT;m++) {
  606. m2=2*m;
  607. aa=m*(b-m)*x/((qam+m2)*(a+m2));
  608. d=1.0+aa*d;
  609. if (fabs(d) < FPMIN) d=FPMIN;
  610. c=1.0+aa/c;
  611. if (fabs(c) < FPMIN) c=FPMIN;
  612. d=1.0/d;
  613. h *= d*c;
  614. aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
  615. d=1.0+aa*d;
  616. if (fabs(d) < FPMIN) d=FPMIN;
  617. c=1.0+aa/c;
  618. if (fabs(c) < FPMIN) c=FPMIN;
  619. d=1.0/d;
  620. del=d*c;
  621. h *= del;
  622. if (fabs(del-1.0) < EPS) break;
  623. }
  624. if (m > MAXIT) nrerror("a or b too big, or MAXIT too small in betacf");
  625. return h;
  626. }
  627. #undef MAXIT
  628. #undef EPS
  629. #undef FPMIN
  631. //---------------------------------------------------------------------------
  633. #define NRANSI
  634. #include "nrutil.h"
  635. #define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
  636. #define M 7
  637. #define NSTACK 50
  639. void sort2i(unsigned long n, int arr[], int brr[])
  640. {
  641. unsigned long i,ir=n,j,k,l=1;
  642. int *istack,jstack=0;
  643. int a,b,temp;
  645. istack=ivector(1,NSTACK);
  646. for (;;) {
  647. if (ir-l < M) {
  648. for (j=l+1;j<=ir;j++) {
  649. a=arr[j];
  650. b=brr[j];
  651. for (i=j-1;i>=l;i--) {
  652. if (arr[i] <= a) break;
  653. arr[i+1]=arr[i];
  654. brr[i+1]=brr[i];
  655. }
  656. arr[i+1]=a;
  657. brr[i+1]=b;
  658. }
  659. if (!jstack) {
  660. free_ivector(istack,1,NSTACK);
  661. return;
  662. }
  663. ir=istack[jstack];
  664. l=istack[jstack-1];
  665. jstack -= 2;
  666. } else {
  667. k=(l+ir) >> 1;
  668. SWAP(arr[k],arr[l+1])
  669. SWAP(brr[k],brr[l+1])
  670. if (arr[l] > arr[ir]) {
  671. SWAP(arr[l],arr[ir])
  672. SWAP(brr[l],brr[ir])
  673. }
  674. if (arr[l+1] > arr[ir]) {
  675. SWAP(arr[l+1],arr[ir])
  676. SWAP(brr[l+1],brr[ir])
  677. }
  678. if (arr[l] > arr[l+1]) {
  679. SWAP(arr[l],arr[l+1])
  680. SWAP(brr[l],brr[l+1])
  681. }
  682. i=l+1;
  683. j=ir;
  684. a=arr[l+1];
  685. b=brr[l+1];
  686. for (;;) {
  687. do i++; while (arr[i] < a);
  688. do j--; while (arr[j] > a);
  689. if (j < i) break;
  690. SWAP(arr[i],arr[j])
  691. SWAP(brr[i],brr[j])
  692. }
  693. arr[l+1]=arr[j];
  694. arr[j]=a;
  695. brr[l+1]=brr[j];
  696. brr[j]=b;
  697. jstack += 2;
  698. if (jstack > NSTACK) nrerror("NSTACK too small in sort2.");
  699. if (ir-i+1 >= j-l) {
  700. istack[jstack]=ir;
  701. istack[jstack-1]=i;
  702. ir=j-1;
  703. } else {
  704. istack[jstack]=j-1;
  705. istack[jstack-1]=l;
  706. l=i;
  707. }
  708. }
  709. }
  710. }
  711. #undef M
  712. #undef NSTACK
  713. #undef SWAP
  714. #undef NRANSI
  716. //---------------------------------------------------------------------------
  717. void svdcmp(float **a, int m, int n, float w[], float **v)
  718. {
  719. float pythag(float a, float b);
  720. int flag,i,its,j,jj,k,l,nm;
  721. float anorm,c,f,g,h,s,scale,x,y,z,*rv1;
  723. rv1=vector(1,n);
  724. g=scale=anorm=0.0;
  725. for (i=1;i<=n;i++) {
  726. l=i+1;
  727. rv1[i]=scale*g;
  728. g=s=scale=0.0;
  729. if (i <= m) {
  730. for (k=i;k<=m;k++) scale += fabs(a[k][i]);
  731. if (scale) {
  732. for (k=i;k<=m;k++) {
  733. a[k][i] /= scale;
  734. s += a[k][i]*a[k][i];
  735. }
  736. f=a[i][i];
  737. g = -SIGN(sqrt(s),f);
  738. h=f*g-s;
  739. a[i][i]=f-g;
  740. for (j=l;j<=n;j++) {
  741. for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
  742. f=s/h;
  743. for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
  744. }
  745. for (k=i;k<=m;k++) a[k][i] *= scale;
  746. }
  747. }
  748. w[i]=scale *g;
  749. g=s=scale=0.0;
  750. if (i <= m && i != n) {
  751. for (k=l;k<=n;k++) scale += fabs(a[i][k]);
  752. if (scale) {
  753. for (k=l;k<=n;k++) {
  754. a[i][k] /= scale;
  755. s += a[i][k]*a[i][k];
  756. }
  757. f=a[i][l];
  758. g = -SIGN(sqrt(s),f);
  759. h=f*g-s;
  760. a[i][l]=f-g;
  761. for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
  762. for (j=l;j<=m;j++) {
  763. for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
  764. for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
  765. }
  766. for (k=l;k<=n;k++) a[i][k] *= scale;
  767. }
  768. }
  769. anorm=FMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
  770. }
  771. for (i=n;i>=1;i--) {
  772. if (i < n) {
  773. if (g) {
  774. for (j=l;j<=n;j++)
  775. v[j][i]=(a[i][j]/a[i][l])/g;
  776. for (j=l;j<=n;j++) {
  777. for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
  778. for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
  779. }
  780. }
  781. for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
  782. }
  783. v[i][i]=1.0;
  784. g=rv1[i];
  785. l=i;
  786. }
  787. for (i=IMIN(m,n);i>=1;i--) {
  788. l=i+1;
  789. g=w[i];
  790. for (j=l;j<=n;j++) a[i][j]=0.0;
  791. if (g) {
  792. g=1.0/g;
  793. for (j=l;j<=n;j++) {
  794. for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
  795. f=(s/a[i][i])*g;
  796. for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
  797. }
  798. for (j=i;j<=m;j++) a[j][i] *= g;
  799. } else for (j=i;j<=m;j++) a[j][i]=0.0;
  800. ++a[i][i];
  801. }
  802. for (k=n;k>=1;k--) {
  803. for (its=1;its<=30;its++) {
  804. flag=1;
  805. for (l=k;l>=1;l--) {
  806. nm=l-1;
  807. if ((float)(fabs(rv1[l])+anorm) == anorm) {
  808. flag=0;
  809. break;
  810. }
  811. if ((float)(fabs(w[nm])+anorm) == anorm) break;
  812. }
  813. if (flag) {
  814. c=0.0;
  815. s=1.0;
  816. for (i=l;i<=k;i++) {
  817. f=s*rv1[i];
  818. rv1[i]=c*rv1[i];
  819. if ((float)(fabs(f)+anorm) == anorm) break;
  820. g=w[i];
  821. h=pythag(f,g);
  822. w[i]=h;
  823. h=1.0/h;
  824. c=g*h;
  825. s = -f*h;
  826. for (j=1;j<=m;j++) {
  827. y=a[j][nm];
  828. z=a[j][i];
  829. a[j][nm]=y*c+z*s;
  830. a[j][i]=z*c-y*s;
  831. }
  832. }
  833. }
  834. z=w[k];
  835. if (l == k) {
  836. if (z < 0.0) {
  837. w[k] = -z;
  838. for (j=1;j<=n;j++) v[j][k] = -v[j][k];
  839. }
  840. break;
  841. }
  842. if (its == 30) nrerror("no convergence in 30 svdcmp iterations");
  843. x=w[l];
  844. nm=k-1;
  845. y=w[nm];
  846. g=rv1[nm];
  847. h=rv1[k];
  848. f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
  849. g=pythag(f,1.0);
  850. f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
  851. c=s=1.0;
  852. for (j=l;j<=nm;j++) {
  853. i=j+1;
  854. g=rv1[i];
  855. y=w[i];
  856. h=s*g;
  857. g=c*g;
  858. z=pythag(f,h);
  859. rv1[j]=z;
  860. c=f/z;
  861. s=h/z;
  862. f=x*c+g*s;
  863. g = g*c-x*s;
  864. h=y*s;
  865. y *= c;
  866. for (jj=1;jj<=n;jj++) {
  867. x=v[jj][j];
  868. z=v[jj][i];
  869. v[jj][j]=x*c+z*s;
  870. v[jj][i]=z*c-x*s;
  871. }
  872. z=pythag(f,h);
  873. w[j]=z;
  874. if (z) {
  875. z=1.0/z;
  876. c=f*z;
  877. s=h*z;
  878. }
  879. f=c*g+s*y;
  880. x=c*y-s*g;
  881. for (jj=1;jj<=m;jj++) {
  882. y=a[jj][j];
  883. z=a[jj][i];
  884. a[jj][j]=y*c+z*s;
  885. a[jj][i]=z*c-y*s;
  886. }
  887. }
  888. rv1[l]=0.0;
  889. rv1[k]=f;
  890. w[k]=x;
  891. }
  892. }
  893. free_vector(rv1,1,n);
  894. }
  896. float pythag(float a, float b)
  897. {
  898. float absa,absb;
  899. absa=fabs(a);
  900. absb=fabs(b);
  901. if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa));
  902. else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
  903. }
  906. void kendl1(float data1[], float data2[], unsigned long n, float *tau,
  907. float *z, float *prob)
  908. {
  909. float erfcc(float x);
  910. unsigned long n2=0,n1=0,k,j;
  911. long is=0;
  912. float svar,aa,a2,a1;
  914. for (j=1;j<n;j++) {
  915. for (k=(j+1);k<=n;k++) {
  916. a1=data1[j]-data1[k];
  917. a2=data2[j]-data2[k];
  918. aa=a1*a2;
  919. if (aa) {
  920. ++n1;
  921. ++n2;
  922. aa > 0.0 ? ++is : --is;
  923. } else {
  924. if (a1) ++n1;
  925. if (a2) ++n2;
  926. }
  927. }
  928. }
  929. *tau=is/(sqrt((double) n1)*sqrt((double) n2));
  930. svar=(4.0*n+10.0)/(9.0*n*(n-1.0));
  931. *z=(*tau)/sqrt(svar);
  932. *prob=erfcc(fabs(*z)/1.4142136);
  933. }