/* Statistical methods library */
/* linear fits, qudratic fits, basic statistics */
/* in a program with x[N], y[N], a,b,sa,sb call with  */
/* linear_fit(&a, &b, &sa, &sb, x,y,N); */
/* chebyshev fit and definite integration */
/* solution to linear equations, determinants */
/* chebyshev derivatives */
/* J. R. Schmidt, 2001/2002/2003  */


/*declarations */

#define MAX 100
#include<gmp.h>

/* prototypes */

void standard (double *mean, double *sd, double *x, int N);
void linear_fit (double *a, double *b, double *sa, double *sb, double *x,
		 double *y, int N);
double quad_fit (double *x0, double *v0, double *a, double *x, double *y,
		 int N);
double sigma (double x0, double v0, double a, double *x, double *y, int N);
double determinant (double (*A)[MAX], int N);
void sums (double *sx, double *sy, double *sx2, double *sxy, double *x, double *y, int N);
/* compute derivative from chbyshev fit */
void derivative (double *xp, double *yp, double *derptr, int M, int NN);
/* chebyshev polynomials */
double T (double x, int n);
/* fit to chebyshev function expansion */
/* NN is order of fit, M is length of data arrays */
/* \int_{xmin}^{xmax} f(x) dx from \{(x_0, y_0),...,(x_{M-1},y_{M-1})\} */
void chebyshev_int (double *integ, double *x_min, double *x_max, double *xp,
		    double *yp, int M, int NN);
/* solve a\cdot x=b */
void solve (double (*A)[MAX], double *B, double *X, int N);
/* chi-squared statistic and probability */
double chi_2_prob(double val, int N);
double chisquared(double *exper, double *theor, int N);
/* error function for normal distributions */
double erf(double z);
double poisson(double mean, int M);
double normal_dist(double xmin, double xmax, double mean, double sig);



/* The routines */

void
linear_fit (double *a, double *b, double *sa, double *sb, double *x,
	    double *y, int N)
{
  double Sx, Sy, Sx2, Sxy, sig, A, B, SA, SB, Den, term;
  int i, j;
  Sx = 0.0;
  Sy = 0.0;
  Sx2 = 0.0;
  Sxy = 0.0;
  sig = 0.0;
  for (i = 0; i < N; i++)
    {
      Sx = Sx + *(x + i);
      Sy = Sy + *(y + i);
      Sx2 = Sx2 + (*(x + i)) * (*(x + i));
      Sxy = Sxy + (*(x + i)) * (*(y + i));
    }
  Den = (double) N *Sx2 - Sx * Sx;
  if (Den != 0.0)
    {
      A = ((double) N * Sxy - Sx * Sy) / Den;
      B = -(Sx * Sxy - Sx2 * Sy) / Den;
    }
  else
    {
      A = 1.0E10;
      B = 1.0E10;
    }

  for (i = 0; i < N; i++)
    {
      term = *(y + i) - (A * (*(x + i)) + B);
      term = term * term;
      sig = sig + term;
    }
  sig = sig / (double) (N - 1);
  

  if (Den != 0.0)
    {
      SB = sig * Sx2 / Den;
      SA = sig * (double) N / Den;
    }
  else
    {
      SA = 1.0E10;
      SB = 1.0E10;
    }
  *a = A;
  *b = B;
  *sa = sqrt (SA);
  *sb = sqrt (SB);


}

double chisquared(double *exper, double *theor, int N){
    int i;
    double sum, term;
    sum=0.0;
    for(i=0;i<N;i++){
	term= *(theor+i) - *(exper+i);
	term=term*term;
	if( *(theor+i) != 0.0){
	    term=term/( *(theor+i) );
	    sum=sum+term;
	}
    }
    return(sum);
}  
	


void
standard (double *mean, double *sd, double *x, int N)
{
  double sum, sum2, ave, standev, term;
  int i;
  sum = 0.0;
  sum2 = 0.0;
  standev = 0.0;
  for (i = 0; i < N; i++)
    {
      sum = sum + *(x + i);
    }
  ave = sum / (double) N;
  for (i = 0; i < N; i++)
    {
      term = *(x + i) - ave;
      term = term * term;
      sum2 = sum2 + term;
    }
  standev = sum2 / (double) (N - 1);
  *mean = ave;
  *sd = sqrt (standev);
}

double
quad_fit (double *x0, double *v0, double *a, double *x, double *y, int N)
{
  /* fits to a quadratic x_0+v_0 t +1/2 a t^2 */
  double Sx, Sy, Sx2, Sx3, Sx4, Sxy, Sx2y;
  double numerator, denominator;
  int i, j;
  Sx = 0.0;
  Sy = 0.0;
  Sx2 = 0.0;
  Sxy = 0.0;
  Sx3 = 0.0;
  Sx4 = 0.0;
  Sx2y = 0.0;
  for (i = 0; i < N; i++)
    {
      Sx = Sx + *(x + i);
      Sy = Sy + *(y + i);
      Sx2 = Sx2 + (*(x + i)) * (*(x + i));
      Sx3 = Sx3 + (*(x + i)) * (*(x + i)) * (*(x + i));
      Sx4 = Sx4 + (*(x + i)) * (*(x + i)) * (*(x + i)) * (*(x + i));
      Sxy = Sxy + (*(x + i)) * (*(y + i));
      Sx2y = Sx2y + (*(x + i)) * (*(y + i)) * (*(x + i));
    }
  denominator =
    (double) N *Sx2 * Sx4 + 2.0 * Sx * Sx3 * Sx2 - Sx2 * Sx2 * Sx2 -
    (double) N *Sx3 * Sx3 - Sx4 * Sx * Sx;
  numerator =
    Sy * Sx2 * Sx4 + Sx * Sx3 * Sx2y + Sxy * Sx3 * Sx2 - Sx2y * Sx2 * Sx2 -
    Sy * Sx3 * Sx3 - Sxy * Sx * Sx4;
  if (denominator != 0.0)
    *x0 = numerator / denominator;
  else
    *x0 = 1.0E10;
  numerator =
    (double) N *Sxy * Sx4 + Sy * Sx3 * Sx2 + Sx2 * Sx * Sx2y -
    Sx2 * Sxy * Sx2 - Sx2y * Sx3 * (double) N - Sx4 * Sx * Sy;
  if (denominator != 0.0)
    *v0 = numerator / denominator;
  else
    *v0 = 1.0E10;
  numerator =
    (double) N *Sx2 * Sx2y + Sx * Sxy * Sx2 + Sy * Sx * Sx3 - Sx2 * Sx2 * Sy -
    Sx3 * Sxy * (double) N - Sx2y * Sx * Sx;
  if (denominator != 0.0)
    *a = 2.0 * numerator / denominator;
  else
    *a = 1.0E10;

}




double
T (double x, int n)
{
  double T0, T1, T2;
  int i;
  T0 = 1.0;
  T1 = x;
  if (n == 0)
    return (T0);
  else
    {
      if (n == 1)
	return (T1);
      else
	{
	  for (i = 2; i <= n; i++)
	    {
	      T2 = 2.0 * x * T1 - T0;
	      T0 = T1;
	      T1 = T2;
	    }
	  return (T2);
	}
    }
}

void
chebyshev_int (double *integ, double *x_min, double *x_max, double *xp,
	       double *yp, int M, int NN)
{

  int k, i, j;
  double X[MAX], Y[MAX], a[MAX][MAX], b[MAX], c[MAX], xmin, xmax, integral,
    scale, offset;
  double (*a_ptr)[MAX], *b_ptr, *c_ptr;
  /* NN is order of polynomial fit */
  /* M data points  */
  k = 0;

  for (k = 0; k < M; k++)
    {
      X[k] = *(xp + k);
      Y[k] = *(yp + k);
    }

  xmin = X[0];
  xmax = X[0];
  for (k = 1; k < M; k++)
    {
      if (X[k] < xmin)
	xmin = X[k];
      if (X[k] > xmax)
	xmax = X[k];
    }
  scale = 2.0 / (xmax - xmin);
  offset = -(xmax + xmin) / (xmax - xmin);
  for (k = 0; k < M; k++)
    {
      X[k] = scale * X[k] + offset;
      /* now -1 <= X[k] <= 1 */
    }

  for (i = 0; i <= NN; i++)
    {
      for (j = 0; j <= NN; j++)
	{
	  a[i][j] = 0.0;
	}
    }
  /* set up a and b for polynomial fit of order N */
  for (j = 0; j <= NN; j++)
    {
      b[j] = 0.0;
      /*   x[j] = 0.0; */
      for (k = 0; k < M; k++)
	{

	  b[j] = b[j] + Y[k] * T (X[k], j);
	}
    }

  for (i = 0; i <= NN; i++)
    {
      for (j = 0; j <= NN; j++)
	{
	  for (k = 0; k < M; k++)
	    a[i][j] = a[i][j] + T (X[k], j) * T (X[k], i);

	}
    }



  b_ptr = &b[0];		/* initialize the first addresses */
  c_ptr = &c[0];
  a_ptr = &a[0];
  solve (a_ptr, b_ptr, c_ptr, NN);	/* solve the equations */

  integral = 2.0 * (*(c_ptr + 0));
  for (k = 1; k <= NN; k++)
    {
      if (2 * (k / 2) == k)
	{
	  integral =
	    integral - 2.0 * (*(c_ptr + k)) / ((double) k * (double) k - 1.0);
	}
    }
  integral = integral / scale;


  /* fill pointers for output */
  *x_min = xmin;
  *x_max = xmax;
  *integ = integral;
}


double
determinant (double (*A)[MAX], int N)
{
  /* the routine for solving the simultaneous equations */
  double temp, factor, product;
  int m, n, u, p;
  int i, j, k, l;
  for (k = 0; k <= N; k++)
    {
      m = k;
      u = 0;
      do
	{
	  u + 1;
	}
      while (*(*(A + k) + u + k) == 0.0);
      m = k + u;
      for (l = m + 1; l <= N; l++)
	{
	  if (*(*(A + k) + m) != 0.0)
	    {
	      factor = *(*(A + k) + l) / *(*(A + k) + m);
	      /* perform row operation on [a] */
	      for (j = 0; j <= N; j++)
		{
		  *(*(A + j) + l) =
		    *(*(A + j) + l) - factor * (*(*(A + j) + m));
		}
	    }
	}
      for (j = 0; j <= N; j++)
	{
	  temp = *(*(A + j) + k);
	  *(*(A + j) + k) = *(*(A + j) + m);
	  *(*(A + j) + m) = temp;
	}
    }

  /* compute the determinant */
  product = 1.0;
  for (i = 0; i <= N; i++)
    product = product * (*(*(A + i) + i));
  return (product);
}


void
solve (double (*A)[MAX], double *B, double *X, int N)
{
  /* the routine for solving the simultaneous equations */
  double temp, factor, sum;
  int m, n, p, i, j, k, l, u;

  for (k = 0; k <= N; k++)
    {
      m = k;
      u = 0;
      do
	{
	  u + 1;
	}
      while (*(*(A + k) + u + k) == 0.0);
      m = k + u;
      for (l = m + 1; l <= N; l++)
	{
	  if (*(*(A + k) + m) != 0.0)
	    {
	      factor = *(*(A + k) + l) / *(*(A + k) + m);
	      /* perform row operation on [a] */
	      for (j = 0; j <= N; j++)
		{
		  *(*(A + j) + l) =
		    *(*(A + j) + l) - factor * (*(*(A + j) + m));
		}

/* perform row operation on [b] */
	      *(B + l) = *(B + l) - factor * (*(B + m));
	    }
	}
      for (j = 0; j <= N; j++)
	{
	  temp = *(*(A + j) + k);
	  *(*(A + j) + k) = *(*(A + j) + m);
	  *(*(A + j) + m) = temp;
	}
      /* do the switch for [b] */
      temp = *(B + k);
      *(B + k) = *(B + m);
      *(B + m) = temp;
    }

  /* compute the solutions x[i] */
  *(X + N) = *(B + N) / (*(*(A + N) + N));
  for (p = 0; p <= N; p++)
    {
      sum = 0.0;
      for (u = N - p + 1; u <= N; u++)
	sum = sum + (*(*(A + u) + N - p)) * (*(X + u));
      *(X + N - p) = (*(B + N - p) - sum) / (*(*(A + N - p) + N - p));
    }

}


void
derivative (double *xp, double *yp, double *derptr, int M, int NN)
{
    double left,right,change,xmin,xmax,X[MAX],Y[MAX],offset,scale;
    double a[MAX][MAX], b[MAX], c[MAX], (*a_ptr)[MAX], *b_ptr,*c_ptr;
    int i,j,k,l;
    
    change=0.001;

  /* NN is order of polynomial fit */
  /* M data points  */
    /* we will fit a chebyshev sum, and differentiate it */
  k = 0;

  for (k = 0; k < M; k++)
    {
      X[k] = *(xp + k);
      Y[k] = *(yp + k);
    }

  xmin = X[0];
  xmax = X[0];
  for (k = 1; k < M; k++)
    {
      if (X[k] < xmin)
	xmin = X[k];
      if (X[k] > xmax)
	xmax = X[k];
    }

 scale = 2.0 / (xmax - xmin);
  offset = -(xmax + xmin) / (xmax - xmin);
  for (k = 0; k < M; k++)
    {
      X[k] = scale * X[k] + offset;
      /* now -1 <= X[k] <= 1 */
    }
 
  for (i = 0; i <= NN; i++)
    {
      for (j = 0; j <= NN; j++)
	{
	  a[i][j] = 0.0;
	}
    }
  /* set up a and b for Hermite polynomial fit of order N */
  for (j = 0; j <= NN; j++)
    {
      b[j] = 0.0;
      c[j] = 0.0;
      for (k = 0; k < M; k++)
	{

	  b[j] = b[j] + Y[k] * T (X[k], j);
	}
    }

  for (i = 0; i <= NN; i++)
    {
      for (j = 0; j <= NN; j++)
	{
	  for (k = 0; k < M; k++)
	    a[i][j] = a[i][j] + T (X[k], j) * T (X[k], i);
	}
    }

  b_ptr = &b[0];		/* initialise the first addresses */
  c_ptr = &c[0];
  a_ptr = &a[0];

  solve (a_ptr, b_ptr, c_ptr, NN);	/* solve the equations */


  /* create the return derivative column */
  /* insufficient data to get slope at first, last points */
  /* so don,t try                       */
  for(i=1; i<M;i++){
      *(derptr+i)=0.0;
      left=0.0;
      right=0.0;
      for(k=1;k<=NN;k++){
	  left=left+(*(c_ptr+k))*T(X[i]-change,k);
	  right=right+(*(c_ptr+k))*T(X[i]+change,k);}
      *(derptr+i)=scale*(right-left)/(2.0*change);
  }
}

void sums (double *sx, double *sy, double *sx2, double *sxy, double *x, double *y, int N){ 
    int n;
    double SX, SX2, SY, SXY;
    SX=0.0;
    SX2=0.0;
    SXY=0.0;
    SY=0.0;
    for(n=0;n<N;n++){
        SX=SX+*(x+n);
        SY=SY+*(y+n);
        SXY=SXY+(*(x+n))*(*(y+n));
        SX2=SX2+(*(x+n))*(*(x+n));
        *sx=SX;
        *sy=SY;
        *sx2=SX2;
        *sxy=SXY;
    }
}


double sigma (double x0, double v0, double a, double *x, double *y, int N){
    /* for linear fit just pass it a=0 */
    /* use to test goodness of a linear or quadratic fit */
    int n;
    double sum, term, exper;
    sum=0.0;
    for(n=0;n<N;n++){
	exper=x0+v0*(*(x+n))+0.5*a*(*(x+n))*(*(x+n));
	term=*(y+n)-exper;
	term=term*term;
	sum=sum+term;
    }
    return(sqrt(sum/(double)(N-1)));
}
	
double chi_2_prob(double val, int N){
double retval;
int m,p, num;
mpf_t w[21];
mpf_t y[21];
mpf_t D[21];
mpf_t factors[21];
mpf_t powers[21];
mpf_t chi, sum, tmp, term, two;

    for(m=0;m<=20;m++){
	mpf_init2(y[m], 256);
	mpf_init2(w[m],256);
	mpf_init2(D[m],256);
	mpf_init2(powers[m],256);
	mpf_init2(factors[m],256);
    }
    mpf_init_set_ui(sum,0);
    mpf_init_set_ui(two,2);
    mpf_init2(chi, 256);
    mpf_init2(term, 256);
    mpf_init2(sum, 256);
    mpf_init2(tmp, 256); 
 
        
 
  /* the abscissa points and weights */
 mpf_set_str(y[1],"0.705398896919887533667e-1",10);
 mpf_set_str(w[1],"0.168746801851113862149223899677e0",10);
 mpf_set_str(y[2],"0.372126818001611443794e0",10);
 mpf_set_str(w[2],"0.291254362006068281716795323806e0",10);
 mpf_set_str(y[3],"0.916582102483273564668e0",10);
 mpf_set_str(w[3],"0.266686102867001288549520868998e0",10);
 mpf_set_str(y[4],"0.170730653102834388068e1",10);
 mpf_set_str(w[4],"0.166002453269506840032140094428e0",10);
 mpf_set_str(y[5],"0.27491992553094321296e1",10);
 mpf_set_str(w[5],"0.748260646687923705422201683691e-1",10);
 mpf_set_str(y[6],"0.404892531385088692234e1",10);
 mpf_set_str(w[6],"0.249644173092832210734541468302e-1",10);
 mpf_set_str(y[7],"0.561517497086161651408e1",10);
 mpf_set_str(w[7],"0.620255084457223684757314111118e-2",10);
 mpf_set_str(y[8],"0.745901745367106330976e1",10);
 mpf_set_str(w[8],"0.114496238647690824204551915748e-2",10);
 mpf_set_str(y[9],"0.959439286958109677244e1",10);
 mpf_set_str(w[9],"0.155741773027811974784485571297e-3",10);
 mpf_set_str(y[10],"0.120388025469643163096e2",10);
 mpf_set_str(w[10],"0.15401440865224915690105126352e-4",10);
 mpf_set_str(y[11],"0.148142934426307399785e2",10);
 mpf_set_str(w[11],"0.108648636651798235149891019091e-5",10);
 mpf_set_str(y[12],"0.179488955205193760174e2",10);
 mpf_set_str(w[12],"0.533012090955671475092792098967e-7",10);
 mpf_set_str(y[13],"0.214787882402850109757e2",10);
 mpf_set_str(w[13],"0.17579811790505820036686240782e-8",10);
 mpf_set_str(y[14],"0.254517027931869055035e2",10);
 mpf_set_str(w[14],"0.372550240251232087264898344503e-10",10);
 mpf_set_str(y[15],"0.299325546317006120067e2",10);
 mpf_set_str(w[15],"0.476752925157819052455714099041e-12",10);
 mpf_set_str(y[16],"0.350134342404790000062e2",10);
 mpf_set_str(w[16],"0.337284424336243841253294925521e-14",10);
 mpf_set_str(y[17],"0.40833057056728571062e2",10);
 mpf_set_str(w[17],"0.115501433950039883096855598093e-16",10);
 mpf_set_str(y[18],"0.476199940473465021399e2",10);
 mpf_set_str(w[18],"0.153952214058234355350101868104e-19",10);
 mpf_set_str(y[19],"0.558107957500638988907e2",10);
 mpf_set_str(w[19],"0.528644272556915782906500664212e-23",10);
 mpf_set_str(y[20],"0.665244165256157538186e2",10);
 mpf_set_str(w[20],"0.165645661249902329596426420137e-27",10);
 /* denominators as function of N */
 mpf_set_str(D[2],"0.1e1", 10);
 mpf_set_str(D[3],"0.125331413731550025120788264241e1",10);
 mpf_set_str(D[4],"0.2e1",10);
 mpf_set_str(D[5],"0.375994241194650075362364792722e1",10); 
 mpf_set_str(D[6],"0.8e1",10);
 mpf_set_str(D[7],"0.187997120597325037681182396361e2",10); 
 mpf_set_str(D[8],"0.48e2",10);
 mpf_set_str(D[9],"0.131597984418127526376827677453e3",10);
 mpf_set_str(D[10],"0.384e3",10);
 mpf_set_str(D[11],"0.118438185976314773739144909707e4",10);


 mpf_set_d(chi, val);

 for(m=1;m<=20;m++){
     mpf_mul(factors[m], y[m], two);
     mpf_add(factors[m], factors[m],chi);
     mpf_sqrt(factors[m], factors[m]);
     mpf_pow_ui(powers[m], factors[m], N-2);
    
 }


 /* denominator=pow(2.0, N/2.0-1.0)*exp(gamma(N/2.0)); */
 /* compute exp(-chi/2) */
 /* need to use at least 1000000 to get 6 decimal place accuracy */
 mpf_div_ui(tmp, chi, 1000000);
 mpf_sub(tmp, two, tmp);
 mpf_div(tmp, tmp, two);
 mpf_pow_ui(tmp, tmp, 1000000);
  /* a=chi^2 */

for(m=1;m<=20;m++){
    mpf_mul(term, w[m], powers[m]);
    mpf_add(sum,sum,term);
}
 mpf_div(sum, sum, D[N]); 
 /* mpf_set_d(tmp, exp(-val/2.0)); */
 
 
 
 mpf_mul(sum, sum, tmp);
 /* mpf_out_str(stdout,10, 30, sum); */

 /* gmp_printf("Prob(Chi^2 >= %f)=%.*Ff <p>\r\n", val, 5, sum); */
 retval=mpf_get_d(sum);


 return(retval);
 for(m=0;m<21;m++){
	mpf_clear(y[m]);
	mpf_clear(w[m]);
	mpf_clear(D[m]);
	mpf_clear(powers[m]);
	mpf_clear(factors[m]);
	mpf_clear(two);
	mpf_clear(sum);
	mpf_clear(tmp);
	mpf_clear(term);
	mpf_clear(chi);

 }
}


double erf(double z){
    double a[7], power[7], sum;
    int n;
power[0]=1.0;
    for(n=1;n<7;n++){
	power[n]=z*power[n-1];
    }
    a[0]=1.0;
    a[1]=0.0705230784;
    a[2]=0.0422820123;
    a[3]=0.0092705272;
    a[4]=0.0001520143;
    a[5]=0.0002765672;
    a[6]=0.0000430638;
    sum=0.0;
    for(n=0;n<7;n++){
	sum=sum+power[n]*a[n];
    }
    sum=sum*sum;
    sum=sum*sum;
    sum=sum*sum;
    sum=sum*sum;
    return(1.0-1.0/sum);
}



double poisson(double mean, int M){
    /* poisson distribution */
    int fact[M+1];
    double powers[M+1];
    int n;
    fact[0]=1;
    powers[0]=1;
    for(n=1;n<=M;n++){
	fact[n]=n*fact[n-1];
	powers[n]=mean*powers[n-1];
    }
    return(exp(-mean)*(powers[M]/(double)fact[M]));
}

double normal_dist(double xmin, double xmax, double mean, double sig){
    double z, val1, val2;
    z=(xmax-mean)/sig;
    if(z >= 0.0)
	val1=0.5*(1.0+erf(z));
    else
	val1=0.5*(1.0-erf(fabs(z)));
    
    z=(xmin-mean)/sig;
    if(z >= 0.0)
	val2=0.5*(1.0+erf(z));
    else
	val2=0.5*(1.0-erf(fabs(z)));

    return(val1-val2);
}