# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <time.h>

int main ( int argc, char *argv[] );
void jacobi_eigenvalue ( int n, double a[], int it_max, double v[], 
  double d[], int *it_num, int *rot_num );
void moment_method ( int n, double moment[], double x[], double w[] );
double *moments_normal ( int m, double mu, double sigma );
double *moments_truncated_normal_ab ( int m, double mu, double sigma,
  double a, double b );
double *moments_truncated_normal_a ( int m, double mu, double sigma,
  double a );
double *moments_truncated_normal_b ( int m, double mu, double sigma,
  double b );
double normal_01_cdf ( double x );
double normal_01_pdf ( double x );
double r8_choose ( int n, int k );
double r8_factorial ( int n );
double r8_factorial2 ( int n );
double r8_mop ( int i );
double *r8mat_cholesky_factor_upper ( int n, double a[], int *flag );
double *r8mat_copy_new ( int m, int n, double a1[] );
void r8mat_diag_get_vector ( int n, double a[], double v[] );
void r8mat_identity  ( int n, double a[] );
void r8mat_print ( int m, int n, double a[], char *title );
void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, char *title );
void r8mat_write ( char *output_filename, int m, int n, double table[] );
void r8vec_print ( int n, double a[], char *title );
void rule_write ( int order, char *filename, double x[], double w[],
  double r[] );
void timestamp ( );
double truncated_normal_ab_moment ( int order, double mu, double s, double a,
  double b );
double truncated_normal_a_moment ( int order, double mu, double s, double a );
double truncated_normal_b_moment ( int order, double mu, double s, double b );

/******************************************************************************/

int main ( int argc, char *argv[] )

/******************************************************************************/
/*
  Purpose:

    MAIN is the main program for TRUNCATED_NORMAL_RULE.

  Discussion:

    This program computes a truncated normal quadrature rule
    and writes it to a file.

    The user specifies:
    * option: 0/1/2/3 for none, lower, upper, double truncation.
    * N, the number of points in the rule;
    * MU, the mean of the original normal distribution;
    * SIGMA, the standard deviation of the original normal distribution,
    * A, the left endpoint (for options 1 or 3)
    * B, the right endpoint (for options 2 or 3);
    * FILENAME, the root name of the output files.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    20 September 2013

  Author:

    John Burkardt
*/
{
  double a;
  double b;
  char filename[255];
  int iarg;
  double *moment;
  double mu;
  int n;
  int option;
  double *r;
  double sigma;
  double *w;
  double *x;

  timestamp ( );
  printf ( "\n" );
  printf ( "TRUNCATED_NORMAL_RULE\n" );
  printf ( "  C version\n" );
  printf ( "\n" );
  printf ( "  For the (truncated) Gaussian probability density function\n" );
  printf ( "    pdf(x) = exp(-0.5*((x-MU)/SIGMA)^2) / SIGMA / sqrt ( 2 * pi )\n" );
  printf ( "  compute an N-point quadrature rule for approximating\n" );
  printf ( "    Integral ( A <= x <= B ) f(x) pdf(x) dx\n" );
  printf ( "\n" );
  printf ( "  The value of OPTION determines the truncation interval [A,B]:\n" );
  printf ( "  0: (-oo,+oo)\n" );
  printf ( "  1: [A,+oo)\n" );
  printf ( "  2: (-oo,B]\n" );
  printf ( "  3: [A,B]\n" );
  printf ( "\n" );
  printf ( "  The user specifies OPTION, N, MU, SIGMA, A, B and FILENAME.\n" );
  printf ( "\n" );
  printf ( "  FILENAME is used to generate 3 files:\n" );
  printf ( "    filename_w.txt - the weight file\n" );
  printf ( "    filename_x.txt - the abscissa file.\n" );
  printf ( "    filename_r.txt - the region file, listing A and B.\n" );

  iarg = 0;
/*
  Get OPTION.
*/
  iarg = iarg + 1;
  if ( iarg < argc )
  {
    option = atoi ( argv[iarg] );
  }
  else
  {
    printf ( "\n" );
    printf ( "  Enter the value of OPTION 0/1/2/3:  " );
    scanf ( "%d", &option );
  }
  if ( option < 0 || 3 < option )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_RULE - Fatal error!\n" );
    fprintf ( stderr, "  0 <= OPTION <= 3 was required.\n" );
    exit ( 1 );
  }
/*
  Get N.
*/
  iarg = iarg + 1;
  if ( iarg < argc )
  {
    n = atoi ( argv[iarg] );
  }
  else
  {
    printf ( "\n" );
    printf ( "  Enter the value of N (1 or greater)\n" );
    scanf ( "%d", &n );
  }
/*
  Get MU.
*/
  iarg = iarg + 1;
  if ( iarg < argc )
  {
    mu = atof ( argv[iarg] );
  }
  else
  {
    printf ( "\n" );
    printf ( "  Enter MU, the mean value of the normal distribution:\n" );
    scanf ( "%lf", &mu );
  }
/*
  Get SIGMA.
*/
  iarg = iarg + 1;
  if ( iarg < argc )
  {
    sigma = atof ( argv[iarg] );
  }
  else
  {
    printf ( "\n" );
    printf ( "  Enter SIGMA, the standard deviation of the normal distribution:\n" );
    scanf ( "%lf", &sigma );
  }
  sigma = fabs ( sigma );
/*
  Get A.
*/
  if ( option == 1 || option == 3 )
  {
    iarg = iarg + 1;
    if ( iarg < argc )
    {
      a = atof ( argv[iarg] );
    }
    else
    {
      printf ( "\n" );
      printf ( "  Enter the left endpoint A:\n" );
      scanf ( "%lf", &a );
    }
  }
  else
  {
    a = - HUGE_VAL;
  }
/*
  Get B.
*/
  if ( option == 2 || option == 3 )
  {
    iarg = iarg + 1;
    if ( iarg < argc )
    {
      b = atof ( argv[iarg] );
    }
    else
    {
      printf ( "\n" );
      printf ( "  Enter the right endpoint B:\n" );
      scanf ( "%lf", &b );
    }
  }
  else
  {
    b = HUGE_VAL;
  }

  if ( b <= a )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_RULE - Fatal error!\n" );
    fprintf ( stderr, "  A < B required!\n" );
    exit ( 1 );
  }
/*
  Get FILENAME:
*/
  iarg = iarg + 1;
  if ( iarg < argc )
  {
    strcpy ( filename, argv[iarg] );
  }
  else
  {
    printf ( "\n" );
    printf ( "  Enter FILENAME, the \"root name\" of the quadrature files).\n" );
    scanf ( "%s", filename );
  }
/*
  Input summary.
*/
  printf ( "\n" );
  printf ( "  OPTION = %d\n", option );
  printf ( "  N = %d\n", n );
  printf ( "  MU = %g\n", mu );
  printf ( "  SIGMA = %g\n", sigma );
  printf ( "  A = %g\n", a );
  printf ( "  B = %g\n", b );
  printf ( "  FILENAME = \"%s\".\n", filename );
/*
  Compute the moments.
*/
  if ( option == 0 )
  {
    moment = moments_normal ( 2 * n + 1, mu, sigma );
  }
  else if ( option == 1 )
  {
    moment = moments_truncated_normal_a ( 2 * n + 1, mu, sigma, a );
  }
  else if ( option == 2 )
  {
    moment = moments_truncated_normal_b ( 2 * n + 1, mu, sigma, b );
  }
  else if ( option == 3 )
  {
    moment = moments_truncated_normal_ab ( 2 * n + 1, mu, sigma, a, b );
  }
/*
  Construct the rule from the moments.
*/
  w = ( double * ) malloc ( n * sizeof ( double ) );
  x = ( double * ) malloc ( n * sizeof ( double ) );

  moment_method ( n, moment, x, w );
/*
  Output the rule.
*/
  r = ( double * ) malloc ( 2 * sizeof ( double ) );
  r[0] = a;
  r[1] = b;

  rule_write ( n, filename, x, w, r );
/*
  Free memory.
*/
  free ( moment );
  free ( r );
  free ( w );
  free ( x );
/*
  Terminate.
*/
  printf ( "\n" );
  printf ( "TRUNCATED_NORMAL_RULE:\n" );
  printf ( "  Normal end of execution.\n" );
  printf ( "\n" );
  timestamp ( );

  return 0;
}
/******************************************************************************/

void jacobi_eigenvalue ( int n, double a[], int it_max, double v[], 
  double d[], int *it_num, int *rot_num )

/******************************************************************************/
/*
  Purpose:

    JACOBI_EIGENVALUE carries out the Jacobi eigenvalue iteration.

  Discussion:

    This function computes the eigenvalues and eigenvectors of a
    real symmetric matrix, using Rutishauser's modfications of the classical
    Jacobi rotation method with threshold pivoting. 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    17 September 2013

  Author:

    C version by John Burkardt

  Parameters:

    Input, int N, the order of the matrix.

    Input, double A[N*N], the matrix, which must be square, real,
    and symmetric.

    Input, int IT_MAX, the maximum number of iterations.

    Output, double V[N*N], the matrix of eigenvectors.

    Output, double D[N], the eigenvalues, in descending order.

    Output, int *IT_NUM, the total number of iterations.

    Output, int *ROT_NUM, the total number of rotations.
*/
{
  double *bw;
  double c;
  double g;
  double gapq;
  double h;
  int i;
  int j;
  int k;
  int l;
  int m;
  int p;
  int q;
  double s;
  double t;
  double tau;
  double term;
  double termp;
  double termq;
  double theta;
  double thresh;
  double w;
  double *zw;

  r8mat_identity ( n, v );

  r8mat_diag_get_vector ( n, a, d );

  bw = ( double * ) malloc ( n * sizeof ( double ) );
  zw = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    bw[i] = d[i];
    zw[i] = 0.0;
  }
  *it_num = 0;
  *rot_num = 0;

  while ( *it_num < it_max )
  {
    *it_num = *it_num + 1;
/*
  The convergence threshold is based on the size of the elements in
  the strict upper triangle of the matrix.
*/
    thresh = 0.0;
    for ( j = 0; j < n; j++ )
    {
      for ( i = 0; i < j; i++ )
      {
        thresh = thresh + a[i+j*n] * a[i+j*n];
      }
    }

    thresh = sqrt ( thresh ) / ( double ) ( 4 * n );

    if ( thresh == 0.0 )
    {
      break;
    }

    for ( p = 0; p < n; p++ )
    {
      for ( q = p + 1; q < n; q++ )
      {
        gapq = 10.0 * fabs ( a[p+q*n] );
        termp = gapq + fabs ( d[p] );
        termq = gapq + fabs ( d[q] );
/*
  Annihilate tiny offdiagonal elements.
*/
        if ( 4 < *it_num &&
             termp == fabs ( d[p] ) &&
             termq == fabs ( d[q] ) )
        {
          a[p+q*n] = 0.0;
        }
/*
  Otherwise, apply a rotation.
*/
        else if ( thresh <= fabs ( a[p+q*n] ) )
        {
          h = d[q] - d[p];
          term = fabs ( h ) + gapq;

          if ( term == fabs ( h ) )
          {
            t = a[p+q*n] / h;
          }
          else
          {
            theta = 0.5 * h / a[p+q*n];
            t = 1.0 / ( fabs ( theta ) + sqrt ( 1.0 + theta * theta ) );
            if ( theta < 0.0 )
            {
              t = - t;
            }
          }
          c = 1.0 / sqrt ( 1.0 + t * t );
          s = t * c;
          tau = s / ( 1.0 + c );
          h = t * a[p+q*n];
/*
  Accumulate corrections to diagonal elements.
*/
          zw[p] = zw[p] - h;                 
          zw[q] = zw[q] + h;
          d[p] = d[p] - h;
          d[q] = d[q] + h;

          a[p+q*n] = 0.0;
/*
  Rotate, using information from the upper triangle of A only.
*/
          for ( j = 0; j < p; j++ )
          {
            g = a[j+p*n];
            h = a[j+q*n];
            a[j+p*n] = g - s * ( h + g * tau );
            a[j+q*n] = h + s * ( g - h * tau );
          }

          for ( j = p + 1; j < q; j++ )
          {
            g = a[p+j*n];
            h = a[j+q*n];
            a[p+j*n] = g - s * ( h + g * tau );
            a[j+q*n] = h + s * ( g - h * tau );
          }

          for ( j = q + 1; j < n; j++ )
          {
            g = a[p+j*n];
            h = a[q+j*n];
            a[p+j*n] = g - s * ( h + g * tau );
            a[q+j*n] = h + s * ( g - h * tau );
          }
/*
  Accumulate information in the eigenvector matrix.
*/
          for ( j = 0; j < n; j++ )
          {
            g = v[j+p*n];
            h = v[j+q*n];
            v[j+p*n] = g - s * ( h + g * tau );
            v[j+q*n] = h + s * ( g - h * tau );
          }
          *rot_num = *rot_num + 1;
        }
      }
    }

    for ( i = 0; i < n; i++ )
    {
      bw[i] = bw[i] + zw[i];
      d[i] = bw[i];
      zw[i] = 0.0;
    }
  }
/*
  Restore upper triangle of input matrix.
*/
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < j; i++ )
    {
      a[i+j*n] = a[j+i*n];
    }
  }
/*
  Ascending sort the eigenvalues and eigenvectors.
*/
  for ( k = 0; k < n - 1; k++ )
  {
    m = k;
    for ( l = k + 1; l < n; l++ )
    {
      if ( d[l] < d[m] )
      {
        m = l;
      }
    }

    if ( m != k )
    {
      t    = d[m];
      d[m] = d[k];
      d[k] = t;
      for ( i = 0; i < n; i++ )
      {
        w        = v[i+m*n];
        v[i+m*n] = v[i+k*n];
        v[i+k*n] = w;
      }
    }
  }

  free ( bw );
  free ( zw );

  return;
}
/******************************************************************************/

void moment_method ( int n, double moment[], double x[], double w[] )

/******************************************************************************/
/*
  Purpose:

    MOMENT_METHOD computes a quadrature rule by the method of moments.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    18 September 2013

  Author:

    John Burkardt

  Reference:

    Gene Golub, John Welsch,
    Calculation of Gaussian Quadrature Rules,
    Mathematics of Computation,
    Volume 23, Number 106, April 1969, pages 221-230.

  Parameters:

    Input, int N, the order of the quadrature rule.

    Input, double MOMENT[2*N+1], moments 0 through 2*N.

    Output, double X[N], W[N], the points and weights of the quadrature rule.
*/
{
  double *alpha;
  double *beta;
  int debug;
  int flag;
  double *h;
  int i;
  int it_max;
  int it_num;
  int j;
  double *jacobi;
  double *r;
  int rot_num;
  double *v;

  debug = 0;

  if ( debug )
  {
    r8vec_print ( 2 * n + 1, moment, "  Moments:" );
  }

/*
  Define the N+1 by N+1 Hankel matrix H(I,J) = moment(I+J).
*/
  h = ( double * ) malloc ( ( n + 1 ) * ( n + 1 ) * sizeof ( double ) );

  for ( i = 0; i <= n; i++ )
  {
    for ( j = 0; j <= n; j++ )
    {
      h[i+j*(n+1)] = moment[i+j];
    }
  }

  if ( debug )
  {
    r8mat_print ( n + 1, n + 1, h, "  Hankel matrix:" );
  }
/*
  Compute R, the upper triangular Cholesky factor of H.
*/
  r = r8mat_cholesky_factor_upper ( n + 1, h, &flag );

  if ( flag != 0 )
  {
    printf ( "\n" );
    printf ( "MOMENT_METHOD - Fatal error!\n" );
    printf ( "  R8MAT_CHOLESKY_FACTOR_UPPER returned FLAG = %d\n", flag );
    exit ( 1 );
  }

  if ( debug )
  {
    r8mat_print ( n + 1, n + 1, r, "  Cholesky factor:" );
  }
/*
  Compute ALPHA and BETA from R, using Golub and Welsch's formula.
*/
  alpha = ( double * ) malloc ( n * sizeof ( double ) );

  alpha[0] = r[0+1*(n+1)] / r[0+0*(n+1)];
  for ( i = 1; i < n; i++ )
  {
    alpha[i] = r[i+(i+1)*(n+1)] / r[i+i*(n+1)] - r[i-1+i*(n+1)] / r[i-1+(i-1)*(n+1)];
  }
  beta = ( double * ) malloc ( ( n - 1 ) * sizeof ( double ) );

  for ( i = 0; i < n - 1; i++ )
  {
    beta[i] = r[i+1+(i+1)*(n+1)] / r[i+i*(n+1)];
  }
/*
  Compute the points and weights from the moments.
*/
  jacobi = ( double * ) malloc ( n * n * sizeof ( double ) );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      jacobi[i+j*n] = 0.0;
    }
  }

  for ( i = 0; i < n; i++ )
  {
    jacobi[i+i*n] = alpha[i];
  }

  for ( i = 0; i < n - 1; i++ )
  {
    jacobi[i+(i+1)*n] = beta[i];
    jacobi[i+1+i*n] = beta[i];
  }

  if ( debug )
  {
    r8mat_print ( n, n, jacobi, "  The Jacobi matrix:" );
  }
/*
  Get the eigendecomposition of the Jacobi matrix.
*/
  it_max = 100;
  v = ( double * ) malloc ( n * n * sizeof ( double ) );

  jacobi_eigenvalue ( n, jacobi, it_max, v, x, &it_num, &rot_num );

  if ( debug )
  {
    r8mat_print ( n, n, v, "  Eigenvector" );
  }

  for ( i = 0; i < n; i++ )
  {
    w[i] = moment[0] * pow ( v[0+i*n], 2 );
  }

  free ( alpha );
  free ( beta );
  free ( h );
  free ( jacobi );
  free ( r );
  free ( v );

  return;
}
/******************************************************************************/

double *moments_normal ( int m, double mu, double sigma )

/******************************************************************************/
/*
  Purpose:

    MOMENTS_NORMAL returns moments of the standard Normal distribution.

  Discussion:

    pdf(x) = exp ( -((x-mu)/sigma)^2/2 ) / sigma / sqrt ( pi * 2 )
    mu(k) = integral ( -oo < x < +oo ) x^k pdf(x) dx

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    17 September 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of moments desired.

    Input, double MU, SIGMA, the mean and standard deviation.

    Output, double W(0:M-1), the weighted integrals of X^0 
    through X^(M-1).
*/
{
  int j;
  int j_hi;
  int k;
  double t;
  double *w;

  w = ( double * ) malloc ( m * sizeof ( double ) );

  for ( k = 0; k < m; k++ )
  {
    t = 0.0;
    j_hi = k / 2;
    for ( j = 0; j <= j_hi; j++ )
    {
      t = t + r8_choose ( k, 2 * j ) * r8_factorial2 ( 2 * j - 1 ) 
        * pow ( sigma, 2 * j ) * pow ( mu, k - 2 * j );
    }
    w[k] = t;
  }

  return w;
}
/******************************************************************************/

double *moments_truncated_normal_ab ( int m, double mu, double sigma,
  double a, double b )

/******************************************************************************/
/*
  Purpose:

    MOMENTS_TRUNCATED_NORMAL_AB: moments of the truncated Normal distribution.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    19 September 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of moments desired.

    Input, double MU, SIGMA, the mean and standard deviation.

    Input, double A, B, the lower and upper truncation limits.

    Output, double W(0:M-1), the weighted integrals of X^0 
    through X^(M-1).
*/
{
  int order;
  double *w;

  w = ( double * ) malloc ( m * sizeof ( double ) );

  for ( order = 0; order < m; order++ )
  {
    w[order] = truncated_normal_ab_moment ( order, mu, sigma, a, b );
  }

  return w;
}
/******************************************************************************/

double *moments_truncated_normal_a ( int m, double mu, double sigma,
  double a )

/******************************************************************************/
/*
  Purpose:

    MOMENTS_TRUNCATED_NORMAL_A: moments of the lower truncated Normal.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    19 September 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of moments desired.

    Input, double MU, SIGMA, the mean and standard deviation.

    Input, double A, the lower truncation limit.

    Output, double W(0:M-1), the weighted integrals of X^0 
    through X^(M-1).
*/
{
  int order;
  double *w;

  w = ( double * ) malloc ( m * sizeof ( double ) );

  for ( order = 0; order < m; order++ )
  {
    w[order] = truncated_normal_a_moment ( order, mu, sigma, a );
  }

  return w;
}
/******************************************************************************/

double *moments_truncated_normal_b ( int m, double mu, double sigma,
  double b )

/******************************************************************************/
/*
  Purpose:

    MOMENTS_TRUNCATED_NORMAL_B: moments of the upper truncated Normal.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    19 September 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of moments desired.

    Input, double MU, SIGMA, the mean and standard deviation.

    Input, double B, the upper truncation limit.

    Output, double W(0:M-1), the weighted integrals of X^0 
    through X^(M-1).
*/
{
  int order;
  double *w;

  w = ( double * ) malloc ( m * sizeof ( double ) );

  for ( order = 0; order < m; order++ )
  {
    w[order] = truncated_normal_b_moment ( order, mu, sigma, b );
  }

  return w;
}
/******************************************************************************/

double normal_01_cdf ( double x )

/******************************************************************************/
/*
  Purpose:

    NORMAL_01_CDF evaluates the Normal 01 CDF.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    10 February 1999

  Author:

    John Burkardt

  Reference:

    A G Adams,
    Areas Under the Normal Curve,
    Algorithm 39,
    Computer j.,
    Volume 12, pages 197-198, 1969.

  Parameters:

    Input, double X, the argument of the CDF.

    Output, double CDF, the value of the CDF.
*/
{
  double a1 = 0.398942280444;
  double a2 = 0.399903438504;
  double a3 = 5.75885480458;
  double a4 = 29.8213557808;
  double a5 = 2.62433121679;
  double a6 = 48.6959930692;
  double a7 = 5.92885724438;
  double b0 = 0.398942280385;
  double b1 = 3.8052E-08;
  double b2 = 1.00000615302;
  double b3 = 3.98064794E-04;
  double b4 = 1.98615381364;
  double b5 = 0.151679116635;
  double b6 = 5.29330324926;
  double b7 = 4.8385912808;
  double b8 = 15.1508972451;
  double b9 = 0.742380924027;
  double b10 = 30.789933034;
  double b11 = 3.99019417011;
  double cdf;
  double q;
  double y;
/*
  |X| <= 1.28.
*/
  if ( fabs ( x ) <= 1.28 )
  {
    y = 0.5 * x * x;

    q = 0.5 - fabs ( x ) * ( a1 - a2 * y / ( y + a3 - a4 / ( y + a5
      + a6 / ( y + a7 ) ) ) );
/*
  1.28 < |X| <= 12.7
*/
  }
  else if ( fabs ( x ) <= 12.7 )
  {
    y = 0.5 * x * x;

    q = exp ( - y ) * b0 / ( fabs ( x ) - b1
      + b2 / ( fabs ( x ) + b3
      + b4 / ( fabs ( x ) - b5
      + b6 / ( fabs ( x ) + b7
      - b8 / ( fabs ( x ) + b9
      + b10 / ( fabs ( x ) + b11 ) ) ) ) ) );
/*
  12.7 < |X|
*/
  }
  else
  {
    q = 0.0;
  }
/*
  Take account of negative X.
*/
  if ( x < 0.0 )
  {
    cdf = q;
  }
  else
  {
    cdf = 1.0 - q;
  }

  return cdf;
}
/******************************************************************************/

double normal_01_pdf ( double x )

/******************************************************************************/
/*
  Purpose:

    NORMAL_01_PDF evaluates the Normal 01 PDF.

  Discussion:

    The Normal 01 PDF is also called the "Standard Normal" PDF, or
    the Normal PDF with 0 mean and variance 1.

    PDF(X) = exp ( - 0.5 * X^2 ) / sqrt ( 2 * PI )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    18 September 2004

  Author:

    John Burkardt

  Parameters:

    Input, double X, the argument of the PDF.

    Output, double PDF, the value of the PDF.
*/
{
  double pdf;
  const double pi = 3.14159265358979323;

  pdf = exp ( -0.5 * x * x ) / sqrt ( 2.0 * pi );

  return pdf;
}
/******************************************************************************/

double r8_choose ( int n, int k )

/******************************************************************************/
/*
  Purpose:

    R8_CHOOSE computes the binomial coefficient C(N,K) as an R8.

  Discussion:

    The value is calculated in such a way as to avoid overflow and
    roundoff.  The calculation is done in R8 arithmetic.

    The formula used is:

      C(N,K) = N! / ( K! * (N-K)! )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 July 2008

  Author:

    John Burkardt

  Reference:

    ML Wolfson, HV Wright,
    Algorithm 160:
    Combinatorial of M Things Taken N at a Time,
    Communications of the ACM,
    Volume 6, Number 4, April 1963, page 161.

  Parameters:

    Input, int N, K, the values of N and K.

    Output, double R8_CHOOSE, the number of combinations of N
    things taken K at a time.
*/
{
  int i;
  int mn;
  int mx;
  double value;

  if ( k < n - k )
  {
    mn = k;
    mx = n - k;
  }
  else
  {
    mn = n - k;
    mx = k;
  }

  if ( mn < 0 )
  {
    value = 0.0;
  }
  else if ( mn == 0 )
  {
    value = 1.0;
  }
  else
  {
    value = ( double ) ( mx + 1 );

    for ( i = 2; i <= mn; i++ )
    {
      value = ( value * ( double ) ( mx + i ) ) / ( double ) i;
    }
  }

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

double r8_factorial ( int n )

/******************************************************************************/
/*
  Purpose:

    R8_FACTORIAL computes the factorial of N.

  Discussion:

    factorial ( N ) = product ( 1 <= I <= N ) I

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2008

  Author:

    John Burkardt

  Parameters:

    Input, int N, the argument of the factorial function.
    If N is less than 1, the function value is returned as 1.

    Output, double R8_FACTORIAL, the factorial of N.
*/
{
  int i;
  double value;

  value = 1.0;

  for ( i = 1; i <= n; i++ )
  {
    value = value * ( double ) ( i );
  }

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

double r8_factorial2 ( int n )

/******************************************************************************/
/*
  Purpose:

    R8_FACTORIAL2 computes the double factorial function.

  Discussion:

    FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 )  (N even)
                    = Product ( N * (N-2) * (N-4) * ... * 1 )  (N odd)

     N Value
    -- -----
     0     1
     1     1
     2     2
     3     3
     4     8
     5    15
     6    48
     7   105
     8   384
     9   945
    10  3840

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 July 2008

  Author:

    John Burkardt

  Parameters:

    Input, int N, the argument of the double factorial
    function.  If N is less than 1, R8_FACTORIAL2 is returned as 1.0.

    Output, double R8_FACTORIAL2, the value of Factorial2(N).
*/
{
  int n_copy;
  double value;

  value = 1.0;

  if ( n < 1 )
  {
    return value;
  }

  n_copy = n;

  while ( 1 < n_copy )
  {
    value = value * ( double ) n_copy;
    n_copy = n_copy - 2;
  }

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

double r8_mop ( int i )

/******************************************************************************/
/*
  Purpose:

    R8_MOP returns the I-th power of -1 as an R8 value.

  Discussion:

    An R8 is an double value.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 July 2008

  Author:

    John Burkardt

  Parameters:

    Input, int I, the power of -1.

    Output, double R8_MOP, the I-th power of -1.
*/
{
  double value;

  if ( ( i % 2 ) == 0 )
  {
    value = + 1.0;
  }
  else
  {
    value = - 1.0;
  }

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

double *r8mat_cholesky_factor_upper ( int n, double a[], int *flag )

/******************************************************************************/
/*
  Purpose:

    R8MAT_CHOLESKY_FACTOR_UPPER: upper Cholesky factor of a symmetric R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

    The matrix must be symmetric and positive semidefinite.

    For a positive semidefinite symmetric matrix A, the Cholesky factorization
    is an upper triangular matrix R such that:

      A = R' * R

    Note that the usual Cholesky factor is a LOWER triangular matrix L
    such that

      A = L * L'

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    03 August 2013

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of rows and columns of the matrix A.

    Input, double A[N*N], the N by N matrix.

    Output, int *FLAG, an error flag.
    0, no error was detected.
    1, the matrix was not positive definite.  A NULL factor was returned.

    Output, double R8MAT_CHOLESKY_FACTOR_UPPER[N*N], the N by N upper triangular
    "Choresky" factor.
*/
{
  double *c;
  int i;
  int j;
  int k;
  double sum2;

  *flag = 0;

  c = r8mat_copy_new ( n, n, a );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < j; i++ )
    {
      c[j+i*n] = 0.0;
    }
    for ( i = j; i < n; i++ )
    {
      sum2 = c[i+j*n];
      for ( k = 0; k < j; k++ )
      {
        sum2 = sum2 - c[k+j*n] * c[k+i*n];
      }
      if ( i == j )
      {
        if ( sum2 <= 0.0 )
        {
          *flag = 1;
          return NULL;
        }
        c[j+i*n] = sqrt ( sum2 );
      }
      else
      {
        if ( c[j+j*n] != 0.0 )
        {
          c[j+i*n] = sum2 / c[j+j*n];
        }
        else
        {
          c[j+i*n] = 0.0;
        }
      }
    }
  }

  return c;
}
/******************************************************************************/

double *r8mat_copy_new ( int m, int n, double a1[] )

/******************************************************************************/
/*
  Purpose:

    R8MAT_COPY_NEW copies one R8MAT to a "new" R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 July 2008

  Author:

    John Burkardt

  Parameters:

    Input, int M, N, the number of rows and columns.

    Input, double A1[M*N], the matrix to be copied.

    Output, double R8MAT_COPY_NEW[M*N], the copy of A1.
*/
{
  double *a2;
  int i;
  int j;

  a2 = ( double * ) malloc ( m * n * sizeof ( double ) );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      a2[i+j*m] = a1[i+j*m];
    }
  }

  return a2;
}
/******************************************************************************/

void r8mat_diag_get_vector ( int n, double a[], double v[] )

/******************************************************************************/
/*
  Purpose:

    R8MAT_DIAG_GET_VECTOR gets the value of the diagonal of an R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    15 July 2013

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of rows and columns of the matrix.

    Input, double A[N*N], the N by N matrix.

    Output, double V[N], the diagonal entries
    of the matrix.
*/
{
  int i;

  for ( i = 0; i < n; i++ )
  {
    v[i] = a[i+i*n];
  }

  return;
}
/******************************************************************************/

void r8mat_identity  ( int n, double a[] )

/******************************************************************************/
/*
  Purpose:

    R8MAT_IDENTITY sets an R8MAT to the identity matrix.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    06 September 2005

  Author:

    John Burkardt

  Parameters:

    Input, int N, the order of A.

    Output, double A[N*N], the N by N identity matrix.
*/
{
  int i;
  int j;
  int k;

  k = 0;
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < n; i++ )
    {
      if ( i == j )
      {
        a[k] = 1.0;
      }
      else
      {
        a[k] = 0.0;
      }
      k = k + 1;
    }
  }

  return;
}
/******************************************************************************/

void r8mat_print ( int m, int n, double a[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8MAT_PRINT prints an R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

    Entry A(I,J) is stored as A[I+J*M]

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    28 May 2008

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows in A.

    Input, int N, the number of columns in A.

    Input, double A[M*N], the M by N matrix.

    Input, char *TITLE, a title.
*/
{
  r8mat_print_some ( m, n, a, 1, 1, m, n, title );

  return;
}
/******************************************************************************/

void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi,
  int jhi, char *title )

/******************************************************************************/
/*
  Purpose:

    R8MAT_PRINT_SOME prints some of an R8MAT.

  Discussion:

    An R8MAT is a doubly dimensioned array of R8 values, stored as a vector
    in column-major order.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    26 June 2013

  Author:

    John Burkardt

  Parameters:

    Input, int M, the number of rows of the matrix.
    M must be positive.

    Input, int N, the number of columns of the matrix.
    N must be positive.

    Input, double A[M*N], the matrix.

    Input, int ILO, JLO, IHI, JHI, designate the first row and
    column, and the last row and column to be printed.

    Input, char *TITLE, a title.
*/
{
# define INCX 5

  int i;
  int i2hi;
  int i2lo;
  int j;
  int j2hi;
  int j2lo;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );

  if ( m <= 0 || n <= 0 )
  {
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  (None)\n" );
    return;
  }
/*
  Print the columns of the matrix, in strips of 5.
*/
  for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX )
  {
    j2hi = j2lo + INCX - 1;
    if ( n < j2hi )
    {
      j2hi = n;
    }
    if ( jhi < j2hi )
    {
      j2hi = jhi;
    }

    fprintf ( stdout, "\n" );
/*
  For each column J in the current range...

  Write the header.
*/
    fprintf ( stdout, "  Col:  ");
    for ( j = j2lo; j <= j2hi; j++ )
    {
      fprintf ( stdout, "  %7d     ", j - 1 );
    }
    fprintf ( stdout, "\n" );
    fprintf ( stdout, "  Row\n" );
    fprintf ( stdout, "\n" );
/*
  Determine the range of the rows in this strip.
*/
    if ( 1 < ilo )
    {
      i2lo = ilo;
    }
    else
    {
      i2lo = 1;
    }
    if ( m < ihi )
    {
      i2hi = m;
    }
    else
    {
      i2hi = ihi;
    }

    for ( i = i2lo; i <= i2hi; i++ )
    {
/*
  Print out (up to) 5 entries in row I, that lie in the current strip.
*/
      fprintf ( stdout, "%5d:", i - 1 );
      for ( j = j2lo; j <= j2hi; j++ )
      {
        fprintf ( stdout, "  %14g", a[i-1+(j-1)*m] );
      }
      fprintf ( stdout, "\n" );
    }
  }

  return;
# undef INCX
}
/******************************************************************************/

void r8mat_write ( char *output_filename, int m, int n, double table[] )

/******************************************************************************/
/*
  Purpose:

    R8MAT_WRITE writes an R8MAT file.

  Discussion:

    An R8MAT is an array of R8's.

  Licensing:

    This code is distributed under the MIT license. 

  Modified:

    01 June 2009

  Author:

    John Burkardt

  Parameters:

    Input, char *OUTPUT_FILENAME, the output filename.

    Input, int M, the spatial dimension.

    Input, int N, the number of points.

    Input, double TABLE[M*N], the data.
*/
{
  int i;
  int j;
  FILE *output;
/*
  Open the file.
*/
  output = fopen ( output_filename, "wt" );

  if ( !output )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "R8MAT_WRITE - Fatal error!\n" );
    fprintf ( stderr, "  Could not open the output file.\n" );
    exit ( 1 );
  }
/*
  Write the data.
*/
  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      fprintf ( output, "  %24.16g", table[i+j*m] );
    }
    fprintf ( output, "\n" );
  }
/*
  Close the file.
*/
  fclose ( output );

  return;
}
/******************************************************************************/

void r8vec_print ( int n, double a[], char *title )

/******************************************************************************/
/*
  Purpose:

    R8VEC_PRINT prints an R8VEC.

  Discussion:

    An R8VEC is a vector of R8's.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    08 April 2009

  Author:

    John Burkardt

  Parameters:

    Input, int N, the number of components of the vector.

    Input, double A[N], the vector to be printed.

    Input, char *TITLE, a title.
*/
{
  int i;

  fprintf ( stdout, "\n" );
  fprintf ( stdout, "%s\n", title );
  fprintf ( stdout, "\n" );
  for ( i = 0; i < n; i++ )
  {
    fprintf ( stdout, "  %8d: %14g\n", i, a[i] );
  }

  return;
}
/******************************************************************************/

void rule_write ( int order, char *filename, double x[], double w[],
  double r[] )

/******************************************************************************/
/*
  Purpose:

    RULE_WRITE writes a quadrature rule to three files.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    01 October 2012

  Author:

    John Burkardt

  Parameters:

    Input, int ORDER, the order of the rule.

    Input, double A, the left endpoint.

    Input, double B, the right endpoint.

    Input, char *FILENAME, specifies the output filenames.
    "filename_w.txt", "filename_x.txt", "filename_r.txt"
    defining weights, abscissas, and region.
*/
{
  char filename_r[80];
  char filename_w[80];
  char filename_x[80];

  strcpy ( filename_r, filename );
  strcat ( filename_r, "_r.txt" );
  strcpy ( filename_w, filename );
  strcat ( filename_w, "_w.txt" );
  strcpy ( filename_x, filename );
  strcat ( filename_x, "_x.txt" );

  printf ( "\n" );
  printf ( "  Creating quadrature files.\n" );
  printf ( "\n" );
  printf ( "  Root file name is     \"%s\".\n", filename );
  printf ( "\n" );
  printf ( "  Weight file will be   \"%s\".\n", filename_w );
  printf ( "  Abscissa file will be \"%s\".\n", filename_x );
  printf ( "  Region file will be   \"%s\".\n", filename_r );

  r8mat_write ( filename_w, 1, order, w );
  r8mat_write ( filename_x, 1, order, x );
  r8mat_write ( filename_r, 1, 2,     r );

  return;
}
/******************************************************************************/

void timestamp ( void )

/******************************************************************************/
/*
  Purpose:

    TIMESTAMP prints the current YMDHMS date as a time stamp.

  Example:

    31 May 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 );

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}
/******************************************************************************/

double truncated_normal_ab_moment ( int order, double mu, double s, double a,
  double b )

/******************************************************************************/
/*
  Purpose:

    TRUNCATED_NORMAL_AB_MOMENT: moments of the truncated Normal PDF.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 September 2013

  Author:

    John Burkardt

  Reference:

    Phoebus Dhrymes,
    Moments of Truncated Normal Distributions,
    May 2005.

  Parameters:

    Input, int ORDER, the order of the moment.
    0 <= ORDER.

    Input, double MU, S, the mean and standard deviation of the
    parent Normal distribution.
    0.0 < S.

    Input, double A, B, the lower and upper truncation limits.
    A < B.

    Output, double TRUNCATED_NORMAL_AB_MOMENT, the moment of the PDF.
*/
{
  double a_h;
  double a_cdf;
  double a_pdf;
  double b_h;
  double b_cdf;
  double b_pdf;
  double ir;
  double irm1;
  double irm2;
  double moment;
  int r;

  if ( order < 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" );
    fprintf ( stderr, "  ORDER < 0.\n" );
    exit ( 1 );
  }

  if ( s <= 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" );
    fprintf ( stderr, "  S <= 0.\n" );
    exit ( 1 );
  }

  if ( b <= a )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" );
    fprintf ( stderr, "  B <= A.\n" );
    exit ( 1 );
  }

  a_h = ( a - mu ) / s;
  a_pdf = normal_01_pdf ( a_h );
  a_cdf = normal_01_cdf ( a_h );

  if ( a_cdf == 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" );
    fprintf ( stderr, "  PDF/CDF ratio fails, A_CDF too small.\n" );
    fprintf ( stderr, "  A_PDF = %g\n", a_pdf );
    fprintf ( stderr, "  A_CDF = %g\n", a_cdf );
    exit ( 1 );
  }

  b_h = ( b - mu ) / s;
  b_pdf = normal_01_pdf ( b_h );
  b_cdf = normal_01_cdf ( b_h );

  if ( b_cdf == 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" );
    fprintf ( stderr, "  PDF/CDF ratio fails, B_CDF too small.\n" );
    fprintf ( stderr, "  B_PDF = %g\n", b_pdf );
    fprintf ( stderr, "  B_CDF = %g\n", b_cdf );
    exit ( 1 );
  }

  moment = 0.0;
  irm2 = 0.0;
  irm1 = 0.0;

  for ( r = 0; r <= order; r++ )
  {
    if ( r == 0 )
    {
      ir = 1.0;
    }
    else if ( r == 1 )
    {
      ir = - ( b_pdf - a_pdf ) / ( b_cdf - a_cdf );
    }
    else
    {
      ir = ( double ) ( r - 1 ) * irm2 
        - ( pow ( b_h, r - 1 ) * b_pdf - pow ( a_h, r - 1 ) * a_pdf )
        / ( b_cdf - a_cdf );
    }

    moment = moment + r8_choose ( order, r ) * pow ( mu, order - r ) 
      * pow ( s, r ) * ir;

    irm2 = irm1;
    irm1 = ir;
  }

  return moment;
}
/******************************************************************************/

double truncated_normal_a_moment ( int order, double mu, double s, double a )

/******************************************************************************/
/*
  Purpose:

    TRUNCATED_NORMAL_A_MOMENT: moments of the lower truncated Normal PDF.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 September 2013

  Author:

    John Burkardt

  Reference:

    Phoebus Dhrymes,
    Moments of Truncated Normal Distributions,
    May 2005.

  Parameters:

    Input, int ORDER, the order of the moment.
    0 <= ORDER.

    Input, double MU, S, the mean and standard deviation of the
    parent Normal distribution.

    Input, double A, the lower truncation limit.

    Output, double TRUNCATED_NORMAL_A_MOMENT, the moment of the PDF.
*/
{
  double moment;

  moment = r8_mop ( order ) 
    * truncated_normal_b_moment ( order, - mu, s, - a );

  return moment;
}
/******************************************************************************/

double truncated_normal_b_moment ( int order, double mu, double s, double b )

/******************************************************************************/
/*
  Purpose:

    TRUNCATED_NORMAL_B_MOMENT: moments of the upper truncated Normal PDF.

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    11 September 2013

  Author:

    John Burkardt

  Reference:

    Phoebus Dhrymes,
    Moments of Truncated Normal Distributions,
    May 2005.

  Parameters:

    Input, int ORDER, the order of the moment.
    0 <= ORDER.

    Input, double MU, S, the mean and standard deviation of the
    parent Normal distribution.

    Input, double B, the upper truncation limit.

    Output, double TRUNCATED_NORMAL_B_MOMENT, the moment of the PDF.
*/
{
  double f;
  double h;
  double h_cdf;
  double h_pdf;
  double ir;
  double irm1;
  double irm2;
  double moment;
  int r;

  if ( order < 0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n" );
    fprintf ( stderr, "  ORDER < 0.\n" );
    exit ( 1 );
  }

  h = ( b - mu ) / s;
  h_pdf = normal_01_pdf ( h );
  h_cdf = normal_01_cdf ( h );

  if ( h_cdf == 0.0 )
  {
    fprintf ( stderr, "\n" );
    fprintf ( stderr, "TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n" );
    fprintf ( stderr, "  CDF((B-MU)/S) = 0.\n" );
    exit ( 1 );
  }

  f = h_pdf / h_cdf;

  moment = 0.0;
  irm2 = 0.0;
  irm1 = 0.0;

  for ( r = 0; r <= order; r++ )
  {
    if ( r == 0 )
    {
      ir = 1.0;
    }
    else if ( r == 1 )
    {
      ir = - f;
    }
    else
    {
      ir = - pow ( h, r - 1 ) * f + ( double ) ( r - 1 ) * irm2;
    }

    moment = moment + r8_choose ( order, r ) * pow ( mu, order - r ) 
      * pow ( s, r ) * ir;

    irm2 = irm1;
    irm1 = ir;
  }

  return moment;
}