/*============================================================================

 Function      ncdf

 Usage         ncdf(q)

 Arguments     q: constant double representing a quantile of the normal
                  density.
                   
 Description   This is a CUDA device function that computes the normal CDF at
               q.

 Return value  double, representing the value of the normal CDF.

 =============================================================================

 Author:       Eric M. Aldrich

 Contact:      ealdrich@gmail.com

 Date:         20 April 2010

 ============================================================================*/

// normal cdf
__device__ double ncdf(const double q)
{

 // quadrature nodes
  const double x1 = 0.9739065285;
  const double x2 = -0.9739065285;
  const double x3 = 0.8650633666;
  const double x4 = -0.8650633666;
  const double x5 = 0.6794095682;
  const double x6 = -0.6794095682;
  const double x7 = 0.4333953941;
  const double x8 = -0.4333953941;
  const double x9 = 0.1488743389;
  const double x10 = -0.1488743389;

  // quadrature weights
  const double w1 = 0.06667134430;
  const double w2 = 0.06667134430;
  const double w3 = 0.1494513491;
  const double w4 = 0.1494513491;
  const double w5 = 0.2190863625;
  const double w6 = 0.2190863625;
  const double w7 = 0.2692667193;
  const double w8 = 0.2692667193;
  const double w9 = 0.2955242247;
  const double w10 = 0.2955242247;

  // symmetry
  double c;
  if(q < 0){
    c = -q;
  } else {
    c = q;
  }

  double integral;
  // bound
  if(q >= 5.2){
    integral = 1.0;
  } else if(q <= -5.2) {
    integral = 0.0;
  } else {

    // function values
    const double pi = 3.14159265;
    const double mult = 1/pow(2*pi, 0.5);
    const double f1 = mult*exp(-0.5*pow(0.5*(x1+1)*c, 2));
    const double f2 = mult*exp(-0.5*pow(0.5*(x2+1)*c, 2));
    const double f3 = mult*exp(-0.5*pow(0.5*(x3+1)*c, 2));
    const double f4 = mult*exp(-0.5*pow(0.5*(x4+1)*c, 2));
    const double f5 = mult*exp(-0.5*pow(0.5*(x5+1)*c, 2));
    const double f6 = mult*exp(-0.5*pow(0.5*(x6+1)*c, 2));
    const double f7 = mult*exp(-0.5*pow(0.5*(x7+1)*c, 2));
    const double f8 = mult*exp(-0.5*pow(0.5*(x8+1)*c, 2));
    const double f9 = mult*exp(-0.5*pow(0.5*(x9+1)*c, 2));
    const double f10 = mult*exp(-0.5*pow(0.5*(x10+1)*c, 2));

    if(q < 0){
      integral = 0.5 - 0.5*c*(w1*f1+w2*f2+w3*f3+w4*f4+w5*f5+w6*f6+w7*f7+w8*f8+w9*f9+w10*f10);
    } else {
      integral = 0.5 + 0.5*c*(w1*f1+w2*f2+w3*f3+w4*f4+w5*f5+w6*f6+w7*f7+w8*f8+w9*f9+w10*f10);
    }

    // enforce bounds
    if(integral > 1) integral = 1;
    if(integral < 0) integral = 0;
  }

  return integral;
}
