# include <cmath>
# include <ctime>
# include <cstdlib>
# include <iostream>
# include <iomanip>

using namespace std;

# include "toms462.hpp"

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

void bivariate_normal_cdf_values ( int &n_data, double &x, double &y,
  double &r, double &fxy )

//****************************************************************************80
//
//  Purpose:
//
//    BIVARIATE_NORMAL_CDF_VALUES returns some values of the bivariate normal CDF.
//
//  Discussion:
//
//    FXY is the probability that two variables A and B, which are
//    related by a bivariate normal distribution with correlation R,
//    respectively satisfy A <= X and B <= Y.
//
//    Mathematica can evaluate the bivariate normal CDF via the commands:
//
//      <<MultivariateStatistics`
//      cdf = CDF[MultinormalDistribution[{0,0}{{1,r},{r,1}}],{x,y}]
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 November 2010
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    National Bureau of Standards,
//    Tables of the Bivariate Normal Distribution and Related Functions,
//    NBS, Applied Mathematics Series, Number 50, 1959.
//
//  Parameters:
//
//    Input/output, int &N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, double &X, &Y, the parameters of the function.
//
//    Output, double &R, the correlation value.
//
//    Output, double &FXY, the value of the function.
//
{
#  define N_MAX 41

  static double fxy_vec[N_MAX] = {
  0.02260327218569867E+00,
  0.1548729518584100E+00,
  0.4687428083352184E+00,
  0.7452035868929476E+00,
  0.8318608306874188E+00,
  0.8410314261134202E+00,
  0.1377019384919464E+00,
  0.1621749501739030E+00,
  0.1827411243233119E+00,
  0.2010067421506235E+00,
  0.2177751155265290E+00,
  0.2335088436446962E+00,
  0.2485057781834286E+00,
  0.2629747825154868E+00,
  0.2770729823404738E+00,
  0.2909261168683812E+00,
  0.3046406378726738E+00,
  0.3183113449213638E+00,
  0.3320262544108028E+00,
  0.3458686754647614E+00,
  0.3599150462310668E+00,
  0.3742210899871168E+00,
  0.3887706405282320E+00,
  0.4032765198361344E+00,
  0.4162100291953678E+00,
  0.6508271498838664E+00,
  0.8318608306874188E+00,
  0.0000000000000000,
  0.1666666666539970,
  0.2500000000000000,
  0.3333333333328906,
  0.5000000000000000,
  0.7452035868929476,
  0.1548729518584100,
  0.1548729518584100,
  0.06251409470431653,
  0.7452035868929476,
  0.1548729518584100,
  0.1548729518584100,
  0.06251409470431653,
  0.6337020457912916 };
  static double r_vec[N_MAX] = {
     0.500,  0.500,  0.500,  0.500,  0.500,
     0.500, -0.900, -0.800, -0.700, -0.600,
    -0.500, -0.400, -0.300, -0.200, -0.100,
     0.000,  0.100,  0.200,  0.300,  0.400,
     0.500,  0.600,  0.700,  0.800,  0.900,
     0.673,  0.500, -1.000, -0.500,  0.000,
     0.500,  1.000,  0.500,  0.500,  0.500,
     0.500,  0.500,  0.500,  0.500,  0.500,
     0.500 };
  static double x_vec[N_MAX] = {
    -2.0, -1.0,  0.0,  1.0,  2.0,
     3.0, -0.2, -0.2, -0.2, -0.2,
    -0.2, -0.2, -0.2, -0.2, -0.2,
    -0.2, -0.2, -0.2, -0.2, -0.2,
    -0.2, -0.2, -0.2, -0.2, -0.2,
     1.0,  2.0,  0.0,  0.0,  0.0,
     0.0,  0.0,  1.0,  1.0, -1.0,
    -1.0,  1.0,  1.0, -1.0, -1.0,
     0.7071067811865475 };
  static double y_vec[N_MAX] = {
     1.0,  1.0,  1.0,  1.0,  1.0,
     1.0,  0.5,  0.5,  0.5,  0.5,
     0.5,  0.5,  0.5,  0.5,  0.5,
     0.5,  0.5,  0.5,  0.5,  0.5,
     0.5,  0.5,  0.5,  0.5,  0.5,
     0.5,  1.0,  0.0,  0.0,  0.0,
     0.0,  0.0,  1.0, -1.0,  1.0,
    -1.0,  1.0, -1.0,  1.0, -1.0,
     0.7071067811865475 };

  if ( n_data < 0 )
  {
    n_data = 0;
  }

  n_data = n_data + 1;

  if ( N_MAX < n_data )
  {
    n_data = 0;
    r = 0.0;
    x = 0.0;
    y = 0.0;
    fxy = 0.0;
  }
  else
  {
    r = r_vec[n_data-1];
    x = x_vec[n_data-1];
    y = y_vec[n_data-1];
    fxy = fxy_vec[n_data-1];
  }

  return;
# undef N_MAX
}
//****************************************************************************80

double bivnor ( double ah, double ak, double r )

//****************************************************************************80
//
//  Purpose:
//
//    BIVNOR computes the bivariate normal CDF.
//
//  Discussion:
//
//    BIVNOR computes the probability for two normal variates X and Y
//    whose correlation is R, that AH <= X and AK <= Y.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 April 2012
//
//  Author:
//
//    Original FORTRAN77 version by Thomas Donnelly.
//    C++ version by John Burkardt.
//
//  Reference:
//
//    Thomas Donnelly,
//    Algorithm 462: Bivariate Normal Distribution,
//    Communications of the ACM,
//    October 1973, Volume 16, Number 10, page 638.
//
//  Parameters:
//
//    Input, double AH, AK, the lower limits of integration.
//
//    Input, double R, the correlation between X and Y.
//
//    Output, double BIVNOR, the bivariate normal CDF.
//
//  Local Parameters:
//
//    Local, int IDIG, the number of significant digits
//    to the right of the decimal point desired in the answer.
//
{
  double a2;
  double ap;
  double b;
  double cn;
  double con;
  double conex;
  double ex;
  double g2;
  double gh;
  double gk;
  double gw;
  double h2;
  double h4;
  int i;
  static int idig = 15;
  int is;
  double rr;
  double s1;
  double s2;
  double sgn;
  double sn;
  double sp;
  double sqr;
  double t;
  static double twopi = 6.283185307179587;
  double w2;
  double wh;
  double wk;

  b = 0.0;

  gh = gauss ( - ah ) / 2.0;
  gk = gauss ( - ak ) / 2.0;

  if ( r == 0.0 )
  {
    b = 4.00 * gh * gk;
    b = fmax ( b, 0.0 );
    b = fmin ( b, 1.0 );
    return b;
  }

  rr = ( 1.0 + r ) * ( 1.0 - r );

  if ( rr < 0.0 )
  {
    cerr << "\n";
    cerr << "BIVNOR - Fatal error!\n";
    cerr << "  1 < |R|.\n";
    exit ( 0 );
  }

  if ( rr == 0.0 )
  {
    if ( r < 0.0 )
    {
      if ( ah + ak < 0.0 )
      {
        b = 2.0 * ( gh + gk ) - 1.0;
      }
    }
    else
    {
      if ( ah - ak < 0.0 )
      {
        b = 2.0 * gk;
      }
      else
      {
        b = 2.0 * gh;
      }
    }
    b = fmax ( b, 0.0 );
    b = fmin ( b, 1.0 );
    return b;
  }

  sqr = sqrt ( rr );

  if ( idig == 15 )
  {
    con = twopi * 1.0E-15 / 2.0;
  }
  else
  {
    con = twopi / 2.0;
    for ( i = 1; i <= idig; i++ )
    {
      con = con / 10.0;
    }
  }
//
//  (0,0)
//
  if ( ah == 0.0 && ak == 0.0 )
  {
    b = 0.25 + asin ( r ) / twopi;
    b = fmax ( b, 0.0 );
    b = fmin ( b, 1.0 );
    return b;
  }
//
//  (0,nonzero)
//
  if ( ah == 0.0 && ak != 0.0 )
  {
    b = gk;
    wh = -ak;
    wk = ( ah / ak - r ) / sqr;
    gw = 2.0 * gk;
    is = 1;
  }
//
//  (nonzero,0)
//
  else if ( ah != 0.0 && ak == 0.0 )
  {
    b = gh;
    wh = -ah;
    wk = ( ak / ah - r ) / sqr;
    gw = 2.0 * gh;
    is = -1;
  }
//
//  (nonzero,nonzero)
//
  else if ( ah != 0.0 && ak != 0.0 )
  {
    b = gh + gk;
    if ( ah * ak < 0.0 )
    {
      b = b - 0.5;
    }
    wh = - ah;
    wk = ( ak / ah - r ) / sqr;
    gw = 2.0 * gh;
    is = -1;
  }

  for ( ; ; )
  {
    sgn = -1.0;
    t = 0.0;

    if ( wk != 0.0 )
    {
      if ( fabs ( wk ) == 1.0 )
      {
        t = wk * gw * ( 1.0 - gw ) / 2.0;
        b = b + sgn * t;
      }
      else
      {
        if ( 1.0 < fabs ( wk ) )
        {
          sgn = -sgn;
          wh = wh * wk;
          g2 = gauss ( wh );
          wk = 1.0 / wk;

          if ( wk < 0.0 )
          {
            b = b + 0.5;
          }
          b = b - ( gw + g2 ) / 2.0 + gw * g2;
        }
        h2 = wh * wh;
        a2 = wk * wk;
        h4 = h2 / 2.0;
        ex = exp ( - h4 );
        w2 = h4 * ex;
        ap = 1.0;
        s2 = ap - ex;
        sp = ap;
        s1 = 0.0;
        sn = s1;
        conex = fabs ( con / wk );

        for ( ; ; )
        {
          cn = ap * s2 / ( sn + sp );
          s1 = s1 + cn;

          if ( fabs ( cn ) <= conex )
          {
            break;
          }
          sn = sp;
          sp = sp + 1.0;
          s2 = s2 - w2;
          w2 = w2 * h4 / sp;
          ap = - ap * a2;
        }
        t = ( atan ( wk ) - wk * s1 ) / twopi;
        b = b + sgn * t;
      }
    }
    if ( 0 <= is )
    {
      break;
    }
    if ( ak == 0.0 )
    {
      break;
    }
    wh = -ak;
    wk = ( ah / ak - r ) / sqr;
    gw = 2.0 * gk;
    is = 1;
  }

  b = fmax ( b, 0.0 );
  b = fmin ( b, 1.0 );

  return b;
}
//****************************************************************************80

double gauss ( double t )

//****************************************************************************80
//
//  Purpose:
//
//    GAUSS returns the area of the lower tail of the normal curve.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    13 April 2012
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double T, the evaluation point.
//
//    Output, double GAUSS, the lower normal tail area.
//
{
  double value;

  value = ( 1.0 + erf ( t / sqrt ( 2.0 ) ) ) / 2.0;

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

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 September 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}