double r8mat_amax ( int m, int n, double a[] )

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

    R8MAT_AMAX returns the maximum absolute value entry 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:

    07 September 2012

  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.

    Output, double R8MAT_AMAX, the maximum absolute value entry of A.
*/
{
  int i;
  int j;
  double value;

  value = r8_abs ( a[0+0*m] );

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      if ( value < r8_abs ( a[i+j*m] ) )
      {
        value = r8_abs ( a[i+j*m] );
      }
    }
  }
  return value;
}

double zabs1 ( _Complex double z )

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

    ZABS1 returns the L1 norm of a number.

  Discussion:

    This routine uses double precision complex arithmetic.

    The L1 norm of a complex number is the sum of the absolute values
    of the real and imaginary components.

    ZABS1 ( Z ) = abs ( real ( Z ) ) + abs ( imaginary ( Z ) )

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    31 March 2007

  Author:

    C version by John Burkardt

  Reference:

    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
    LINPACK User's Guide,
    SIAM, 1979.

    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
    Basic Linear Algebra Subprograms for FORTRAN usage,
    ACM Transactions on Mathematical Software,
    Volume 5, Number 3, pages 308-323, 1979.

  Parameters:

    Input, _Complex double Z, the number whose norm is desired.

    Output, double ZABS1, the L1 norm of Z.
*/
{
  double value;

  value = r8_abs ( creal ( z ) ) + r8_abs ( cimagf ( z ) );

  return value;
}

int r8_nint ( double x )
{
  int value;

  if ( x < 0.0 )
  {
    value = - ( int ) ( r8_abs ( x ) + 0.5 );
  }
  else
  {
    value =   ( int ) ( r8_abs ( x ) + 0.5 );
  }

  return value;
}

void test05 ( )

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

    TEST05 tests P00_GAUSS_HERMITE against the polynomials.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    11 September 2012

  Author:

    John Burkardt
*/
{
  double error;
  double estimate;
  double exact;
  int m;
  int order;
  int order_log;
  int problem;

  printf ( "\n" );
  printf ( "TEST05\n" );
  printf ( "  P00_GAUSS_HERMITE applies a Gauss-Hermite rule to\n" );
  printf ( "  estimate the integral x^m exp(-x*x) over (-oo,+oo).\n" );

  problem = 6;

  printf ( "\n" );
  printf ( "         M     Order      Estimate        Exact           Error\n" );

  for ( m = 0; m <= 6; m++ )
  {
    p06_param ( 'S', 'M', &m );

    exact = p00_exact ( problem );

    printf ( "\n" );

    for ( order = 1; order <= 3 + ( m / 2 ); order++ )
    {
      estimate = p00_gauss_hermite ( problem, order );

      error = r8_abs ( exact - estimate );

      printf ( "  %8d  %8d  %14g  %14g  %14g\n", 
        m, order, estimate, exact, error );
    }
  }

  return;
}

void test02 ( void )

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

    TEST02 demonstrates the use of BIVNOR.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    13 April 2012

  Author:

    John Burkardt
*/
{
  double fxy1;
  double fxy2;
  int n_data;
  double r;
  double x;
  double y;

  printf ( "\n" );
  printf ( "TEST02:\n" );
  printf ( "  BIVNOR computes the bivariate normal probability.\n" );
  printf ( "  Compare to tabulated values.\n" );
  printf ( "\n" );
  printf ( "          X               Y               " );
  printf ( "R           P                         P                       DIFF\n" );
  printf ( "                                          " );
  printf ( "           (Tabulated)               (BIVNOR)\n" );
  printf ( "\n" );

  n_data = 0;

  for ( ; ; )
  {
    bivariate_normal_cdf_values ( &n_data, &x, &y, &r, &fxy1 );

    if ( n_data == 0 )
    {
      break;
    }

    fxy2 = bivnor ( - x, - y, r );

    printf ( "  %12.4f  %12.4f  %12.4f  %24.16f  %24.16f  %10.4g\n",
      x, y, r, fxy1, fxy2, r8_abs ( fxy1 - fxy2 ) );
  }

  return;
}

void test01 ( void )

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

    TEST01 demonstrates the use of T.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    21 November 2010

  Author:

    John Burkardt
*/
{
  double a;
  double h;
  int n_data;
  double t1;
  double t2;

  printf ( "\n" );
  printf ( "TEST01:\n" );
  printf ( "  T computes Owen's T function.\n" );
  printf ( "  Compare to tabulated values.\n" );
  printf ( "\n" );
  printf ( "             H             A      " );
  printf ( "    T                         T  \n" );
  printf ( "                                  " );
  printf ( "    (Tabulated)               (TFN)               DIFF\n" );
  printf ( "\n" );

  n_data = 0;

  for ( ; ; )
  {
    owen_values ( &n_data, &h, &a, &t1 );

    if ( n_data == 0 )
    {
      break;
    }

    t2 = t ( h, a );

    printf ( "  %12.4f  %12.4f  %24.16f  %24.16f  %10.4g\n",
      h, a, t1, t2, r8_abs ( t1 - t2 ) );
  }

  return;
}

void test04 ( void )

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

    TEST04 demonstrates the use of ZNORM2.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    23 November 2010

  Author:

    John Burkardt
*/
{
  double fx1;
  double fx2;
  int n_data;
  double x;

  printf ( "\n" );
  printf ( "TEST04:\n" );
  printf ( "  ZNORM2 computes the complementary normal CDF.\n" );
  printf ( "  Compare to tabulated values.\n" );
  printf ( "\n" );
  printf ( "          X           P                         P                       DIFF\n" );
  printf ( "                     (Tabulated)               (ZNORM2)\n" );
  printf ( "\n" );

  n_data = 0;

  for ( ; ; )
  {
    normal_01_cdf_values ( &n_data, &x, &fx1 );

    if ( n_data == 0 )
    {
      break;
    }

    fx1 = 1.0 - fx1;

    fx2 = znorm2 ( x );

    printf ( "  %12.4f  %24.16f  %24.16f  %10.4g\n",
      x, fx1, fx2, r8_abs ( fx1 - fx2 ) );
  }

  return;
}

void test02 ( void )

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

    TEST02 demonstrates the use of MDBETA.

  Modified:

    30 January 2008

  Author:

    John Burkardt
*/
{
  double fx;
  double fx2;
  int ier;
  int n_data;
  double p;
  double q;
  double x;

  printf ( "\n" );
  printf ( "TEST02:\n" );
  printf ( "  MDBETA estimates the value of th modified Beta function.\n" );
  printf ( "  Compare with tabulated values.\n" );
  printf ( "\n" );
  printf ( "         X         P         Q         Beta                       Beta                  DIFF\n" );
  printf ( "                                       (Tabulated)                (MDBETA)\n" );
  printf ( "\n" );

  n_data = 0;

  for ( ; ; )
  {
    beta_cdf_values ( &n_data, &p, &q, &x, &fx );

    if ( n_data == 0 )
    {
      break;
    }

    fx2 = mdbeta ( x, p, q, &ier );

    printf ( "  %8.4f  %8.4f  %8.4f  %24.16e  %24.16e %10.4e\n",
      x, p, q, fx, fx2, r8_abs ( fx - fx2 ) );
  }

  return;
}

void test01 ( void )

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

    TEST01 demonstrates the use of ALOGAM.

  Modified:

    30 January 2008

  Author:

    John Burkardt
*/
{
  double fx;
  double fx2;
  int ifault;
  int n_data;
  double x;

  printf ( "\n" );
  printf ( "TEST01:\n" );
  printf ( "  ALOGAM computes the logarithm of the \n" );
  printf ( "  Gamma function.  We compare the result\n" );
  printf ( "  to tabulated values.\n" );
  printf ( "\n" );
  printf ( "          X                     FX                        FX2\n" );
  printf ( "                                (Tabulated)               (ALOGAM)                DIFF\n" );
  printf ( "\n" );

  n_data = 0;

  for ( ; ; )
  {
    gamma_log_values ( &n_data, &x, &fx );

    if ( n_data == 0 )
    {
      break;
    }

    fx2 = alogam ( x, &ifault );

    printf ( "  %24.16e  %24.16e  %24.16e  %10.4e\n",
      x, fx, fx2, r8_abs ( fx - fx2 ) );
  }

  return;
}

void dqrank(double a[], int lda, int m, int n, double tol, int* kr,
            int jpvt[], double qraux[])

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

    DQRANK computes the QR factorization of a rectangular matrix.

  Discussion:

    This routine is used in conjunction with DQRLSS to solve
    overdetermined, underdetermined and singular linear systems
    in a least squares sense.

    DQRANK uses the LINPACK subroutine DQRDC to compute the QR
    factorization, with column pivoting, of an M by N matrix A.
    The numerical rank is determined using the tolerance TOL.

    Note that on output, ABS ( A(1,1) ) / ABS ( A(KR,KR) ) is an estimate
    of the condition number of the matrix of independent columns,
    and of R.  This estimate will be <= 1/TOL.

  Licensing:

    This code is distributed under the GNU LGPL license.

  Modified:

    21 April 2012

  Author:

    C version by John Burkardt.

  Reference:

    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
    LINPACK User's Guide,
    SIAM, 1979,
    ISBN13: 978-0-898711-72-1,
    LC: QA214.L56.

  Parameters:

    Input/output, double A[LDA*N].  On input, the matrix whose
    decomposition is to be computed.  On output, the information from DQRDC.
    The triangular matrix R of the QR factorization is contained in the
    upper triangle and information needed to recover the orthogonal
    matrix Q is stored below the diagonal in A and in the vector QRAUX.

    Input, int LDA, the leading dimension of A, which must
    be at least M.

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

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

    Input, double TOL, a relative tolerance used to determine the
    numerical rank.  The problem should be scaled so that all the elements
    of A have roughly the same absolute accuracy, EPS.  Then a reasonable
    value for TOL is roughly EPS divided by the magnitude of the largest
    element.

    Output, int *KR, the numerical rank.

    Output, int JPVT[N], the pivot information from DQRDC.
    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
    independent to within the tolerance TOL and the remaining columns
    are linearly dependent.

    Output, double QRAUX[N], will contain extra information defining
    the QR factorization.
*/
{
  double work[n];

  for (int i = 0; i < n; i++)
    jpvt[i] = 0;

  int job = 1;

  dqrdc(a, lda, m, n, qraux, jpvt, work, job);

  *kr = 0;
  int k = i4_min(m, n);
  for (int j = 0; j < k; j++) {
    if (r8_abs(a[j + j * lda]) <= tol * r8_abs(a[0 + 0 * lda]))
      return;
    *kr = j + 1;
  }
}

void r8_test ( int logn_max )

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

    R8_TEST estimates the value of PI using double.

  Discussion:

    PI is estimated using N terms.  N is increased from 10^2 to 10^LOGN_MAX.
    The calculation is repeated using both sequential and Open MP enabled code.
    Wall clock time is measured by calling SYSTEM_CLOCK.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    13 November 2007

  Author:

    John Burkardt
*/
{
  double error;
  double estimate;
  int logn;
  char mode[4];
  int n;
  double r8_pi = 3.141592653589793;
  double wtime;

  printf ( "\n" );
  printf ( "R8_TEST:\n" );
  printf ( "  Estimate the value of PI,\n" );
  printf ( "  using double arithmetic.\n" );
  printf ( "\n" );
  printf ( "  N = number of terms computed and added;\n" );
  printf ( "\n" );
  printf ( "  MODE = SEQ for sequential code;\n" );
  printf ( "  MODE = OMP for Open MP enabled code;\n" );
  printf ( "  (performance depends on whether Open MP is used,\n" );
  printf ( "  and how many processes are available)\n" );
  printf ( "\n" );
  printf ( "  ESTIMATE = the computed estimate of PI;\n" );
  printf ( "\n" );
  printf ( "  ERROR = ( the computed estimate - PI );\n" );
  printf ( "\n" );
  printf ( "  TIME = elapsed wall clock time;\n" );
  printf ( "\n" );
  printf ( "  Note that you can''t increase N forever, because:\n" );
  printf ( "  A) ROUNDOFF starts to be a problem, and\n" );
  printf ( "  B) maximum integer size is a problem.\n" );
  printf ( "\n" );
  printf ( "             N Mode    Estimate        Error           Time\n" );
  printf ( "\n" );

  n = 1;

  for ( logn = 1; logn <= logn_max; logn++ )
  {
/*
  Note that when I set N = 10**LOGN directly, rather than using
  recursion, I got inaccurate values of N when LOGN was "large",
  that is, for LOGN = 10, despite the fact that N itself was
  a KIND = 8 integer!  
 
  Sequential calculation.
*/
    strcpy ( mode, "SEQ" );

    wtime = omp_get_wtime ( );

    estimate = r8_pi_est_seq ( n );

    wtime = omp_get_wtime ( ) - wtime;

    error = r8_abs ( estimate - r8_pi );

    printf ( "%14d  %s  %14f  %14g  %14f\n", n, mode, estimate, error, wtime );
/*
  Open MP enabled calculation.
*/
    strcpy ( mode, "OMP" );

    wtime = omp_get_wtime ( );

    estimate = r8_pi_est_omp ( n );

    wtime = omp_get_wtime ( ) - wtime;

    error = r8_abs ( estimate - r8_pi );

    printf ( "%14d  %s  %14f  %14g  %14f\n", n, mode, estimate, error, wtime );

    n = n * 10;
  }
//.........这里部分代码省略.........

void test03 ( )

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

    TEST03 tests P00_GAUSS_HERMITE.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 September 2012

  Author:

    John Burkardt
*/
{
  double error;
  double estimate;
  double exact;
  int m;
  int order;
  int order_log;
  int problem;
  int problem_num;

  printf ( "\n" );
  printf ( "TEST03\n" );
  printf ( "  P00_GAUSS_HERMITE applies a Gauss-Hermite rule\n" );
  printf ( "  to estimate an integral on (-oo,+oo).\n" );

  problem_num = p00_problem_num ( );

  m = 4;
  p06_param ( 'S', 'M', &m );

  printf ( "\n" );
  printf ( "   Problem     Order          Estimate        Exact          Error\n" );

  for ( problem = 1; problem <= problem_num; problem++ )
  {
    exact = p00_exact ( problem );

    order = 1;

    printf ( "\n" );

    for ( order_log = 0; order_log <= 6; order_log++ )
    {
      estimate = p00_gauss_hermite ( problem, order );

      error = r8_abs ( exact - estimate );

      printf ( "  %8d  %8d  %14.6g  %14.6g  %14.6g\n", 
        problem, order, estimate, exact, error );

      order = order * 2;
    }
  }
  return;
}

double class_matrix ( int kind, int m, double alpha, double beta, double aj[], 
  double bj[] )

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

    CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.

  Discussion:

    This routine computes the diagonal AJ and sub-diagonal BJ
    elements of the order M tridiagonal symmetric Jacobi matrix
    associated with the polynomials orthogonal with respect to
    the weight function specified by KIND.

    For weight functions 1-7, M elements are defined in BJ even
    though only M-1 are needed.  For weight function 8, BJ(M) is
    set to zero.

    The zero-th moment of the weight function is returned in ZEMU.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    11 January 2010

  Author:

    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
    C version by John Burkardt.

  Reference:

    Sylvan Elhay, Jaroslav Kautsky,
    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of 
    Interpolatory Quadrature,
    ACM Transactions on Mathematical Software,
    Volume 13, Number 4, December 1987, pages 399-415.

  Parameters:

    Input, int KIND, the rule.
    1, Legendre,             (a,b)       1.0
    2, Chebyshev Type 1,     (a,b)       ((b-x)*(x-a))^(-0.5)
    3, Gegenbauer,           (a,b)       ((b-x)*(x-a))^alpha
    4, Jacobi,               (a,b)       (b-x)^alpha*(x-a)^beta
    5, Generalized Laguerre, (a,inf)     (x-a)^alpha*exp(-b*(x-a))
    6, Generalized Hermite,  (-inf,inf)  |x-a|^alpha*exp(-b*(x-a)^2)
    7, Exponential,          (a,b)       |x-(a+b)/2.0|^alpha
    8, Rational,             (a,inf)     (x-a)^alpha*(x+b)^beta

    Input, int M, the order of the Jacobi matrix.

    Input, double ALPHA, the value of Alpha, if needed.

    Input, double BETA, the value of Beta, if needed.

    Output, double AJ[M], BJ[M], the diagonal and subdiagonal
    of the Jacobi matrix.

    Output, double CLASS_MATRIX, the zero-th moment.
*/
{
  double a2b2;
  double ab;
  double aba;
  double abi;
  double abj;
  double abti;
  double apone;
  int i;
  const double pi = 3.14159265358979323846264338327950;
  double temp;
  double temp2;
  double zemu;

  temp = r8_epsilon ( );

  parchk ( kind, 2 * m - 1, alpha, beta );

  temp2 = 0.5;

  if ( 500.0 * temp < r8_abs ( pow ( r8_gamma ( temp2 ), 2 ) - pi ) )
  {
    printf ( "\n" );
    printf ( "CLASS_MATRIX - Fatal error!\n" );
    printf ( "  Gamma function does not match machine parameters.\n" );
    exit ( 1 );
  }

  if ( kind == 1 )
  {
    ab = 0.0;

    zemu = 2.0 / ( ab + 1.0 );

//.........这里部分代码省略.........

//.........这里部分代码省略.........
  c = xx - 1.0;
/*
  Starting approximation for small chi-squared
*/
  if ( v < - c5 * log ( p ) )
  {
    ch = pow ( p * xx * exp ( g + xx * aa ), 1.0 / xx );

    if ( ch < e )
    {
      value = ch;
      return value;
    }
  }
/*
  Starting approximation for V less than or equal to 0.32
*/
  else if ( v <= c3 )
  {
    ch = c4;
    a = log ( 1.0 - p );

    for ( ; ; )
    {
      q = ch;
      p1 = 1.0 + ch * ( c7 + ch );
      p2 = ch * (c9 + ch * ( c8 + ch ) );

      t = - 0.5 + (c7 + 2.0 * ch ) / p1 - ( c9 + ch * ( c10 + 
        3.0 * ch ) ) / p2;

      ch = ch - ( 1.0 - exp ( a + g + 0.5 * ch + c * aa ) * p2 / p1) / t;

      if ( r8_abs ( q / ch - 1.0 ) <= c1 )
      {
        break;
      }
    }
  }
  else
  {
/*
  Call to algorithm AS 111 - note that P has been tested above.
  AS 241 could be used as an alternative.
*/
    x = ppnd ( p, ifault );
/*
  Starting approximation using Wilson and Hilferty estimate
*/
    p1 = c2 / v;
    ch = v * pow ( x * sqrt ( p1 ) + 1.0 - p1, 3 );
/*
  Starting approximation for P tending to 1.
*/
    if ( c6 * v + 6.0 < ch )
    {
      ch = - 2.0 * ( log ( 1.0 - p ) - c * log ( 0.5 * ch ) + g );
    }
  }
/*
  Call to algorithm AS 239 and calculation of seven term
  Taylor series
*/
  for ( i = 1; i <= maxit; i++ )
  {
    q = ch;

double *r8ge_fs_new ( int n, double a[], double b[] )

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

    R8GE_FS_NEW factors and solves a R8GE system.

  Discussion:

    The R8GE storage format is used for a "general" M by N matrix.  
    A physical storage space is made for each logical entry.  The two 
    dimensional logical array is mapped to a vector, in which storage is 
    by columns.

    The function does not save the LU factors of the matrix, and hence cannot
    be used to efficiently solve multiple linear systems, or even to
    factor A at one time, and solve a single linear system at a later time.

    The function uses partial pivoting, but no pivot vector is required.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    10 February 2012

  Author:

    John Burkardt

  Parameters:

    Input, int N, the order of the matrix.
    N must be positive.

    Input/output, double A[N*N].
    On input, A is the coefficient matrix of the linear system.
    On output, A is in unit upper triangular form, and
    represents the U factor of an LU factorization of the
    original coefficient matrix.

    Input, double B[N], the right hand side of the linear system.

    Output, double R8GE_FS_NEW[N], the solution of the linear system.
*/
{
  int i;
  int ipiv;
  int j;
  int jcol;
  double piv;
  double t;
  double *x;

  x = ( double * ) malloc ( n * sizeof ( double ) );

  for ( i = 0; i < n; i++ )
  {
    x[i] = b[i];
  }

  for ( jcol = 1; jcol <= n; jcol++ )
  {
/*
  Find the maximum element in column I.
*/
    piv = r8_abs ( a[jcol-1+(jcol-1)*n] );
    ipiv = jcol;
    for ( i = jcol+1; i <= n; i++ )
    {
      if ( piv < r8_abs ( a[i-1+(jcol-1)*n] ) )
      {
        piv = r8_abs ( a[i-1+(jcol-1)*n] );
        ipiv = i;
      }
    }

    if ( piv == 0.0 )
    {
      fprintf ( stderr, "\n" );
      fprintf ( stderr, "R8GE_FS_NEW - Fatal error!\n" );
      fprintf ( stderr, "  Zero pivot on step %d\n", jcol );
      exit ( 1 );
    }
/*
  Switch rows JCOL and IPIV, and X.
*/
    if ( jcol != ipiv )
    {
      for ( j = 1; j <= n; j++ )
      {
        t                 = a[jcol-1+(j-1)*n];
        a[jcol-1+(j-1)*n] = a[ipiv-1+(j-1)*n];
        a[ipiv-1+(j-1)*n] = t;
      }
      t         = x[jcol-1];
      x[jcol-1] = x[ipiv-1];
//.........这里部分代码省略.........

void test01 ( void )

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

    TEST01 demonstrates the use of XINBTA.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    05 November 2010

  Author:

    John Burkardt
*/
{
  double a;
  double b;
  double beta_log;
  double fx;
  int ifault;
  int n_data;
  double x;
  double x2;

  printf ( "\n" );
  printf ( "TEST01:\n" );
  printf ( "  XINBTA inverts the incomplete Beta function.\n" );
  printf ( "  Given CDF, it computes an X.\n" );
  printf ( "\n" );
  printf ( "           A           B           CDF    " );
  printf ( "    X                         X\n" );
  printf ( "                                          " );
  printf ( "    (Tabulated)               (XINBTA)            DIFF\n" );
  printf ( "\n" );

  n_data = 0;

  for ( ; ; )
  {
    beta_inc_values ( &n_data, &a, &b, &x, &fx );

    if ( n_data == 0 )
    {
      break;
    }

    beta_log = alngam ( a, &ifault )
             + alngam ( b, &ifault )
             - alngam ( a + b, &ifault );

    x2 = xinbta ( a, b, beta_log, fx, &ifault );

    printf ( "  %10.4f  %10.4f  %10.4f  %24.16g  %24.16g  %10.4e\n",
      a, b, fx, x, x2, r8_abs ( x - x2 ) );
  }

  return;
}

//.........这里部分代码省略.........
/*
  Use Pearson's series expansion.
  (Note that P is not large enough to force overflow in ALOGAM).
  No need to test IFAULT on exit since P > 0.
*/
  if ( x <= 1.0 || x < p )
  {
    arg = p * log ( x ) - x - alngam ( p + 1.0, ifault );
    c = 1.0;
    value = 1.0;
    a = p;

    for ( ; ; )
    {
      a = a + 1.0;
      c = c * x / a;
      value = value + c;

      if ( c <= tol )
      {
        break;
      }
    }

    arg = arg + log ( value );

    if ( elimit <= arg )
    {
      value = exp ( arg );
    }
    else
    {
      value = 0.0;
    }
  }
/*
  Use a continued fraction expansion.
*/
  else 
  {
    arg = p * log ( x ) - x - alngam ( p, ifault );
    a = 1.0 - p;
    b = a + x + 1.0;
    c = 0.0;
    pn1 = 1.0;
    pn2 = x;
    pn3 = x + 1.0;
    pn4 = x * b;
    value = pn3 / pn4;

    for ( ; ; )
    {
      a = a + 1.0;
      b = b + 2.0;
      c = c + 1.0;
      an = a * c;
      pn5 = b * pn3 - an * pn1;
      pn6 = b * pn4 - an * pn2;

      if ( pn6 != 0.0 )
      {
        rn = pn5 / pn6;

        if ( r8_abs ( value - rn ) <= r8_min ( tol, tol * rn ) )
        {
          break;
        }
        value = rn;
      }

      pn1 = pn3;
      pn2 = pn4;
      pn3 = pn5;
      pn4 = pn6;
/*
  Re-scale terms in continued fraction if terms are large.
*/
      if ( oflo <= abs ( pn5 ) )
      {
        pn1 = pn1 / oflo;
        pn2 = pn2 / oflo;
        pn3 = pn3 / oflo;
        pn4 = pn4 / oflo;
      }
    }

    arg = arg + log ( value );

    if ( elimit <= arg )
    {
      value = 1.0 - exp ( arg );
    }
    else
    {
      value = 1.0;
    }
  }

  return value;
}

//.........这里部分代码省略.........
  double g2 = 0.1591549431;
  double g3 = 2.5066282746;
  int ioe;
  int k;
  double rb;
  double sum;
  double value;

  f = ( double ) ( idf );
/*
  For very large IDF, use the normal approximation.
*/
  if ( 100 < idf )
  {
    *ifault = 1;

    a = sqrt ( 0.5 * f ) 
    * exp ( alngam ( 0.5 * ( f - 1.0 ), &k ) 
    - alngam ( 0.5 * f, &k ) ) * d;

    value = alnorm ( ( st - a ) / sqrt ( f * ( 1.0 + d * d ) 
    / ( f - 2.0 ) - a * a ), 0 );
    return value;
  }

  *ifault = 0;
  ioe = ( idf % 2 );
  a = st / sqrt ( f );
  b = f / ( f + st * st );
  rb = sqrt ( b );
  da = d * a;
  drb = d * rb;

  if ( idf == 1 )
  {
    value = alnorm ( drb, 1 ) + 2.0 * tfn ( drb, a );
    return value;
  }

  sum = 0.0;

  if ( r8_abs ( drb ) < emin )
  {
    fmkm2 = a * rb * exp ( - 0.5 * drb * drb ) 
    * alnorm ( a * drb, 0 ) * g1;
  }
  else
  {
    fmkm2 = 0.0;
  }

  fmkm1 = b * da * fmkm2;

  if ( r8_abs ( d ) < emin )
  {
    fmkm1 = fmkm1 + b * a * g2 * exp ( - 0.5 * d * d );
  }

  if ( ioe == 0 )
  {
    sum = fmkm2;
  }
  else
  {
    sum = fmkm1;
  }

  ak = 1.0;
  fk = 2.0;

  for ( k = 2; k <= idf - 2; k = k + 2 )
  {
    fkm1 = fk - 1.0;
    fmkm2 = b * ( da * ak * fmkm1 + fmkm2 ) * fkm1 / fk;
    ak = 1.0 / ( ak * fkm1 );
    fmkm1 = b * ( da * ak * fmkm2 + fmkm1 ) * fk / ( fk + 1.0 );

    if ( ioe == 0 )
    {
      sum = sum + fmkm2;
    }
    else
    {
      sum = sum + fmkm1;
    }
    ak = 1.0 / ( ak * fk );
    fk = fk + 2.0;
  }

  if ( ioe == 0 )
  {
    value = alnorm ( d, 1 ) + sum * g3;
  }
  else
  {
    value = alnorm ( drb, 1 ) + 2.0 * ( sum + tfn ( drb, a ) );
  }

  return value;
}

double ppnd ( double p, int *ifault )

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

    PPND produces the normal deviate value corresponding to lower tail area = P.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    03 November 2010

  Author:

    Original FORTRAN77 version by J Beasley, S Springer.
    C version by John Burkardt.

  Reference:

    J Beasley, S Springer,
    Algorithm AS 111:
    The Percentage Points of the Normal Distribution,
    Applied Statistics,
    Volume 26, Number 1, 1977, pages 118-121.

  Parameters:

    Input, double P, the value of the cumulative probability
    densitity function.  0 < P < 1.

    Output, integer *IFAULT, error flag.
    0, no error.
    1, P <= 0 or P >= 1.  PPND is returned as 0.

    Output, double PPND, the normal deviate value with the property that
    the probability of a standard normal deviate being less than or
    equal to PPND is P.
*/
{
  double a0 = 2.50662823884;
  double a1 = -18.61500062529;
  double a2 = 41.39119773534;
  double a3 = -25.44106049637;
  double b1 = -8.47351093090;
  double b2 = 23.08336743743;
  double b3 = -21.06224101826;
  double b4 = 3.13082909833;
  double c0 = -2.78718931138;
  double c1 = -2.29796479134;
  double c2 = 4.85014127135;
  double c3 = 2.32121276858;
  double d1 = 3.54388924762;
  double d2 = 1.63706781897;
  double r;
  double split = 0.42;
  double value;

  *ifault = 0;
/*
  0.08 < P < 0.92
*/
  if ( r8_abs ( p - 0.5 ) <= split )
  {
    r = ( p - 0.5 ) * ( p - 0.5 );

    value = ( p - 0.5 ) * ( ( ( 
        a3   * r 
      + a2 ) * r 
      + a1 ) * r 
      + a0 ) / ( ( ( ( 
        b4   * r 
      + b3 ) * r 
      + b2 ) * r 
      + b1 ) * r 
      + 1.0 );
  }
/*
  P < 0.08 or P > 0.92,
  R = min ( P, 1-P )
*/
  else if ( 0.0 < p && p < 1.0 )
  {
    if ( 0.5 < p )
    {
      r = sqrt ( - log ( 1.0 - p ) );
    }
    else
    {
      r = sqrt ( - log ( p ) );
    }

    value = ( ( ( 
        c3   * r 
      + c2 ) * r 
      + c1 ) * r 
      + c0 ) / ( ( 
//.........这里部分代码省略.........

double tfn ( double x, double fx )

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

    TFN calculates the T-function of Owen.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    16 January 2008

  Author:

    Original FORTRAN77 version by JC Young, Christoph Minder.
    C version by John Burkardt.

  Reference:

    MA Porter, DJ Winstanley,
    Remark AS R30:
    A Remark on Algorithm AS76:
    An Integral Useful in Calculating Noncentral T and Bivariate
    Normal Probabilities,
    Applied Statistics,
    Volume 28, Number 1, 1979, page 113.

    JC Young, Christoph Minder,
    Algorithm AS 76: 
    An Algorithm Useful in Calculating Non-Central T and 
    Bivariate Normal Distributions,
    Applied Statistics,
    Volume 23, Number 3, 1974, pages 455-457.

  Parameters:

    Input, double X, FX, the parameters of the function.

    Output, double TFN, the value of the T-function.
*/
{
# define NG 5

  double fxs;
  int i;
  double r[NG] = {
    0.1477621, 
    0.1346334, 
    0.1095432, 
    0.0747257, 
    0.0333357 };
  double r1;
  double r2;
  double rt;
  double tp = 0.159155;
  double tv1 = 1.0E-35;
  double tv2 = 15.0;
  double tv3 = 15.0;
  double tv4 = 1.0E-05;
  double u[NG] = {
    0.0744372, 
    0.2166977, 
    0.3397048,
    0.4325317, 
    0.4869533 };
  double value;
  double x1;
  double x2;
  double xs;
/*
  Test for X near zero.
*/
  if ( r8_abs ( x ) < tv1 )
  {
    value = tp * atan ( fx );
    return value;
  }
/*
  Test for large values of abs(X).
*/
  if ( tv2 < r8_abs ( x ) )
  {
    value = 0.0;
    return value;
  }
/*
  Test for FX near zero.
*/
  if ( r8_abs ( fx ) < tv1 )
  {
    value = 0.0;
    return value;
  }
/*
  Test whether abs ( FX ) is so large that it must be truncated.
*/
//.........这里部分代码省略.........