# include <cmath>
# include <cstdlib>

using namespace std;

# include "gaussian.hpp"

//****************************************************************************80

double gaussian_antideriv ( double mu, double sigma, double x )

//****************************************************************************80
//
//  Purpose:
//
//    gaussian_antideriv() evaluates the antiderivative of a Gaussian function.
//
//  Discussion:
//
//    g_anti(mu,sigmax) = 1/sigma/sqrt(2*pi) * 
//      integral ( 0 <= z <= x ) exp ( - 0.5 * (( z - mu )/sigma)^2 ) dz
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 July 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    real MU, the mean.
//
//    real SIGMA: the standard deviation.
//
//    real X: the evaluation point.
//
//  Output:
//
//    real gaussian_antideriv: values of the antiderivative function.
//
{
  double value;
  double z1;
  double z2;

  z1 =   mu       / sqrt ( 2.0 ) / sigma;
  z2 = ( mu - x ) / sqrt ( 2.0 ) / sigma;

  value = r8_erf ( z1 ) - r8_erf ( z2 );
  
  return value;
}
//****************************************************************************80

double gaussian_value ( int n, double mu, double sigma, double x )

//****************************************************************************80
//
//  Purpose:
//
//    gaussian_value() evaluates a Gaussian function or a derivative.
//
//  Discussion:
//
//    g(mu,sigmax) = 1/sigma/sqrt(2*pi) * exp ( - 0.5 * (( x - mu )/sigma).^2 )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 July 2025
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    integer N, the derivative order.
//    0 <= N.
//
//    real MU, the mean.
//
//    real SIGMA: the standard deviation.
//
//    real X: the evaluation point.
//
//  Output:
//
//    real gaussian_value: the value of the N-th derivative of the 
//    Gaussian function.
//
{
  double g;
  double hn;
  double root2;
  double value;

  root2 = sqrt ( 2.0 );

  g = exp ( - 0.5 * pow ( ( x - mu ) / sigma, 2 ) );

  hn = hermite_polynomial_value ( n, ( x - mu ) / sqrt ( sigma ) / root2 );

  value = pow ( -1.0 / sigma / root2, n ) * hn * g;
  
  return value;
}
//****************************************************************************80

double hermite_polynomial_value ( int n, double x )

//****************************************************************************80
//
//  Purpose:
//
//    hermite_polynomial_value() evaluates He(i,x).
//
//  Discussion:
//
//    He(i,x) represents the probabilist's Hermite polynomial.
//
//  Differential equation:
//
//    ( exp ( - 0.5 * x^2 ) * He(n,x)' )' + n * exp ( - 0.5 * x^2 ) * He(n,x) = 0
//
//  First terms:
//
//   1
//   X
//   X^2  -  1
//   X^3  -  3 X
//   X^4  -  6 X^2 +   3
//   X^5  - 10 X^3 +  15 X
//   X^6  - 15 X^4 +  45 X^2 -   15
//   X^7  - 21 X^5 + 105 X^3 -  105 X
//   X^8  - 28 X^6 + 210 X^4 -  420 X^2 +  105
//   X^9  - 36 X^7 + 378 X^5 - 1260 X^3 +  945 X
//   X^10 - 45 X^8 + 630 X^6 - 3150 X^4 + 4725 X^2 - 945
//
//  Recursion:
//
//    He(0,X) = 1,
//    He(1,X) = X,
//    He(N,X) = X * He(N-1,X) - (N-1) * He(N-2,X)
//
//  Orthogonality:
//
//    Integral ( -oo < X < +oo ) exp ( - 0.5 * X^2 ) * He(M,X) He(N,X) dX
//    = sqrt ( 2 * pi ) * N! * delta ( N, M )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 February 2012
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//    Frank Olver, Daniel Lozier, Ronald Boisvert, Charles Clark,
//    NIST Handbook of Mathematical Functions,
//    Cambridge University Press, 2010,
//    ISBN: 978-0521192255,
//    LC: QA331.N57.
//
//  Input:
//
//    int N, the highest order polynomial to compute.
//    Note that polynomials 0 through N will be computed.
//
//    double X, the evaluation point.
//
//  Output:
//
//    double HE_POLYNOMIAL_VALUE: the value of the 
//    probabilist's Hermite polynomials of index N.
//
{
  int j;
  double *p;
  double value;

  p = new double[n+1];

  p[0] = 1.0;

  if ( n == 0 )
  {
    value = p[0];
    delete [] p;
    return value;
  }

  p[1] = x;

  for ( j = 2; j <= n; j++ )
  {
    p[j] = x * p[j-1] - ( double ) ( j - 1 ) * p[j-2];
  }

  value = p[n];

  delete [] p;

  return value;
}
//****************************************************************************80

double r8_erf ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    r8_erf() evaluates the error function ERF(X).
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    24 May 2007
//
//  Author:
//
//    Original FORTRAN77 version by William Cody.
//    This version by John Burkardt.
//
//  Reference:
//
//    William Cody,
//    "Rational Chebyshev approximations for the error function",
//    Mathematics of Computation, 
//    1969, pages 631-638.
//
//  Input:
//
//    double X, the argument of the error function.
//
//  Output:
//
//    double R8_ERF, the value of the error function.
//
{
  double a[5] = {
    3.16112374387056560, 
    1.13864154151050156E+02, 
    3.77485237685302021E+02, 
    3.20937758913846947E+03, 
    1.85777706184603153E-01 };
  double b[4] = {
    2.36012909523441209E+01, 
    2.44024637934444173E+02, 
    1.28261652607737228E+03, 
    2.84423683343917062E+03 };
  double c[9] = {
    5.64188496988670089E-01, 
    8.88314979438837594, 
    6.61191906371416295E+01, 
    2.98635138197400131E+02, 
    8.81952221241769090E+02, 
    1.71204761263407058E+03, 
    2.05107837782607147E+03, 
    1.23033935479799725E+03, 
    2.15311535474403846E-08 };
  double d[8] = {
    1.57449261107098347E+01, 
    1.17693950891312499E+02, 
    5.37181101862009858E+02, 
    1.62138957456669019E+03, 
    3.29079923573345963E+03, 
    4.36261909014324716E+03, 
    3.43936767414372164E+03, 
    1.23033935480374942E+03 };
  double del;
  double erfx;
  int i;
  double p[6] = {
    3.05326634961232344E-01, 
    3.60344899949804439E-01, 
    1.25781726111229246E-01, 
    1.60837851487422766E-02, 
    6.58749161529837803E-04, 
    1.63153871373020978E-02 };
  double q[5] = {
    2.56852019228982242, 
    1.87295284992346047, 
    5.27905102951428412E-01, 
    6.05183413124413191E-02, 
    2.33520497626869185E-03 };
  double sqrpi = 0.56418958354775628695;
  double thresh = 0.46875;
  double xbig = 26.543;
  double xabs;
  double xden;
  double xnum;
  double xsmall = 1.11E-16;
  double xsq;

  xabs = fabs ( x );
/*
  Evaluate ERF(X) for |X| <= 0.46875.
*/
  if ( xabs <= thresh )
  {
    if ( xsmall < xabs )
    {
      xsq = xabs * xabs;
    }
    else
    {
      xsq = 0.0;
    }

    xnum = a[4] * xsq;
    xden = xsq;
    for ( i = 0; i < 3; i++ )
    {
      xnum = ( xnum + a[i] ) * xsq;
      xden = ( xden + b[i] ) * xsq;
    }

    erfx = x * ( xnum + a[3] ) / ( xden + b[3] );
  }
/*
  Evaluate ERFC(X) for 0.46875 <= |X| <= 4.0.
*/
  else if ( xabs <= 4.0 )
  {
    xnum = c[8] * xabs;
    xden = xabs;
    for ( i = 0; i < 7; i++ )
    {
      xnum = ( xnum + c[i] ) * xabs;
      xden = ( xden + d[i] ) * xabs;
    }

    erfx = ( xnum + c[7] ) / ( xden + d[7] );
    xsq = ( double ) ( ( int ) ( ( xabs * 16.0 ) / 16.0 ) );
    del = ( xabs - xsq ) * ( xabs + xsq );
    erfx = exp ( - xsq * xsq ) * exp ( - del ) * erfx;

    erfx = ( 0.5 - erfx ) + 0.5;

    if ( x < 0.0 )
    {
      erfx = - erfx;
    }
  }
/*
  Evaluate ERFC(X) for 4.0 < |X|.
*/
  else
  {
    if ( xbig <= xabs )
    {
      if ( 0.0 < x )
      {
        erfx = 1.0;
      }
      else
      {
        erfx = - 1.0;
      }
    }
    else
    {
      xsq = 1.0 / ( xabs * xabs );

      xnum = p[5] * xsq;
      xden = xsq;
      for ( i = 0; i < 4; i++ )
      {
        xnum = ( xnum + p[i] ) * xsq;
        xden = ( xden + q[i] ) * xsq;
      }

      erfx = xsq * ( xnum + p[4] ) / ( xden + q[4] );
      erfx = ( sqrpi - erfx ) / xabs;
      xsq = ( double ) ( ( int ) ( ( xabs * 16.0 ) / 16.0 ) );
      del = ( xabs - xsq ) * ( xabs + xsq );
      erfx = exp ( - xsq * xsq ) * exp ( - del ) * erfx;

      erfx = ( 0.5 - erfx ) + 0.5;
      if ( x < 0.0 )
      {
        erfx = - erfx;
      }

    }

  }

  return erfx;
}