# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <time.h>
# include <string.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_huge ( );
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 GNU LGPL 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 = - r8_huge ( );
  }
/*
  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 = r8_huge ( );
  }

  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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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_huge ( void )

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

    R8_HUGE returns a "huge" R8.

  Discussion:

    The value returned by this function is NOT required to be the
    maximum representable R8.  This value varies from machine to machine,
    from compiler to compiler, and may cause problems when being printed.
    We simply want a "very large" but non-infinite number.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    06 October 2007

  Author:

    John Burkardt

  Parameters:

    Output, double R8_HUGE, a "huge" R8 value.
*/
{
  double value;

  value = 1.0E+30;

  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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL 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;
}