double choose(int n,int k);
double loggam(double n);
float factrl(int n);
double betaFcn(double alpha,double beta);
double gammp(double a, double x);
void gser(double *gamser,double a,double x,double *gln);
void gcf(double *gammcf,double a,double x,double *gln);
double qgaus(double (*func)(double),double a,double b);
double trapzd(double (*func)(double),double a,double b,int n);
double qtrap(double (*func)(double),double a,double b);
double qromb(double (*func)(double),double a,double b);
double quad2d(double (*func)(double,double,double,double,double,
   double,double,double,double,double), double x1,
   double x2,double W,double y2,double n2,double y1,
   double n1,double m0,double C0,double tt);
void polint(double *xa,double *ya,int n,double x,double *y,double *dy);
void spline(double *x, double *y, int n,double yp1, double ypn,
   double *y2);
void splint(double *xa, double *ya, double *y2a, int n, double x,
   double *y);
void splie2(double *x1a,double *x2a,double **ya,int m,int n,double **y2a);
void splin2(double *x1a,double *x2a,double **ya,double **y2a,int m,int n,
   double x1,double x2,double *y);
double zbrent(double (*func)(double,double **),double **rhs,double x1,
   double x2,double tol);
void zbrak(double (*fx)(double,double **),double **rhs,double x1,
   double x2,int n,double *xb1,double *xb2,int *nb);
void tridag(double *a,double *b,double *c,double *r,double *u,long n);
double gausCDF(double x);
double erfcc(double x);
double kIntegrate(int k,double (*func)(double *));
double func1(double *x);

double a=-3.0,b=3.0,kGlob=12;

double choose(int n,int k)
{
   if((k<0) | (n<0)){
      nrerror("Error in choose: Negative arg.s");
   }
   return factrl(n)/(factrl(k)*factrl(n-k));
}

double loggam(double xx)
{
   double x,y,tmp,ser;
   static double cof[6]={76.18009172947146, -86.50532032941677,
      24.01409824083091,-1.231739572450155, 0.001208650973866179,
      -5.395239384953e-006};
   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.506628274631*ser/x);
}

float factrl(int n)
{
   double loggam(double xx);
   static int ntop=4;
   static double aaa[33]={1.0,1.0,2.0,6.0,24.0};
   int j;
   if(n < 0) nrerror("Negative arg to factrl");
   if(n > 32) return exp(loggam(n+1.0));
   while (ntop<n){
      j=ntop++;
      aaa[ntop]=aaa[j]*ntop;
   }
   return aaa[n];
}

double betaFcn(double alpha,double beta)
{
   if(alpha<=0)
      nrerror("Error in betaFcn: alpha<0");
   if(beta<=0)
      nrerror("Error in betaFcn: beta<0");
   return exp(loggam(alpha)+loggam(beta)-loggam(alpha+beta));
}

double gammp(double a, double x)
{
/* high level function for incomplete gamma function */
   void gcf(double *gammcf,double a,double x,double *gln);
   void gser(double *gamser,double a,double x,double *gln);
   double gamser,gammcf,gln;
   if(x < 0.0 || a <= 0.0) nrerror("Invalid arg in gammp");
   if(x < (a+1.0)){
/* here I change routine so that it returns \gamma(a,x)
   or P(a,x)-just take out comments to get P(a,x) vs \gamma(a,x)-
   to get latter use the exp(log(.)+gln) expression */
      gser(&gamser,a,x,&gln);
//      return exp(log(gamser)+gln);
      return gamser;
   }
   else{
      gcf(&gammcf,a,x,&gln);
//      return exp(log(1.0-gammcf)+gln);
      return 1.0-gammcf;
   }
}

#define ITMAX 100
#define EPSW 3.0e-7

void gser(double *gamser,double a,double x,double *gln)
{
   double loggam(double xx);
   int n;
   double sum,del,ap;
   *gln=loggam(a);
   if(x <= 0.0){
      if(x < 0.0) nrerror("x less than 0 in routine gser");
      *gamser=0.0;
      return;
   }
   else{
      ap=a;
      del=sum=1.0/a;
      for(n=1;n<=ITMAX;n++){
         ++ap;
         del *= x/ap;
         sum += del;
         if(fabs(del) < fabs(sum)*EPSW){
            *gamser=sum*exp(-x+a*log(x)-(*gln));
            return;
         }
      }
      nrerror("a too large, ITMAX too small in routine gser");
      return;
   }
}


#define FPMIN 1.0e-30

void gcf(double *gammcf,double a,double x,double *gln)
{
   int i;
   double an,b,c,d,del,h;
   *gln=loggam(a);
   b=x+1.0-a;
   c=1.0/FPMIN;
   d=1.0/b;
   h=d;
   for(i=1;i<=ITMAX;i++){
      an = -i*(i-a);
      b += 2.0;
      d=an*d+b;
      if(fabs(d) < FPMIN) d=FPMIN;
      c=b+an/c;
      if(fabs(c) < FPMIN) c=FPMIN;
      d=1.0/d;
      del=d*c;
      h *= del;
      if(fabs(del-1.0) < EPSW) break;
   }
   if(i > ITMAX) nrerror("a too large, ITMAX too small in gcf");
   *gammcf=exp(-x+a*log(x)-(*gln))*h;
}
#undef ITMAX

double qgaus(double (*func)(double),double a,double b)
{
   int j;
   double xr,xm,dx,s;
   static double x[]={0.0,0.1488743389,0.4333953941,
      0.6794095682,0.8650633666,0.9739065285};
   static double w[]={0.0,0.2955242247,0.2692667193,
      0.2190863625,0.1494513491,0.0666713443};
   xm=0.5*(b+a);
   xr=0.5*(b-a);
   s=0;
   for(j=1;j<=5;j++){
      dx=xr*x[j];
      s += w[j]*((*func)(xm+dx)+
      (*func)(xm-dx));
   }
   return s *= xr;
}

/* functions for integration using Romberg's method */
double qromb(double (*func)(double),double a,double b)
{
   double trapzd(double (*func)(double),double a,double b,int n);
   void polint(double *xa,double *ya,int n,double x,double *y,double *dy);
   int j,K=5;
   const int JMAX=20,JMAXP=JMAX+1;
   /* this EPS controls accuarcy of numerical integral */
   double ss,dss,EPS=1.0e-6;
   double s[JMAXP],h[JMAXP+1];
   h[1]=1.0;
   for(j=1;j<=JMAX;j++){
      s[j]=trapzd(func,a,b,j);
      if(j >= K){
         polint(&h[j-K],&s[j-K],K,0.0,&ss,&dss);
         if(fabs(dss) <= EPS*fabs(ss)) return ss;
      }
      h[j+1]=0.25*h[j];
   }
   nrerror("too many steps in qromb");
   return 0.0;
}

/* functions for integration using the trapezoidal rule */
double qtrap(double (*func)(double),double a,double b)
{
   double trapzd(double (*func)(double),double a,double b,int n);
   int j,JMAX=20;
   /* this EPS controls accuarcy of numerical integral */
   double s,olds,EPS=1.0e-3;
   olds= -1.0e30;
   for(j=1;j<=JMAX;j++){
      s=trapzd(func,a,b,j);
      if(fabs(s-olds) < EPS*fabs(olds)) return s;
      if(s==0.0 && olds == 0.0 && j>6) return s;
      olds=s;
   }
   nrerror("too many steps in routine qtrap");
   return 0.0;
}

#define FUNC(x) ((*func)(x))

double trapzd(double (*func)(double),double a,double b,int n)
{
   double x,tnm,sum,del;
   static double s;
   int it,j;
   if(n==1){
      return (s=0.5*(b-a)*(FUNC(a)+FUNC(b)));
   } else{
      for(it=1,j=1;j<n-1;j++) it <<= 1;
      tnm=it;
      del=(b-a)/tnm;
      x=a+0.5*del;
      for(sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x);
      s=0.5*(s+(b-a)*sum/tnm);
      return s;
   }
}

/* routines for integrating over the square [0,1]^2 in 2 dimensions */
static double xsav;
static double (*nrfuncq)(double,double);

double quad2d(double (*func)(double,double), double x1,double x2)
{
   double qgaus(double (*func)(double),double a,double b);
   double f1(double x);
   nrfuncq=func;
   return qgaus(f1,x1,x2);
}

double f1(double x)
{
   double qgaus(double (*func)(double),double a,double b);
   double f2(double y);
   xsav=x;
   return qgaus(f2,0.0,1.0);
}

double f2(double y)
{
   return (*nrfuncq)(xsav,y);
}

void polint(double *xa,double *ya,int n,double x,double *y,double *dy)
{
   int i,m,ns=1;
   double den,dif,dift,ho,hp,w;
   double *c,*d;
   dif=fabs(x-xa[1]);
   c=dvector(1,n);
   d=dvector(1,n);
   for(i=1;i<=n;i++){
      if((dift=fabs(x-xa[i]))<dif){
         ns=i;
         dif=dift;
      }
      c[i]=ya[i];
      d[i]=ya[i];
   }
   *y=ya[ns--];
   for(m=1;m<n;m++){
      for(i=1;i<=n-m;i++){
         ho=xa[i]-x;
         hp=xa[i+m]-x;
         w=c[i+1]-d[i];
         if((den=ho-hp)==0.0) nrerror("error in polint");
         den=w/den;
         d[i]=hp*den;
         c[i]=ho*den;
      }
      *y += (*dy=(2*ns < (n-m) ? c[ns+1] : d[ns--]));
   }
   free_dvector(d,1,n);
   free_dvector(c,1,n);
}

void spline(double *x, double *y, int n, double yp1, double ypn,
   double *y2)
{
   int i,k;
   double p,qn,sig,un,*u;
   u=dvector(1,n-1);
   if(yp1 > 0.99e30)
      y2[1]=u[1]=0.0;
   else{
      y2[1]=-0.5;
      u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
   }
   for(i=2;i<=n-1;i++){
      sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
      p=sig*y2[i-1]+2.0;
      y2[i]=(sig-1.0)/p;
      u[i]=(y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
      u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
   }
   if(ypn > 0.99e30)
      qn=un=0.0;
   else{
      qn=0.5;
      un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
   }
   y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
   for(k=n-1;k>=1;k--)
      y2[k]=y2[k]*y2[k+1]+u[k];
   free_dvector(u,1,n-1);
}

void splint(double *xa, double *ya, double *y2a, int n, double x,
   double *y)
{
   int klo,khi,k;
   double h,b,a;
   klo=1;
   khi=n;
   while(khi-klo > 1){
      k=(khi+klo) >> 1;
      if(xa[k] > x) khi=k;
      else klo=k;
   }
   h=xa[khi]-xa[klo];
   if(h == 0.0) nrerror("Bad xa input to routine splint:ties!");
   a=(xa[khi]-x)/h;
   b=(x-xa[klo])/h;
   *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+
      (b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

void splie2(double *x1a,double *x2a,double **ya,int m,int n,double **y2a)
{
   int j;
   for(j=1;j<=m;j++)
      spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]);
}

void splin2(double *x1a,double *x2a,double **ya,double **y2a,int m,int n,
   double x1,double x2,double *y)
{
   int j;
   double *ytmp,*yytmp;
   ytmp=dvector(1,m);
   yytmp=dvector(1,m);
   for(j=1;j<=m;j++)
      splint(x2a,ya[j],y2a[j],n,x2,&yytmp[j]);
   spline(x1a,yytmp,m,1.0e30,1.0e30,ytmp);
   splint(x1a,yytmp,ytmp,m,x1,y);
   free_dvector(yytmp,1,m);
   free_dvector(ytmp,1,m);
}

/* a 1-D root bracketer */
void zbrak(double (*fx)(double,double **),double **rhs,double x1,
   double x2,int n,double *xb1,double *xb2,int *nb)
{
   int nbb,i;
   double x,fp,fc,dx;
   nbb=0;
   dx=(x2-x1)/n;
   fp=(*fx)(x=x1,rhs);
   for(i=1;i<=n;i++){
      fc=(*fx)(x += dx,rhs);
      if(fc*fp <= 0.0){
         xb1[++nbb]=x-dx;
         xb2[nbb]=x;
         if(*nb == nbb) return;
      }
      fp=fc;
   }
   *nb=nbb;
}
#include "/misc/radon2/xiaoxiak/thesis/header/nrutil.h"
/* a 1-D nonlinear equation root finder */
double zbrent(double (*func)(double,double **),double **rhs,double x1,
   double x2,double tol)
{
   int itmax=100,iter;
   double EPS=3.0e-8,a=x1,b=x2,c=x2,d,e,min1,min2;
   double fa=(*func)(a,rhs),fb=(*func)(b,rhs),fc,p,q,r,s,tol1,xm;
   if((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
      nrerror("Root must be bracketed in zbrent");
   fc=fb;
   for(iter=1;iter<=itmax;iter++){
      if((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)){
         c=a;
         fc=fa;
         e=d=b-a;
      }
      if(fabs(fc) < fabs(fb)){
         a=b;
         b=c;
         c=a;
         fa=fb;
         fb=fc;
         fc=fa;
      }
      tol1=2.0*EPS*fabs(b)+0.5*tol;
      xm=0.5*(c-b);
      if(fabs(xm) <= tol1 || fb == 0.0) return b;
      if(fabs(e) >= tol1 && fabs(fa) > fabs(fb)){
         s=fb/fa;
         if(a == c){
            p=2.0*xm*s;
            q=1.0-s;
         } else{
            q=fa/fc;
            r=fb/fc;
            p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
            q=(q-1.0)*(r-1.0)*(s-1.0);
         }
         if(p > 0.0) q = -q;
         fabs(p);
         min1=3.0*xm*q-fabs(tol1*q);
         min2=fabs(e*q);
         if(2.0*p < (min1 < min2 ? min1 : min2)){
            e=d;
            d=p/q;
         } else{
            d=xm;
            e=d;
         }
      } else{
         d=xm;
         e=d;
      }
      a=b;
      fa=fb;
      if(fabs(d) > tol1)
         b += d;
      else
         b += SIGN(tol1,xm);
      fb=(*func)(b,rhs);
   }
   nrerror("Maximum number of iterations exceeded in zbrent");
   return 0.0;
}

void tridag(double *a,double *b,double *c,double *r,double *u,long n)
{
   long j;
   double bet,*gam;
   gam=dvector(1,n);
   if(b[1]==0.0) nrerror("error 1 in tridag");
   u[1]=r[1]/(bet=b[1]);
   for(j=2;j<=n;j++){
      gam[j]=c[j-1]/bet;
      bet=b[j]-a[j]*gam[j];
      if(bet==0.0) nrerror("error 2 in tridag");
      u[j]=(r[j]-a[j]*u[j-1])/bet;
   }
   for(j=(n-1);j>=1;j--)
      u[j] -= gam[j+1]*u[j+1];
   free_dvector(gam,1,n);
}

double gausCDF(double x)
{
   double res;
//   if(x>=0.0) res=0.5*(1.0+gammp(0.5,0.5*x*x));
//   else res=0.5*(1.0-gammp(0.5,0.5*x*x));
   if(x>=0.0) res=1.0-0.5*erfcc(x/sqrt(2.0));
   else res=0.5*erfcc(-x/sqrt(2.0));
   return res;
} 

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;
}

double func1(double *x)
{
   int i;
   double s=0.0;
   for(i=1;i<=kGlob;i++) s += x[i]*x[i];
   return exp(-0.5*s);
//   return 1.0;
}


double kIntegrate(int k,double (*func)(double *))
{
   int i,j,*indx;
   double temp,xm,xr,dx,res=0.0,*xVec,*wVec,*arg;
   xVec=dvector(1,10);
   wVec=dvector(1,10);
   arg=dvector(1,k);
   indx=ivector(1,k);
   double x[]={0.0,0.1488743389,0.4333953941,
      0.6794095682,0.8650633666,0.9739065285};
   double w[]={0.0,0.2955242247,0.2692667193,
      0.2190863625,0.1494513491,0.0666713443};
   xm=0.5*(b+a);
   xr=0.5*(b-a);
   for(i=1;i<=5;i++) xVec[i]=xm+xr*x[i];
   for(i=6;i<=10;i++) xVec[i]=xm-xr*x[i-5];
   for(i=1;i<=5;i++) wVec[i]=wVec[i+5]=w[i];
   for(i=1;i<=pow(10.0,k);i++){
      indx[1]=(i-1)%10+1;
      arg[1]=xVec[indx[1]];
      for(j=2;j<=k;j++){
         indx[j]=int((i-1)/pow(10,j-1))%10+1;
         arg[j]=xVec[indx[j]];
      }
      temp=func(arg);
      for(j=1;j<=k;j++) temp *= wVec[indx[j]];
      res+=temp;
   }
   free_dvector(arg,1,k);
   free_dvector(xVec,1,10);
   free_dvector(wVec,1,10);
   free_ivector(indx,1,k);
   return pow(xr,k)*res;
}