// Noel Belcourt
// April 26, 1999

#include <cmath>
#include <cstdlib>
#include <string>

#include <normal.hpp>

using std::atan;
using std::log;
using std::sqrt;
using std::string;

normal::normal(string s, double mean, double variance) :
    distribution(s), mu(mean), sigma(sqrt(variance))
  , c0(2.515517), c1(0.802853), c2(0.010328)
  , d1(1.432788), d2(0.189269), d3(0.001308)
  , epsilon(4.5e-4), pi(4.0*atan(1.0))
{
}

double normal::density(double x)
{
  return exp(-(x-mu)*(x-mu)/(2.0*sigma*sigma))/(sigma*sqrt(2*pi));
}

double normal::inverse(double p)
{
  double value = 0;
  if (0 <= p && p <= 1) {

    double xp;
    if (p > 0.5) {
      xp = solve(1-p);
    }
    else {
      xp = -solve(p);
    }

    // Transform from standard normal to this mean and variance
    value = mu + xp * sqrt(sigma);
  }
  return value;
}

double normal::solve(double p)
{
  if (0 < p && p <= 0.5) {
    // Transform 
    double t = sqrt(log(1/(p*p)));

    // Compute standard normal estimate for this probability
    return t - (c0 + c1*t + c2*t*t) / (1 + d1*t + d2*t*t + d3*t*t*t) + epsilon;
  }
  else {
    return 0;
  }
}