/* deint.c -- integration of function by double exponential transform
 *  from Applied Numerical Methods with Software, by Nakamura.
 *
 *  (C) 2001, C. Bond. All rights reserved.
 */
#include <float.h>
#include <cmath>
#include <stdio.h>


/*'gaus_integrate' function comes from:
 http://finance.bi.no/~bernt/gcc_prog/recipes/recipes/node23.html*/ 
double gaus_integrate_old(double z)
{
  double ret=1.;
  if (z>6.0) 
    {
      cout <<"OVERFLOW!"<<endl;
      //return 1.0;//*sqrt(2*3.141592653589); //sqrt(2pi) now in get_pdfs()
    } // this guards against overflow
  else if (z<-6.0) 
  {
    cout <<"UNDERFLOW!"<<endl;
    ret=0.;
  }
  else {

  static const double b1 = 0.31938153;
  static const double b2 =-0.356563782;
  static const double b3 = 1.781477937;
  static const double b4 =-1.821255978;
  static const double b5 = 1.330274429;
  static const double p  = 0.2316419;
  static const double c2 = 0.3989423;
  
  //double a = fabs(z);
  double t = 1.0/(1.0+fabs(z)*p);
  double b = c2*exp((-z)*(z/2.0));
  double n= ((((b5*t+b4)*t+b3)*t+b2)*t+b1)*t;
  n=1.0-b*n;
  if (z<0.0) n=1.0-n;
  ret=n; //sqrt(2*3.141592653589)*n; //sqrt(2pi) now in get_pdfs()
  }
  return ret;
}

double gauss(double x)
{
  return exp(-0.5*pow(x,2));
}

double deint(double (*f)(double),int n,double a,double b)
{
    double h,ss,pi,z,exz,hcos,hsin,s,dxdz,x,w;
	int k;

	n = n/2;
    h = 5.0/n;    /* 5.0 is rough limit of K in exp(exp(K)) */

    pi = M_PI_4;
	ss = 0.0;
	for (k=-n;k<=n;k++) {
		z = h * (double)k;
		exz = exp(z);
		hcos = exz+1.0/exz;
		hsin = exz-1.0/exz;
        s = exp(pi*hsin);
		w = s + 1.0/s;
		x = (b*s+a/s)/w;     /* transformed abscissa  */
		if (x != a && x != b) {
            dxdz = 2*(b-a)*pi*hcos/(w*w); /* transformed weight */
			ss += h * f(x)*dxdz;
		}
	}
	return ss;
}

double deint2(double (*f)(double),int n,double a,double b)
{
    double h,ss,pi,z,exz,hcos,hsin,s,dxdz,x1,x2,w;
	int k;

	n = n/2;
    h = 5.0/n;    /* 5.0 is rough limit of K in exp(exp(K)) */

    pi = M_PI_4;
	ss = 0.5*f((a+b)/2);
	for (k=-n;k<0;k++) {
		z = h*(double)k;
		exz = exp(z);
		hcos = exz+1.0/exz;
		hsin = exz-1.0/exz;
        s = exp(pi*hsin);
		w = s + 1.0/s;
		dxdz = hcos/(w*w);      /* transformed weight   */  
		x1 = (b*s+a/s)/w;       /* transformed abscissa */
		x2 = (a*s+b/s)/w;       /*     "         "      */
		if (x1 != a && x1 != b) 
			ss += dxdz * f(x1);
		if (x2 != a && x2 != b)
			ss += dxdz * f(x2);
	 }
    return 2*(b-a)*pi*h*ss;
}

double deint3(double (*f)(double),int n,double a,double b)
{
    double h,ss,pi,z,exz,hcos,hsin,s,dxdz,x,w;
	int k;

	n = n/2;
    h = 5.0/n;    /* 5.0 is rough limit of K in exp(exp(K)) */

    pi = M_PI_4;
	ss = 0.0;
	for (k=-n;k<=n;k++) {
		z = h * (double)k;
		exz = exp(z);
		hcos = exz+1.0/exz;
		hsin = exz-1.0/exz;
        s = exp(pi*hsin);
		w = s + 1.0/s;
		x = (b*s+a/s)/w;        /* transformed abscissa  */
		if (x != a && x != b) {
			dxdz = hcos/(w*w);  /* transformed weight    */
			ss += f(x) * dxdz; 
		}
	}
    printf ("1st subtotal: %15.15lf\n", h*2*(b-a)*pi*ss);

	for (k=-n;k<=n;k++) {
		z = h * ((double)(k)+0.5);      /* use in-between points this pass */
		exz = exp(z);
		hcos = exz+1.0/exz;
		hsin = exz-1.0/exz;
        s = exp(pi*hsin);
		w = s + 1.0/s;
		x = (b*s+a/s)/w;        /* transformed abscissa  */
		if (x != a && x != b) {
			dxdz = hcos/(w*w);  /* transformed weight    */
			ss += f(x) * dxdz; 
		}
	}
     printf ("2nd subtotal: %15.15lf\n", h*(b-a)*pi*ss);

    return h*(b-a)*pi*ss;
}
/* Progressive version */
double deint4(double (*f)(double),int n,double a,double b)
{
    double h,ss,ss1,ss2,pi,z,exz,hcos,hsin,s,dxdz,x1,x2,w;
	int k,l;

	n = n/2;
    h = 5.0/n;    /* 5.0 is rough limit of K in exp(exp(K)) */

    pi = M_PI_4;
	ss = 0.5*f((a+b)/2);
	ss1 = 0.0;
	ss2 = 0.0;
	for (k=-n;k<0;k++) {
		z = h*(double)k;
		exz = exp(z);
		hcos = exz+1.0/exz;
		hsin = exz-1.0/exz;
        s = exp(pi*hsin);
		w = s + 1.0/s;
		dxdz = hcos/(w*w);      /* transformed weight   */  
		x1 = (b*s+a/s)/w;       /* transformed abscissa */
		x2 = (a*s+b/s)/w;       /*     "         "      */
		if (x1 != a && x1 != b) 
			ss1 += dxdz * f(x1);
		if (x2 != a && x2 != b)
			ss2 += dxdz * f(x2);
	}
    printf("1, %15.15lf, %15.15lf\n",2*(b-a)*pi*h*ss1,2*(b-a)*pi*h*ss2);

	for (l=0;l<3;l++) {

		for (k=-n;k<0;k++) {
			z = h*((double)k+0.5);
			exz = exp(z);
			hcos = exz+1.0/exz;
			hsin = exz-1.0/exz;
            s = exp(pi*hsin);
			w = s + 1.0/s;
			dxdz = hcos/(w*w);      /* transformed weight   */ 
			x1 = (b*s+a/s)/w;       /* transformed abscissa */
			x2 = (a*s+b/s)/w;       /*     "         "      */
			if (x1 != a && x1 != b) 
				ss1 += dxdz * f(x1);
			if (x2 != a && x2 != b)
				ss2 += dxdz * f(x2);
		}
        printf("%d, %15.15lf, %15.15lf\n",l+2,(b-a)*pi*h*ss1,(b-a)*pi*h*ss2);
		n *= 2;
		h /= 2;
	}
    return 2*(b-a)*pi*h*(ss+ss1+ss2);
}

/////////////////////////////////////////////////////////////////////////////
double gaus_integrate(double z) {
  //return gaus_integrate_old(z);
  //lookup table for values calculated by gaus_integrate_old(z) for z=0 to 3 in steps of 0.002
static const float value[] = { 0.5, 0.500798, 0.501596, 0.502394, 0.503191, 0.503989, 0.504787, 0.505585, 0.506383, 0.507181,  \
0.507978, 0.508776, 0.509574, 0.510371, 0.511169, 0.511967, 0.512764, 0.513561, 0.514359, 0.515156,  \
0.515953, 0.516751, 0.517548, 0.518345, 0.519142, 0.519939, 0.520736, 0.521532, 0.522329, 0.523126,  \
0.523922, 0.524719, 0.525515, 0.526311, 0.527107, 0.527903, 0.528699, 0.529495, 0.53029, 0.531086,  \
0.531881, 0.532677, 0.533472, 0.534267, 0.535062, 0.535856, 0.536651, 0.537445, 0.53824, 0.539034,  \
0.539828, 0.540622, 0.541415, 0.542209, 0.543002, 0.543795, 0.544588, 0.545381, 0.546174, 0.546966,  \
0.547758, 0.548551, 0.549342, 0.550134, 0.550926, 0.551717, 0.552508, 0.553299, 0.554089, 0.55488,  \
0.55567, 0.55646, 0.55725, 0.558039, 0.558829, 0.559618, 0.560407, 0.561195, 0.561983, 0.562772,  \
0.563559, 0.564347, 0.565134, 0.565922, 0.566708, 0.567495, 0.568281, 0.569067, 0.569853, 0.570638,  \
0.571424, 0.572209, 0.572993, 0.573778, 0.574562, 0.575345, 0.576129, 0.576912, 0.577695, 0.578477,  \
0.57926, 0.580042, 0.580823, 0.581604, 0.582385, 0.583166, 0.583946, 0.584726, 0.585506, 0.586285,  \
0.587064, 0.587843, 0.588621, 0.589399, 0.590177, 0.590954, 0.591731, 0.592507, 0.593284, 0.594059,  \
0.594835, 0.59561, 0.596385, 0.597159, 0.597933, 0.598706, 0.599479, 0.600252, 0.601025, 0.601796,  \
0.602568, 0.603339, 0.60411, 0.60488, 0.60565, 0.60642, 0.607189, 0.607958, 0.608726, 0.609494,  \
0.610261, 0.611028, 0.611795, 0.612561, 0.613327, 0.614092, 0.614857, 0.615621, 0.616385, 0.617148,  \
0.617911, 0.618674, 0.619436, 0.620198, 0.620959, 0.621719, 0.62248, 0.623239, 0.623999, 0.624757,  \
0.625516, 0.626274, 0.627031, 0.627788, 0.628544, 0.6293, 0.630055, 0.63081, 0.631564, 0.632318,  \
0.633072, 0.633824, 0.634577, 0.635329, 0.63608, 0.636831, 0.637581, 0.63833, 0.63908, 0.639828,  \
0.640576, 0.641324, 0.642071, 0.642817, 0.643563, 0.644309, 0.645054, 0.645798, 0.646541, 0.647285,  \
0.648027, 0.648769, 0.649511, 0.650252, 0.650992, 0.651732, 0.652471, 0.653209, 0.653947, 0.654685,  \
0.655422, 0.656158, 0.656894, 0.657629, 0.658363, 0.659097, 0.65983, 0.660563, 0.661295, 0.662026,  \
0.662757, 0.663487, 0.664217, 0.664946, 0.665674, 0.666402, 0.667129, 0.667856, 0.668582, 0.669307,  \
0.670031, 0.670755, 0.671479, 0.672201, 0.672923, 0.673645, 0.674365, 0.675086, 0.675805, 0.676524,  \
0.677242, 0.677959, 0.678676, 0.679392, 0.680108, 0.680822, 0.681537, 0.68225, 0.682963, 0.683675,  \
0.684386, 0.685097, 0.685807, 0.686516, 0.687225, 0.687933, 0.68864, 0.689347, 0.690053, 0.690758,  \
0.691462, 0.692166, 0.692869, 0.693572, 0.694273, 0.694974, 0.695674, 0.696374, 0.697073, 0.697771,  \
0.698468, 0.699165, 0.699861, 0.700556, 0.70125, 0.701944, 0.702637, 0.703329, 0.704021, 0.704711,  \
0.705401, 0.706091, 0.706779, 0.707467, 0.708154, 0.70884, 0.709526, 0.710211, 0.710895, 0.711578,  \
0.71226, 0.712942, 0.713623, 0.714303, 0.714983, 0.715661, 0.716339, 0.717016, 0.717692, 0.718368,  \
0.719043, 0.719717, 0.72039, 0.721062, 0.721734, 0.722405, 0.723075, 0.723744, 0.724412, 0.72508,  \
0.725747, 0.726413, 0.727078, 0.727743, 0.728406, 0.729069, 0.729731, 0.730392, 0.731053, 0.731712,  \
0.732371, 0.733029, 0.733686, 0.734343, 0.734998, 0.735653, 0.736307, 0.73696, 0.737612, 0.738263,  \
0.738914, 0.739563, 0.740212, 0.74086, 0.741508, 0.742154, 0.742799, 0.743444, 0.744088, 0.744731,  \
0.745373, 0.746014, 0.746655, 0.747295, 0.747933, 0.748571, 0.749208, 0.749844, 0.75048, 0.751114,  \
0.751748, 0.752381, 0.753012, 0.753644, 0.754274, 0.754903, 0.755531, 0.756159, 0.756786, 0.757411,  \
0.758036, 0.75866, 0.759284, 0.759906, 0.760527, 0.761148, 0.761768, 0.762386, 0.763004, 0.763621,  \
0.764238, 0.764853, 0.765467, 0.766081, 0.766693, 0.767305, 0.767916, 0.768526, 0.769135, 0.769743,  \
0.77035, 0.770956, 0.771562, 0.772166, 0.77277, 0.773373, 0.773975, 0.774575, 0.775175, 0.775775,  \
0.776373, 0.77697, 0.777566, 0.778162, 0.778756, 0.77935, 0.779943, 0.780535, 0.781126, 0.781716,  \
0.782305, 0.782893, 0.78348, 0.784066, 0.784652, 0.785236, 0.78582, 0.786402, 0.786984, 0.787565,  \
0.788145, 0.788724, 0.789302, 0.789879, 0.790455, 0.79103, 0.791604, 0.792178, 0.79275, 0.793321,  \
0.793892, 0.794462, 0.79503, 0.795598, 0.796165, 0.796731, 0.797296, 0.79786, 0.798423, 0.798985,  \
0.799546, 0.800106, 0.800665, 0.801224, 0.801781, 0.802337, 0.802893, 0.803448, 0.804001, 0.804554,  \
0.805106, 0.805656, 0.806206, 0.806755, 0.807303, 0.80785, 0.808396, 0.808941, 0.809485, 0.810028,  \
0.81057, 0.811112, 0.811652, 0.812191, 0.81273, 0.813267, 0.813804, 0.814339, 0.814874, 0.815407,  \
0.81594, 0.816472, 0.817002, 0.817532, 0.818061, 0.818589, 0.819116, 0.819642, 0.820167, 0.820691,  \
0.821214, 0.821736, 0.822257, 0.822777, 0.823296, 0.823814, 0.824332, 0.824848, 0.825363, 0.825878,  \
0.826391, 0.826904, 0.827415, 0.827926, 0.828435, 0.828944, 0.829452, 0.829958, 0.830464, 0.830969,  \
0.831472, 0.831975, 0.832477, 0.832978, 0.833478, 0.833977, 0.834475, 0.834972, 0.835468, 0.835963,  \
0.836457, 0.83695, 0.837442, 0.837933, 0.838424, 0.838913, 0.839401, 0.839889, 0.840375, 0.84086,  \
0.841345, 0.841828, 0.842311, 0.842792, 0.843273, 0.843752, 0.844231, 0.844709, 0.845185, 0.845661,  \
0.846136, 0.84661, 0.847082, 0.847554, 0.848025, 0.848495, 0.848964, 0.849432, 0.849899, 0.850365,  \
0.85083, 0.851294, 0.851757, 0.852219, 0.852681, 0.853141, 0.8536, 0.854059, 0.854516, 0.854972,  \
0.855428, 0.855882, 0.856336, 0.856788, 0.85724, 0.85769, 0.85814, 0.858589, 0.859036, 0.859483,  \
0.859929, 0.860374, 0.860818, 0.86126, 0.861702, 0.862143, 0.862583, 0.863022, 0.863461, 0.863898,  \
0.864334, 0.864769, 0.865203, 0.865637, 0.866069, 0.8665, 0.866931, 0.86736, 0.867789, 0.868216,  \
0.868643, 0.869069, 0.869493, 0.869917, 0.87034, 0.870762, 0.871183, 0.871603, 0.872022, 0.87244,  \
0.872857, 0.873273, 0.873688, 0.874102, 0.874516, 0.874928, 0.875339, 0.87575, 0.876159, 0.876568,  \
0.876976, 0.877382, 0.877788, 0.878193, 0.878597, 0.878999, 0.879401, 0.879802, 0.880202, 0.880602,  \
0.881, 0.881397, 0.881793, 0.882189, 0.882583, 0.882977, 0.883369, 0.883761, 0.884152, 0.884541,  \
0.88493, 0.885318, 0.885705, 0.886091, 0.886476, 0.88686, 0.887244, 0.887626, 0.888007, 0.888388,  \
0.888767, 0.889146, 0.889524, 0.889901, 0.890276, 0.890651, 0.891025, 0.891398, 0.891771, 0.892142,  \
0.892512, 0.892882, 0.89325, 0.893618, 0.893984, 0.89435, 0.894715, 0.895079, 0.895442, 0.895804,  \
0.896165, 0.896526, 0.896885, 0.897243, 0.897601, 0.897958, 0.898313, 0.898668, 0.899022, 0.899375,  \
0.899727, 0.900079, 0.900429, 0.900778, 0.901127, 0.901475, 0.901821, 0.902167, 0.902512, 0.902856,  \
0.903199, 0.903542, 0.903883, 0.904224, 0.904563, 0.904902, 0.90524, 0.905577, 0.905913, 0.906248,  \
0.906582, 0.906916, 0.907248, 0.90758, 0.907911, 0.908241, 0.90857, 0.908898, 0.909225, 0.909552,  \
0.909877, 0.910202, 0.910526, 0.910849, 0.911171, 0.911492, 0.911812, 0.912132, 0.91245, 0.912768,  \
0.913085, 0.913401, 0.913716, 0.91403, 0.914344, 0.914656, 0.914968, 0.915279, 0.915589, 0.915898,  \
0.916207, 0.916514, 0.916821, 0.917126, 0.917431, 0.917736, 0.918039, 0.918341, 0.918643, 0.918943,  \
0.919243, 0.919542, 0.919841, 0.920138, 0.920434, 0.92073, 0.921025, 0.921319, 0.921612, 0.921905,  \
0.922196, 0.922487, 0.922777, 0.923066, 0.923354, 0.923641, 0.923928, 0.924214, 0.924499, 0.924783,  \
0.925066, 0.925349, 0.92563, 0.925911, 0.926191, 0.926471, 0.926749, 0.927027, 0.927304, 0.92758,  \
0.927855, 0.928129, 0.928403, 0.928676, 0.928948, 0.929219, 0.92949, 0.929759, 0.930028, 0.930296,  \
0.930563, 0.93083, 0.931096, 0.93136, 0.931625, 0.931888, 0.93215, 0.932412, 0.932673, 0.932933,  \
0.933193, 0.933451, 0.933709, 0.933966, 0.934223, 0.934478, 0.934733, 0.934987, 0.93524, 0.935493,  \
0.935744, 0.935995, 0.936246, 0.936495, 0.936744, 0.936992, 0.937239, 0.937485, 0.937731, 0.937976,  \
0.93822, 0.938463, 0.938706, 0.938948, 0.939189, 0.939429, 0.939669, 0.939908, 0.940146, 0.940383,  \
0.94062, 0.940856, 0.941091, 0.941326, 0.941559, 0.941792, 0.942025, 0.942256, 0.942487, 0.942717,  \
0.942947, 0.943175, 0.943403, 0.94363, 0.943857, 0.944083, 0.944308, 0.944532, 0.944756, 0.944979,  \
0.945201, 0.945422, 0.945643, 0.945863, 0.946082, 0.946301, 0.946519, 0.946736, 0.946953, 0.947169,  \
0.947384, 0.947598, 0.947812, 0.948025, 0.948238, 0.948449, 0.94866, 0.948871, 0.94908, 0.949289,  \
0.949497, 0.949705, 0.949912, 0.950118, 0.950324, 0.950529, 0.950733, 0.950936, 0.951139, 0.951341,  \
0.951543, 0.951744, 0.951944, 0.952143, 0.952342, 0.95254, 0.952738, 0.952935, 0.953131, 0.953326,  \
0.953521, 0.953716, 0.953909, 0.954102, 0.954294, 0.954486, 0.954677, 0.954867, 0.955057, 0.955246,  \
0.955435, 0.955622, 0.955809, 0.955996, 0.956182, 0.956367, 0.956552, 0.956736, 0.956919, 0.957102,  \
0.957284, 0.957465, 0.957646, 0.957826, 0.958006, 0.958185, 0.958363, 0.958541, 0.958718, 0.958895,  \
0.959071, 0.959246, 0.95942, 0.959595, 0.959768, 0.959941, 0.960113, 0.960285, 0.960456, 0.960626,  \
0.960796, 0.960965, 0.961134, 0.961302, 0.96147, 0.961636, 0.961803, 0.961968, 0.962134, 0.962298,  \
0.962462, 0.962625, 0.962788, 0.96295, 0.963112, 0.963273, 0.963434, 0.963593, 0.963753, 0.963912,  \
0.96407, 0.964227, 0.964384, 0.964541, 0.964697, 0.964852, 0.965007, 0.965161, 0.965315, 0.965468,  \
0.965621, 0.965773, 0.965924, 0.966075, 0.966225, 0.966375, 0.966524, 0.966673, 0.966821, 0.966969,  \
0.967116, 0.967262, 0.967408, 0.967554, 0.967699, 0.967843, 0.967987, 0.96813, 0.968273, 0.968416,  \
0.968557, 0.968699, 0.968839, 0.968979, 0.969119, 0.969258, 0.969397, 0.969535, 0.969672, 0.969809,  \
0.969946, 0.970082, 0.970218, 0.970353, 0.970487, 0.970621, 0.970755, 0.970888, 0.97102, 0.971152,  \
0.971284, 0.971414, 0.971545, 0.971675, 0.971804, 0.971933, 0.972062, 0.97219, 0.972318, 0.972445,  \
0.972571, 0.972697, 0.972823, 0.972948, 0.973073, 0.973197, 0.97332, 0.973444, 0.973566, 0.973688,  \
0.97381, 0.973932, 0.974052, 0.974173, 0.974293, 0.974412, 0.974531, 0.974649, 0.974767, 0.974885,  \
0.975002, 0.975119, 0.975235, 0.975351, 0.975466, 0.975581, 0.975695, 0.975809, 0.975923, 0.976036,  \
0.976148, 0.97626, 0.976372, 0.976483, 0.976594, 0.976705, 0.976815, 0.976924, 0.977033, 0.977142,  \
0.97725, 0.977358, 0.977465, 0.977572, 0.977678, 0.977784, 0.97789, 0.977995, 0.9781, 0.978204,  \
0.978308, 0.978412, 0.978515, 0.978618, 0.97872, 0.978822, 0.978923, 0.979024, 0.979125, 0.979225,  \
0.979325, 0.979424, 0.979523, 0.979622, 0.97972, 0.979818, 0.979915, 0.980012, 0.980109, 0.980205,  \
0.980301, 0.980396, 0.980491, 0.980586, 0.98068, 0.980774, 0.980867, 0.98096, 0.981053, 0.981145,  \
0.981237, 0.981329, 0.98142, 0.981511, 0.981601, 0.981691, 0.981781, 0.98187, 0.981959, 0.982047,  \
0.982136, 0.982223, 0.982311, 0.982398, 0.982485, 0.982571, 0.982657, 0.982742, 0.982828, 0.982913,  \
0.982997, 0.983081, 0.983165, 0.983248, 0.983332, 0.983414, 0.983497, 0.983579, 0.98366, 0.983742,  \
0.983823, 0.983903, 0.983984, 0.984064, 0.984143, 0.984222, 0.984301, 0.98438, 0.984458, 0.984536,  \
0.984614, 0.984691, 0.984768, 0.984844, 0.984921, 0.984997, 0.985072, 0.985147, 0.985222, 0.985297,  \
0.985371, 0.985445, 0.985519, 0.985592, 0.985665, 0.985738, 0.98581, 0.985882, 0.985954, 0.986025,  \
0.986097, 0.986167, 0.986238, 0.986308, 0.986378, 0.986447, 0.986517, 0.986586, 0.986654, 0.986723,  \
0.986791, 0.986858, 0.986926, 0.986993, 0.98706, 0.987126, 0.987193, 0.987259, 0.987324, 0.98739,  \
0.987455, 0.987519, 0.987584, 0.987648, 0.987712, 0.987776, 0.987839, 0.987902, 0.987965, 0.988027,  \
0.988089, 0.988151, 0.988213, 0.988274, 0.988335, 0.988396, 0.988457, 0.988517, 0.988577, 0.988637,  \
0.988696, 0.988755, 0.988814, 0.988873, 0.988931, 0.988989, 0.989047, 0.989105, 0.989162, 0.989219,  \
0.989276, 0.989332, 0.989389, 0.989445, 0.9895, 0.989556, 0.989611, 0.989666, 0.989721, 0.989775,  \
0.98983, 0.989884, 0.989937, 0.989991, 0.990044, 0.990097, 0.99015, 0.990202, 0.990254, 0.990306,  \
0.990358, 0.99041, 0.990461, 0.990512, 0.990563, 0.990613, 0.990664, 0.990714, 0.990764, 0.990813,  \
0.990863, 0.990912, 0.990961, 0.991009, 0.991058, 0.991106, 0.991154, 0.991202, 0.991249, 0.991297,  \
0.991344, 0.991391, 0.991437, 0.991484, 0.99153, 0.991576, 0.991622, 0.991667, 0.991712, 0.991758,  \
0.991802, 0.991847, 0.991892, 0.991936, 0.99198, 0.992024, 0.992067, 0.992111, 0.992154, 0.992197,  \
0.99224, 0.992282, 0.992325, 0.992367, 0.992409, 0.992451, 0.992492, 0.992534, 0.992575, 0.992616,  \
0.992656, 0.992697, 0.992737, 0.992777, 0.992817, 0.992857, 0.992897, 0.992936, 0.992975, 0.993014,  \
0.993053, 0.993092, 0.99313, 0.993168, 0.993206, 0.993244, 0.993282, 0.99332, 0.993357, 0.993394,  \
0.993431, 0.993468, 0.993504, 0.993541, 0.993577, 0.993613, 0.993649, 0.993684, 0.99372, 0.993755,  \
0.99379, 0.993825, 0.99386, 0.993895, 0.993929, 0.993963, 0.993998, 0.994031, 0.994065, 0.994099,  \
0.994132, 0.994165, 0.994199, 0.994232, 0.994264, 0.994297, 0.994329, 0.994362, 0.994394, 0.994426,  \
0.994457, 0.994489, 0.99452, 0.994552, 0.994583, 0.994614, 0.994645, 0.994675, 0.994706, 0.994736,  \
0.994766, 0.994796, 0.994826, 0.994856, 0.994886, 0.994915, 0.994944, 0.994973, 0.995002, 0.995031,  \
0.99506, 0.995088, 0.995117, 0.995145, 0.995173, 0.995201, 0.995229, 0.995257, 0.995284, 0.995312,  \
0.995339, 0.995366, 0.995393, 0.99542, 0.995446, 0.995473, 0.995499, 0.995526, 0.995552, 0.995578,  \
0.995603, 0.995629, 0.995655, 0.99568, 0.995706, 0.995731, 0.995756, 0.995781, 0.995805, 0.99583,  \
0.995855, 0.995879, 0.995903, 0.995927, 0.995951, 0.995975, 0.995999, 0.996023, 0.996046, 0.99607,  \
0.996093, 0.996116, 0.996139, 0.996162, 0.996185, 0.996207, 0.99623, 0.996252, 0.996275, 0.996297,  \
0.996319, 0.996341, 0.996363, 0.996384, 0.996406, 0.996427, 0.996449, 0.99647, 0.996491, 0.996512,  \
0.996533, 0.996554, 0.996574, 0.996595, 0.996615, 0.996636, 0.996656, 0.996676, 0.996696, 0.996716,  \
0.996736, 0.996756, 0.996775, 0.996795, 0.996814, 0.996833, 0.996852, 0.996871, 0.99689, 0.996909,  \
0.996928, 0.996947, 0.996965, 0.996984, 0.997002, 0.99702, 0.997038, 0.997056, 0.997074, 0.997092,  \
0.99711, 0.997128, 0.997145, 0.997163, 0.99718, 0.997197, 0.997214, 0.997231, 0.997248, 0.997265,  \
0.997282, 0.997299, 0.997315, 0.997332, 0.997348, 0.997365, 0.997381, 0.997397, 0.997413, 0.997429,  \
0.997445, 0.997461, 0.997476, 0.997492, 0.997507, 0.997523, 0.997538, 0.997553, 0.997569, 0.997584,  \
0.997599, 0.997614, 0.997629, 0.997643, 0.997658, 0.997673, 0.997687, 0.997701, 0.997716, 0.99773,  \
0.997744, 0.997758, 0.997772, 0.997786, 0.9978, 0.997814, 0.997828, 0.997841, 0.997855, 0.997868,  \
0.997882, 0.997895, 0.997908, 0.997921, 0.997935, 0.997948, 0.997961, 0.997973, 0.997986, 0.997999,  \
0.998012, 0.998024, 0.998037, 0.998049, 0.998061, 0.998074, 0.998086, 0.998098, 0.99811, 0.998122,  \
0.998134, 0.998146, 0.998158, 0.99817, 0.998181, 0.998193, 0.998204, 0.998216, 0.998227, 0.998239,  \
0.99825, 0.998261, 0.998272, 0.998283, 0.998294, 0.998305, 0.998316, 0.998327, 0.998338, 0.998348,  \
0.998359, 0.998369, 0.99838, 0.99839, 0.998401, 0.998411, 0.998421, 0.998432, 0.998442, 0.998452,  \
0.998462, 0.998472, 0.998482, 0.998491, 0.998501, 0.998511, 0.998521, 0.99853, 0.99854, 0.998549,  \
0.998559, 0.998568, 0.998577, 0.998587, 0.998596, 0.998605, 0.998614, 0.998623, 0.998632, 0.998641,  \
			       0.99865 };
 double ret=0.;

 if (z>=0.) {
   if (z<=3) {
     ret=value[int(500.*z)];
   }
   else {
     ret=1.;
   }
 }
 else {
   if (z>=-3) {
     ret=1.-value[int(-500.*z)];
   }
   //   else {
   //  ret=0.;
   //}
 }
 // for cross check:
 //if (fabs(ret-gaus_integrate_old(z))>0.002) {
 // cout<<"z: "<<z<<"  gaus_integrate: "<<ret<<"  gaus_integrate_old: "<<gaus_integrate_old(z)<<endl;
   //exit(1);
   //}
 return ret;
}