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

using namespace std;

# include "polpak.H"

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

double agm ( double a, double b )

//****************************************************************************80
//
//  Purpose:
//
//    AGM computes the arithmetic-geometric mean of A and B.
//
//  Discussion:
//
//    The AGM is defined for nonnegative A and B.
//
//    The AGM of numbers A and B is defined by setting
//
//      A(0) = A,
//      B(0) = B
//
//      A(N+1) = ( A(N) + B(N) ) / 2
//      B(N+1) = sqrt ( A(N) * B(N) )
//
//    The two sequences both converge to AGM(A,B).
//
//    In Mathematica, the AGM can be evaluated by
//
//      ArithmeticGeometricMean [ a, b ]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 February 2008
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Parameters:
//
//    Input, double A, B, the arguments whose AGM is to be computed.
//
//    Output, double AGM, the arithmetic-geometric mean of A and B.
//
{
  double c;
  double d;
  int it;
  int it_max = 1000;
  double tol;
  double value;

  if ( a < 0.0 )
  {
    exit ( 1 );
  }

  if ( b < 0.0 )
  {
    exit ( 1 );
  }

  if ( a == 0.0 || b == 0.0 )
  {
    value = 0.0;
    return value;
  }

  it = 0;
  tol = 100.0 * r8_epsilon ( );

  for ( ; ; )
  {
    it = it + 1;

    c = ( a + b ) / 2.0;
    d = sqrt ( a * b );

    if ( r8_abs ( c - d ) <= tol * ( c + d ) )
    {
      break;
    }

    if ( it_max < it )
    {
      break;
    }

    a = c;
    b = d;
  }
  value = c;

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

void agm_values ( int *n_data, double *a, double *b, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    AGM_VALUES returns some values of the AGM.
//
//  Discussion:
//
//    The AGM is defined for nonnegative A and B.
//
//    The AGM of numbers A and B is defined by setting
//
//      A(0) = A,
//      B(0) = B
//
//      A(N+1) = ( A(N) + B(N) ) / 2
//      B(N+1) = sqrt ( A(N) * B(N) )
//
//    The two sequences both converge to AGM(A,B).
//
//    In Mathematica, the AGM can be evaluated by
//
//      ArithmeticGeometricMean [ a, b ]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 February 2008
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Parameters:
//
//    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, double *A, *B, the argument ofs the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 15

  double a_vec[N_MAX] = {
        22.0, 
        83.0, 
        42.0, 
        26.0, 
         4.0, 
         6.0, 
        40.0,
        80.0,
        90.0,
         9.0,
        53.0,
         1.0,
         1.0,
         1.0,
         1.5 }; 
  double b_vec[N_MAX] = {
        96.0,
        56.0,
         7.0,
        11.0,
        63.0,
        45.0,
        75.0,
         0.0,
        35.0,
         1.0,
        53.0,
         2.0,
         4.0,
         8.0,
         8.0 };    
  double fx_vec[N_MAX] = {
        52.274641198704240049,
        68.836530059858524345,
        20.659301196734009322,
        17.696854873743648823,
        23.867049721753300163,
        20.717015982805991662,
        56.127842255616681863,
         0.000000000000000000,
        59.269565081229636528,
        3.9362355036495554780,
        53.000000000000000000,
        1.4567910310469068692,
        2.2430285802876025701,
        3.6157561775973627487,
        4.0816924080221632670 };
 
  if ( *n_data < 0 )
  {
    *n_data = 0;
  }

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *a = 0.0;
    *b = 0.0;
    *fx = 0.0;
  }
  else
  {
    *a = a_vec[*n_data-1];
    *b = b_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

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

double agud ( double gamma )

//****************************************************************************80
//
//  Purpose:
//
//    AGUD evaluates the inverse Gudermannian function.
//
//  Definition:
//
//    The Gudermannian function relates the hyperbolic and trigonomentric
//    functions.  For any argument X, there is a corresponding value
//    GAMMA so that
//
//      SINH(X) = TAN(GAMMA).
//
//    This value GAMMA(X) is called the Gudermannian of X.  The inverse
//    Gudermannian function is given as input a value GAMMA and computes
//    the corresponding value X.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double GAMMA, the value of the Gudermannian.
//
//    Output, double AGUD, the argument of the Gudermannian.
//
{
  double value;

  value = log ( tan ( 0.25 * r8_pi ( ) + 0.5 * gamma ) );

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

int align_enum ( int m, int n )

//****************************************************************************80
//
//  Purpose:
//
//    ALIGN_ENUM counts the alignments of two sequences of M and N elements.
//
//  Discussion:
//
//    We assume that we have sequences A and B of M and N characters each.
//    An alignment of the two sequences is a rule matching corresponding
//    elements of one sequence to another.  Some elements of either sequence
//    can be matched to a null element.  If A(I1) and A(I2) are matched
//    to B(J1) and B(J2), and I1 < I2, then it must be the case that J1 < J2.
//
//    The 5 alignments of a sequence of 1 to a sequence of 2 are:
//
//          _1_   _2_   __3__   __4__   __5__
//
//      A:  1 -   - 1   - 1 -   - - 1   1 - -
//      B:  1 2   1 2   1 - 2   1 2 -   - 1 2
//
//    The formula is:
//
//      F(0,0) = 1
//      F(1,0) = 1
//      F(0,1) = 1
//      F(M,N) = F(M-1,N) + F(M-1,N-1) + F(M,N-1)
//
//    To compute F(M,N), it is not necessary to keep an M+1 by N+1
//    array in memory.  A vector of length N will do.
//
//    F(N,N) is approximately ( 1 + sqrt(2) )**(2*N+1) / sqrt ( N )
//
//  Example:
//
//    The initial portion of the table is:
//
//
//  M/N   0    1    2    3    4       5       6       7       8       9      10
//
//   0    1    1    1    1    1       1       1       1       1       1       1
//   1    1    3    5    7    9      11      13      15      17      19      21
//   2    1    5   13   25   41      61      85     113     145     181     221
//   3    1    7   25   63  129     231     377     575     833    1159    1561
//   4    1    9   41  129  321     681    1289    2241    3649    5641    8361
//   5    1   11   61  231  681    1683    3653    7183   13073   22363   36365
//   6    1   13   85  377 1289    3653    8989   19825   40081   75517  134245
//   7    1   15  113  575 2241    7183   19825   48639  108545  224143  433905
//   8    1   17  145  833 3649   13073   40081  108545  265729  598417 1256465
//   9    1   19  181 1159 5641   22363   75517  224143  598417 1462563 3317445
//  10    1   21  221 1561 8361   36365  134245  433905 1256465 3317445 8097453
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Michael Waterman,
//    Introduction to Computational Biology,
//    Chapman and Hall, 1995, pages 186-190.
//
//  Parameters:
//
//    Input, int M, N, the number of elements of the two sequences.
//
//    Output, int ALIGN_ENUM, the number of possible alignments of the
//    sequences.
//
{
  int *fi;
  int fim1j;
  int fim1jm1;
  int i;
  int j;
  int value;

  if ( m < 0 )
  {
    return 0;
  }
  else if ( n < 0 )
  {
    return 0;
  }
  else if ( m == 0 )
  {
    return 1;
  }
  else if ( n == 0 )
  {
    return 1;
  }

  fi = new int[n+1];
  for ( i = 0; i <= n; i++ )
  {
    fi[i] = 1;
  }

  for ( i = 1; i <= m; i++ )
  {
    fim1jm1 = 1;

    for ( j = 1; j <= n; j++ )
    {

      fim1j = fi[j];

      fi[j] = fi[j] + fi[j-1] + fim1jm1;

      fim1jm1 = fim1j;

    }
  }

  value = fi[n];
  delete [] fi;

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

double arc_cosine ( double c )

//****************************************************************************80
//
//  Purpose:
//
//    ARC_COSINE computes the arc cosine function, with argument truncation.
//
//  Discussion:
//
//    If you call your system ACOS routine with an input argument that is
//    outside the range [-1.0, 1.0 ], you may get an unpleasant surprise.
//
//    In particular, you may get the value NaN returned.
//
//    This routine truncates arguments outside the range, avoiding the problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double C, the argument.
//
//    Output, double ARC_COSINE, an angle whose cosine is C.
//
{
  double value;

  c = r8_max ( c, -1.0 );
  c = r8_min ( c, +1.0 );

  value = acos ( c );

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

double arc_sine ( double s )

//****************************************************************************80
//
//  Purpose:
//
//    ARC_SINE computes the arc sine function, with argument truncation.
//
//  Discussion:
//
//    If you call your system ASIN routine with an input argument that is
//    outside the range [-1.0, 1.0 ], you may get an unpleasant surprise.
//
//    In particular, you may get the value NaN returned.
//
//    This routine truncates arguments outside the range, avoiding the problem.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    04 March 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double S, the argument.
//
//    Output, double ARC_SINE, an angle whose sine is S.
//
{
  double value;

  s = r8_max ( s, -1.0 );
  s = r8_min ( s, +1.0 );

  value = asin ( s );

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

double asinh2 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    ASINH2 returns the inverse hyperbolic sine of a number.
//
//  Definition:
//
//    The assertion that:
//
//      Y = ASINH2(X)
//
//    implies that
//
//      X = SINH(Y) = 0.5 * ( EXP(Y) - EXP(-Y) ).
//
//  Discussion:
//
//    Since a library function ASINH may be available on some systems,
//    this routine is named ASINH2 to avoid name conflicts.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the number whose inverse hyperbolic sine is desired.
//
//    Output, double ASINH2, the inverse hyperbolic sine of X.
//
{
  double value;

  value = log ( x + sqrt ( x * x + 1.0 ) );

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

double atan4 ( double y, double x )

//****************************************************************************80
//
//  Purpose:
//
//    ATAN4 computes the inverse tangent of the ratio Y / X.
//
//  Discussion:
//
//    ATAN4 returns an angle whose tangent is ( Y / X ), a job which
//    the built in functions ATAN and ATAN2 already do.
//
//    However:
//
//    * ATAN4 always returns a positive angle, between 0 and 2 PI,
//      while ATAN and ATAN2 return angles in the interval [-PI/2,+PI/2]
//      and [-PI,+PI] respectively;
//
//    * ATAN4 accounts for the signs of X and Y, (as does ATAN2).  The ATAN
//     function by contrast always returns an angle in the first or fourth
//     quadrants.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double Y, X, two quantities which represent the tangent of
//    an angle.  If Y is not zero, then the tangent is (Y/X).
//
//    Output, double ATAN4, an angle between 0 and 2 * PI, whose tangent is
//    (Y/X), and which lies in the appropriate quadrant so that the signs
//    of its cosine and sine match those of X and Y.
//
{
  double abs_x;
  double abs_y;
  double theta;
  double theta_0;
//
//  Special cases:
//
  if ( x == 0.0 )
  {
    if ( 0.0 < y )
    {
      theta = r8_pi ( ) / 2.0;
    }
    else if ( y < 0.0 )
    {
      theta = 3.0 * r8_pi ( ) / 2.0;
    }
    else if ( y == 0.0 )
    {
      theta = 0.0;
    }
  }
  else if ( y == 0.0 )
  {
    if ( 0.0 < x )
    {
      theta = 0.0;
    }
    else if ( x < 0.0 )
    {
      theta = r8_pi ( );
    }
  }
//
//  We assume that ATAN2 is correct when both arguments are positive.
//
  else
  {
    abs_y = fabs ( y );
    abs_x = fabs ( x );

    theta_0 = atan2 ( abs_y, abs_x );

    if ( 0.0 < x && 0.0 < y )
    {
      theta = theta_0;
    }
    else if ( x < 0.0 && 0.0 < y )
    {
      theta = r8_pi ( ) - theta_0;
    }
    else if ( x < 0.0 && y < 0.0 )
    {
      theta = r8_pi ( ) + theta_0;
    }
    else if ( 0.0 < x && y < 0.0 )
    {
      theta = 2.0 * r8_pi ( ) - theta_0;
    }

  }

  return theta;
}
//****************************************************************************80

double atanh2 ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    ATANH2 returns the inverse hyperbolic tangent of a number.
//
//  Definition:
//
//    Y = ATANH2(X) implies that
//    X = TANH(Y) = ( EXP(Y) - EXP(-Y) ) / ( EXP(Y) + EXP(-Y) )
//
//  Discussion:
//
//    Since a library function ATANH may be available on some systems,
//    this routine is named ATANH2 to avoid name conflicts.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the number whose inverse hyperbolic tangent is desired.
//    The absolute value of X should be less than or equal to 1.
//
//    Output, double ATANH2, the inverse hyperbolic tangent of X.
//
{
  double value;

  if ( x <= -1.0 )
  {
    value = -HUGE_VAL;
  }
  else if ( 1.0 <= x )
  {
    value = HUGE_VAL;
  }
  else
  {
    value = 0.5 * log ( ( 1.0 + x ) / ( 1.0 - x ) );
  }
  return value;
}
//****************************************************************************80

void bell ( int n, int b[] )

//****************************************************************************80
//
//  Purpose:
//
//    BELL returns the Bell numbers from 0 to N.
//
//  Discussion:
//
//    The Bell number B(N) is the number of restricted growth functions
//    on N.
//
//    Note that the Stirling numbers of the second kind, S^m_n, count the
//    number of partitions of N objects into M classes, and so it is
//    true that
//
//      B(N) = S^1_N + S^2_N + ... + S^N_N.
//
//  Definition:
//
//    The Bell number B(N) is defined as the number of partitions (of
//    any size) of a set of N distinguishable objects.
//
//    A partition of a set is a division of the objects of the set into
//    subsets.
//
//  Examples:
//
//    There are 15 partitions of a set of 4 objects:
//
//      (1234), (123)(4), (124)(3), (12)(34), (12)(3)(4),
//      (134)(2), (13)(24), (13)(2)(4), (14)(23), (1)(234),
//      (1)(23)(4), (14)(2)(3), (1)(24)(3), (1)(2)(34), (1)(2)(3)(4)
//
//    and so B(4) = 15.
//
//  First values:
//
//     N         B(N)
//     0           1
//     1           1
//     2           2
//     3           5
//     4          15
//     5          52
//     6         203
//     7         877
//     8        4140
//     9       21147
//    10      115975
//
//  Recursion:
//
//    B(I) = sum ( 1 <= J <= I ) Binomial ( I-1, J-1 ) * B(I-J)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of Bell numbers desired.
//
//    Output, int B[N+1], the Bell numbers from 0 to N.
//
{
  int i;
  int j;

  b[0] = 1;

  for ( i = 1; i <= n; i++ )
  {
    b[i] = 0;
    for ( j = 1; j <= i; j++ )
    {
      b[i] = b[i] + b[i-j] * i4_choose ( i-1, j-1 );
    }
  }

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

void bell_values ( int *n_data, int *n, int *c )

//****************************************************************************80
//
//  Purpose:
//
//    BELL_VALUES returns some values of the Bell numbers.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and N_DATA
//    is set to 1.  On each subsequent call, the input value of N_DATA is
//    incremented and that test data item is returned, if available.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, int *N, the order of the Bell number.
//
//    Output, int *C, the value of the Bell number.
//
{
# define N_MAX 11

  int c_vec[N_MAX] = { 1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 };
  int n_vec[N_MAX] = { 0,  1,  2,  3,  4, 5,  6,  7,  8,  9,  10};

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *n = 0;
    *c = 0;
  }
  else
  {
    *n = n_vec[*n_data];
    *c = c_vec[*n_data];
    *n_data = *n_data + 1;
  }

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

double benford ( int ival )

//****************************************************************************80
//
//  Purpose:
//
//    BENFORD returns the Benford probability of one or more significant digits.
//
//  Discussion:
//
//    Benford's law is an empirical formula explaining the observed
//    distribution of initial digits in lists culled from newspapers,
//    tax forms, stock market prices, and so on.  It predicts the observed
//    high frequency of the initial digit 1, for instance.
//
//    Note that the probabilities of digits 1 through 9 are guaranteed
//    to add up to 1, since
//      LOG10 ( 2/1 ) + LOG10 ( 3/2) + LOG10 ( 4/3 ) + ... + LOG10 ( 10/9 )
//      = LOG10 ( 2/1 * 3/2 * 4/3 * ... * 10/9 ) = LOG10 ( 10 ) = 1.
//
//  Formula:
//
//    Prob ( First significant digits are IVAL ) =
//      LOG10 ( ( IVAL + 1 ) / IVAL ).
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Frank Benford,
//    The Law of Anomalous Numbers,
//    Proceedings of the American Philosophical Society,
//    Volume 78, pages 551-572, 1938.
//
//    T P Hill,
//    The First Digit Phenomenon,
//    American Scientist,
//    Volume 86, July/August 1998, pages 358 - 363.
//
//    R Raimi,
//    The Peculiar Distribution of First Digits,
//    Scientific American,
//    December 1969, pages 109-119.
//
//  Parameters:
//
//    Input, int IVAL, the string of significant digits to be checked.
//    If IVAL is 1, then we are asking for the Benford probability that
//    a value will have first digit 1.  If IVAL is 123, we are asking for
//    the probability that the first three digits will be 123, and so on.
//    Note that IVAL must not be 0 or negative.
//
//    Output, double BENFORD, the Benford probability that an item taken
//    from a real world distribution will have the initial digits IVAL.
//
{
  double value;

  if ( ival <= 0 )
  {
    cout << "\n";
    cout << "BENFORD - Fatal error!\n";
    cout << "  The input argument must be positive.\n";
    cout << "  Your value was " << ival << "\n";
    exit ( 1 );
  }

  value = log10 ( ( double ) ( ival + 1 ) / ( double ) ival  );

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

void bernoulli_number ( int n, double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    BERNOULLI_NUMBER computes the value of the Bernoulli numbers B(0) through B(N).
//
//  Discussion:
//
//    The Bernoulli numbers are rational.
//
//    If we define the sum of the M-th powers of the first N integers as:
//
//      SIGMA(M,N) = sum ( 0 <= I <= N ) I**M
//
//    and let C(I,J) be the combinatorial coefficient:
//
//      C(I,J) = I! / ( ( I - J )! * J! )
//
//    then the Bernoulli numbers B(J) satisfy:
//
//      SIGMA(M,N) = 1/(M+1) * sum ( 0 <= J <= M ) C(M+1,J) B(J) * (N+1)**(M+1-J)
//
//  First values:
//
//   B0  1                   =         1.00000000000
//   B1 -1/2                 =        -0.50000000000
//   B2  1/6                 =         1.66666666666
//   B3  0                   =         0
//   B4 -1/30                =        -0.03333333333
//   B5  0                   =         0
//   B6  1/42                =         0.02380952380
//   B7  0                   =         0
//   B8 -1/30                =        -0.03333333333
//   B9  0                   =         0
//  B10  5/66                =         0.07575757575
//  B11  0                   =         0
//  B12 -691/2730            =        -0.25311355311
//  B13  0                   =         0
//  B14  7/6                 =         1.16666666666
//  B15  0                   =         0
//  B16 -3617/510            =        -7.09215686274
//  B17  0                   =         0
//  B18  43867/798           =        54.97117794486
//  B19  0                   =         0
//  B20 -174611/330          =      -529.12424242424
//  B21  0                   =         0
//  B22  854,513/138         =      6192.123
//  B23  0                   =         0
//  B24 -236364091/2730      =    -86580.257
//  B25  0                   =         0
//  B26  8553103/6           =   1425517.16666
//  B27  0                   =         0
//  B28 -23749461029/870     = -27298231.0678
//  B29  0                   =         0
//  B30  8615841276005/14322 = 601580873.901
//
//  Recursion:
//
//    With C(N+1,K) denoting the standard binomial coefficient,
//
//    B(0) = 1.0
//    B(N) = - ( sum ( 0 <= K < N ) C(N+1,K) * B(K) ) / C(N+1,N)
//
//  Warning:
//
//    This recursion, which is used in this routine, rapidly results
//    in significant errors.
//
//  Special Values:
//
//    Except for B(1), all Bernoulli numbers of odd index are 0.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the highest Bernoulli number to compute.
//
//    Output, double B[N+1], B(I) contains the I-th Bernoulli number.
//
{
  double b_sum;
  int *combo;
  int i;
  int j;
  bool next;

  if ( n < 0 )
  {
    return;
  }

  b[0] = 1.0;

  if ( n < 1 )
  {
    return;
  }

  b[1] = -0.5;

  next = false;

  combo = new int[n+2];

  for ( i = 2; i <= n; i++ )
  {
    comb_row ( next, i+1, combo );
    next = true;

    if ( ( i % 2 ) == 1 )
    {
      b[i] = 0.0;
    }
    else
    {
      b_sum = 0.0;
      for ( j = 0; j <= i-1; j++ )
      {
        b_sum = b_sum + b[j] * ( double ) ( combo[j] );
      }

      b[i] = - b_sum / ( double ) ( combo[i] );

    }

  }

  delete [] combo;

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

void bernoulli_number2 ( int n, double b[] )

//****************************************************************************80
//
//  Purpose:
//
//    BERNOULLI_NUMBER2 evaluates the Bernoulli numbers.
//
//  Discussion:
//
//    The Bernoulli numbers are rational.
//
//    If we define the sum of the M-th powers of the first N integers as:
//
//      SIGMA(M,N) = sum ( 0 <= I <= N ) I**M
//
//    and let C(I,J) be the combinatorial coefficient:
//
//      C(I,J) = I! / ( ( I - J )! * J! )
//
//    then the Bernoulli numbers B(J) satisfy:
//
//      SIGMA(M,N) = 1/(M+1) * sum ( 0 <= J <= M ) C(M+1,J) B(J) * (N+1)**(M+1-J)
//
//    Note that the Bernoulli numbers grow rapidly.  Bernoulli number
//    62 is probably the last that can be computed on the VAX without
//    overflow.
//
//    A different method than that used in BERN is employed.
//
//  First values:
//
//   B0  1                   =         1.00000000000
//   B1 -1/2                 =        -0.50000000000
//   B2  1/6                 =         1.66666666666
//   B3  0                   =         0
//   B4 -1/30                =        -0.03333333333
//   B5  0                   =         0
//   B6  1/42                =         0.02380952380
//   B7  0                   =         0
//   B8 -1/30                =        -0.03333333333
//   B9  0                   =         0
//  B10  5/66                =         0.07575757575
//  B11  0                   =         0
//  B12 -691/2730            =        -0.25311355311
//  B13  0                   =         0
//  B14  7/6                 =         1.16666666666
//  B15  0                   =         0
//  B16 -3617/510            =        -7.09215686274
//  B17  0                   =         0
//  B18  43867/798           =        54.97117794486
//  B19  0                   =         0
//  B20 -174611/330          =      -529.12424242424
//  B21  0                   =         0
//  B22  854,513/138         =      6192.123
//  B23  0                   =         0
//  B24 -236364091/2730      =    -86580.257
//  B25  0                   =         0
//  B26  8553103/6           =   1425517.16666
//  B27  0                   =         0
//  B28 -23749461029/870     = -27298231.0678
//  B29  0                   =         0
//  B30  8615841276005/14322 = 601580873.901
//
//  Recursion:
//
//    With C(N+1,K) denoting the standard binomial coefficient,
//
//    B(0) = 1.0
//    B(N) = - ( sum ( 0 <= K < N ) C(N+1,K) * B(K) ) / C(N+1,N)
//
//  Special Values:
//
//    Except for B(1), all Bernoulli numbers of odd index are 0.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the highest order Bernoulli number to compute.
//
//    Output, double B[N+1], the requested Bernoulli numbers.
//
{
  double altpi;
  int i;
  int k;
  int kmax = 400;
  double sgn;
  double sum2;
  double t;
  double term;
  double tol = 1.0E-06;

  if ( n < 0 )
  {
    return;
  }

  b[0] = 1.0;

  if ( n < 1 )
  {
    return;
  }

  b[1] = -0.5;

  if ( n < 2 )
  {
    return;
  }

  altpi = log ( 2.0 * r8_pi ( ) );
//
//  Initial estimates for B(I), I = 2 to N
//
  b[2] = log ( 2.0 );

  for ( i = 3; i <= n; i++ )
  {
    if ( ( i % 2 ) == 1 )
    {
      b[i] = 0.0;
    }
    else
    {
      b[i] = log ( ( double ) ( i * ( i - 1 ) ) ) + b[i-2];
    }
  }

  b[2] = 1.0 / 6.0;

  if ( n <= 3 )
  {
    return;
  }

  b[4] = -1.0 / 30.0;

  sgn = -1.0;

  for ( i = 6; i <= n; i = i + 2 )
  {
    sgn = -sgn;
    t = 2.0 * sgn * exp ( b[i] - ( double ) ( i ) * altpi );

    sum2 = 1.0;

    for ( k = 2; k <= kmax; k++ )
    {
      term = pow ( ( double ) k, ( double ) -i );

      sum2 = sum2 + term;

      if ( term <= tol * sum2 )
      {
        break;
      }

    }

    b[i] = t * sum2;

  }

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

double bernoulli_number3 ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    BERNOULLI_NUMBER3 computes the value of the Bernoulli number B(N).
//
//  Discussion:
//
//    The Bernoulli numbers are rational.
//
//    If we define the sum of the M-th powers of the first N integers as:
//
//      SIGMA(M,N) = sum ( 0 <= I <= N ) I**M
//
//    and let C(I,J) be the combinatorial coefficient:
//
//      C(I,J) = I! / ( ( I - J )! * J! )
//
//    then the Bernoulli numbers B(J) satisfy:
//
//      SIGMA(M,N) = 1/(M+1) * sum ( 0 <= J <= M ) 
//        C(M+1,J) B(J) * (N+1)**(M+1-J)
//
//  First values:
//
//     B0  1                   =         1.00000000000
//     B1 -1/2                 =        -0.50000000000
//     B2  1/6                 =         1.66666666666
//     B3  0                   =         0
//     B4 -1/30                =        -0.03333333333
//     B5  0                   =         0
//     B6  1/42                =         0.02380952380
//     B7  0                   =         0
//     B8 -1/30                =        -0.03333333333
//     B9  0                   =         0
//    B10  5/66                =         0.07575757575
//    B11  0                   =         0
//    B12 -691/2730            =        -0.25311355311
//    B13  0                   =         0
//    B14  7/6                 =         1.16666666666
//    B15  0                   =         0
//    B16 -3617/510            =        -7.09215686274
//    B17  0                   =         0
//    B18  43867/798           =        54.97117794486
//    B19  0                   =         0
//    B20 -174611/330          =      -529.12424242424
//    B21  0                   =         0
//    B22  854513/138          =      6192.123
//    B23  0                   =         0
//    B24 -236364091/2730      =    -86580.257
//    B25  0                   =         0
//    B26  8553103/6           =   1425517.16666
//    B27  0                   =         0
//    B28 -23749461029/870     = -27298231.0678
//    B29  0                   =         0
//    B30  8615841276005/14322 = 601580873.901
//
//  Recursion:
//
//    With C(N+1,K) denoting the standard binomial coefficient,
//
//    B(0) = 1.0
//    B(N) = - ( sum ( 0 <= K < N ) C(N+1,K) * B(K) ) / C(N+1,N)
//
//  Special Values:
//
//    Except for B(1), all Bernoulli numbers of odd index are 0.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the Bernoulli number to compute.
//
//    Output, double BERNOULLI_NUMBER3, the desired Bernoulli number.
//
{
  int i;
  int itmax = 1000;
  double sum2;
  double term;
  double tol = 5.0E-07;
  double value;

  if ( n < 0 )
  {
    value = 0.0;
  }
  else if ( n == 0 )
  {
    value = 1.0;
  }
  else if ( n == 1 )
  {
    value = -0.5;
  }
  else if ( n == 2 )
  {
    value = 1.0 / 6.0;
  }
  else if ( ( n % 2 ) == 1 )
  {
    value = 0.0;
  }
  else
  {
    sum2 = 0.0;

    for ( i = 1; i <= itmax; i++ )
    {
      term = 1.0 / pow ( ( double ) i, n );
      sum2 = sum2 + term;

      if ( fabs ( term ) < tol || fabs ( term ) < tol * fabs ( sum2 ) )
      {
        break;
      }

    }

    value = 2.0 * sum2 * r8_factorial ( n ) 
      / pow ( ( 2.0 * r8_pi ( ) ), n );

    if ( ( n % 4 ) == 0 )
    {
      value = -value;
    }

  }

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

double bernoulli_poly ( int n, double x )

//****************************************************************************80
//
//  Purpose:
//
//    BERNOULLI_POLY evaluates the Bernoulli polynomial of order N at X.
//
//  Discussion:
//
//    Thanks to Bart Vandewoestyne for pointing out an error in the previous
//    documentation, 31 January 2008.
//
//    Special values of the Bernoulli polynomial include:
//
//      B(N,0) = B(N,1) = B(N), the N-th Bernoulli number.
//
//      B'(N,X) = N * B(N-1,X)
//
//      B(N,X+1) - B(N,X) = N * X**(N-1)
//      B(N,X) = (-1)**N * B(N,1-X)
//
//    A formula for the Bernoulli polynomial in terms of the Bernoulli
//    numbers is:
//
//      B(N,X) = sum ( 0 <= K <= N ) B(K) * C(N,K) * X**(N-K)
//
//    The first few polynomials include:
//
//      B(0,X) = 1
//      B(1,X) = X    - 1/2
//      B(2,X) = X**2 -   X      +  1/6
//      B(3,X) = X**3 - 3/2*X**2 +  1/2*X
//      B(4,X) = X**4 - 2*X**3   +      X**2 - 1/30
//      B(5,X) = X**5 - 5/2*X**4 +  5/3*X**3 - 1/6*X
//      B(6,X) = X**6 - 3*X**5   +  5/2*X**4 - 1/2*X**2 + 1/42
//      B(7,X) = X**7 - 7/2*X**6 +  7/2*X**5 - 7/6*X**3 + 1/6*X
//      B(8,X) = X**8 - 4*X**7   + 14/3*X**6 - 7/3*X**4 + 2/3*X**2 - 1/30
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    31 January 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the Bernoulli polynomial to
//    be evaluated.  N must be 0 or greater.
//
//    Input, double X, the value of X at which the polynomial is to
//    be evaluated.
//
//    Output, double BERNOULLI_POLY, the value of B(N,X).
//
{
  int i;
  int *iwork;
  bool next;
  double value;
  double *work;

  work = new double[n+1];
  bernoulli_number ( n, work );

  iwork = new int[n+1];
  next = false;
  comb_row ( next, n, iwork );

  value = 1.0;
  for ( i = 1; i <= n; i++ )
  {
    value = value * x + work[i] * ( double ) iwork[i];
  }

  delete [] iwork;
  delete [] work;

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

double bernoulli_poly2 ( int n, double x )

//****************************************************************************80
//
//  Purpose:
//
//    BERNOULLI_POLY2 evaluates the N-th Bernoulli polynomial at X.
//
//  Discussion:
//
//    Thanks to Bart Vandewoestyne for pointing out an error in the previous
//    documentation, 31 January 2008.
//
//    Special values of the Bernoulli polynomial include:
//
//      B(N,0) = B(N,1) = B(N), the N-th Bernoulli number.
//
//      B'(N,X) = N * B(N-1,X)
//
//      B(N,X+1) - B(N,X) = N * X**(N-1)
//      B(N,X) = (-1)**N * B(N,1-X)
//
//    A formula for the Bernoulli polynomial in terms of the Bernoulli
//    numbers is:
//
//      B(N,X) = sum ( 0 <= K <= N ) B(K) * C(N,K) * X**(N-K)
//
//    The first few polynomials include:
//
//      B(0,X) = 1
//      B(1,X) = X    - 1/2
//      B(2,X) = X**2 -   X      +  1/6
//      B(3,X) = X**3 - 3/2*X**2 +  1/2*X
//      B(4,X) = X**4 - 2*X**3   +      X**2 - 1/30
//      B(5,X) = X**5 - 5/2*X**4 +  5/3*X**3 - 1/6*X
//      B(6,X) = X**6 - 3*X**5   +  5/2*X**4 - 1/2*X**2 + 1/42
//      B(7,X) = X**7 - 7/2*X**6 +  7/2*X**5 - 7/6*X**3 + 1/6*X
//      B(8,X) = X**8 - 4*X**7   + 14/3*X**6 - 7/3*X**4 + 2/3*X**2 - 1/30
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    31 January 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the Bernoulli polynomial to
//    be evaluated.  N must be 0 or greater.
//
//    Input, double X, the value at which the polynomial is to
//    be evaluated.
//
//    Output, double BERNOULLI_POLY2, the value of B(N,X).
//
{
  double fact;
  int i;
  double value;

  fact = 1.0;

  value = bernoulli_number3 ( 0 );

  for ( i = 1; i <= n; i++ )
  {
    fact = fact * ( double ) ( n + 1 - i ) / ( double ) i;
    value = value * x + fact * bernoulli_number3 ( i );
  }

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

void bernoulli_number_values ( int *n_data, int *n, double *c )

//****************************************************************************80
//
//  Purpose:
//
//    BERNOULLI_NUMBER_VALUES returns some values of the Bernoulli numbers.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and N_DATA
//    is set to 1.  On each subsequent call, the input value of N_DATA is
//    incremented and that test data item is returned, if available.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, int *N, the order of the Bernoulli number.
//
//    Output, double *C, the value of the Bernoulli number.
//
{
# define N_MAX 10

  double c_vec[N_MAX] = {
    1.0000000000, -0.5000000000,  0.1666666667, 
    0.0000000000, -0.0333333333, -0.02380952381, 
   -0.0333333333,  0.0757575757, -529.1242424, 
    601580873.9 };

  int n_vec[N_MAX] = {
     0,  1,  2,  3,  4, 6,  8, 10, 20, 30 };

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *n = 0;
    *c = 0;
  }
  else
  {
    *n = n_vec[*n_data];
    *c = c_vec[*n_data];
    *n_data = *n_data + 1;
  }

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

void bernstein_poly ( int n, double x, double bern[] )

//*****************************************************************************
//
//  Purpose:
//
//    BERNSTEIN_POLY evaluates the Bernstein polynomials at a point X.
//
//  Discussion:
//
//    The Bernstein polynomials are assumed to be based on [0,1].
//
//  Formula:
//
//    B(N,I)(X) = [N!/(I!*(N-I)!)] * (1-X)**(N-I) * X**I
//
//  First values:
//
//    B(0,0)(X) = 1
//
//    B(1,0)(X) =      1-X
//    B(1,1)(X) =                X
//
//    B(2,0)(X) =     (1-X)**2
//    B(2,1)(X) = 2 * (1-X)    * X
//    B(2,2)(X) =                X**2
//
//    B(3,0)(X) =     (1-X)**3
//    B(3,1)(X) = 3 * (1-X)**2 * X
//    B(3,2)(X) = 3 * (1-X)    * X**2
//    B(3,3)(X) =                X**3
//
//    B(4,0)(X) =     (1-X)**4
//    B(4,1)(X) = 4 * (1-X)**3 * X
//    B(4,2)(X) = 6 * (1-X)**2 * X**2
//    B(4,3)(X) = 4 * (1-X)    * X**3
//    B(4,4)(X) =                X**4
//
//  Special values:
//
//    B(N,I)(X) has a unique maximum value at X = I/N.
//
//    B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1.
//
//    B(N,I)(1/2) = C(N,K) / 2**N
//
//    For a fixed X and N, the polynomials add up to 1:
//
//      Sum ( 0 <= I <= N ) B(N,I)(X) = 1
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Parameters:
//
//    Input, int N, the degree of the Bernstein polynomials to be
//    used.  For any N, there is a set of N+1 Bernstein polynomials,
//    each of degree N, which form a basis for polynomials on [0,1].
//
//    Input, double X, the point at which the polynomials are to be evaluated.
//
//    Output, double BERN[N+1], the values of the Bernstein polynomials 
//    of orders 0 through N at X.
//
{
  int i;
  int j;

  if ( n == 0 )
  {
    bern[0] = 1.0;
  }
  else if ( 0 < n )
  {
    bern[0] = 1.0 - x;
    bern[1] = x;

    for ( i = 2; i <= n; i++ )
    {
      bern[i] = x * bern[i-1];
      for ( j = i-1; 1 <= j; j-- )
      {
        bern[j] = x * bern[j-1] + ( 1.0 - x ) * bern[j];
      }
      bern[0] = ( 1.0 - x ) * bern[0];
    }

  }

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

void bernstein_poly_values ( int *n_data, int *n, int *k, double *x, double *b )

//****************************************************************************80
//
//  Purpose:
//
//    BERNSTEIN_POLY_VALUES returns some values of the Bernstein polynomials.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and
//    N_DATA is set to the index of the test data.  On each subsequent
//    call, N_DATA is incremented and that test data is returned.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, int *N, the degree of the polynomial.
//
//    Output, int *K, the index of the polynomial.
//
//    Output, double *X, the argument of the polynomial.
//
//    Output, double *B, the value of the polynomial B(N,K)(X).
//
{
# define N_MAX 15
 
  double b_vec[N_MAX] = {
    1.0, 
    0.75,       0.25, 
    0.5625,     0.3750,   0.0625, 
    0.421875,   0.421875, 0.140625,  0.015625, 
    0.31640625, 0.421875, 0.2109375, 0.046875, 0.00390625 };
  int k_vec[N_MAX] = {
    0, 
    0, 1, 
    0, 1, 2, 
    0, 1, 2, 3, 
    0, 1, 2, 3, 4 };
  int n_vec[N_MAX] = {
    0, 
    1, 1, 
    2, 2, 2, 
    3, 3, 3, 3, 
    4, 4, 4, 4, 4 };
  double x_vec[N_MAX] = {
    0.25, 
    0.25, 0.25, 
    0.25, 0.25, 0.25, 
    0.25, 0.25, 0.25, 0.25, 
    0.25, 0.25, 0.25, 0.25, 0.25 };

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *n = 0;
    *k = 0;
    *x = 0.0;
    *b = 0.0;
  }
  else
  {
    *n = n_vec[*n_data];
    *k = k_vec[*n_data];
    *x = x_vec[*n_data];
    *b = b_vec[*n_data];
    *n_data = *n_data + 1;
  }

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

double beta ( double x, double y )

//****************************************************************************80
//
//  Purpose:
//
//    BETA returns the value of the Beta function.
//
//  Formula:
//
//    BETA(X,Y) = ( GAMMA(X) * GAMMA(Y) ) / GAMMA(X+Y)
//
//  Properties:
//
//    BETA(X,Y) = BETA(Y,X).
//    BETA(X,Y) = Integral ( 0 <= T <= 1 ) T**(X-1) (1-T)**(Y-1) dT.
//    BETA(X,Y) = GAMMA(X) * GAMMA(Y) / GAMMA(X+Y)
//
//  Restrictions:
//
//    Both X and Y must be greater than 0.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the two parameters that define the Beta function.
//    X and Y must be greater than 0.
//
//    Output, double BETA, the value of the Beta function.
//
{
  double value;

  if ( x <= 0.0 || y <= 0.0 )
  {
    cout << "\n";
    cout << "BETA - Fatal error!\n";
    cout << "  Both X and Y must be greater than 0.\n";
    exit ( 1 );
  }

  value = exp ( gamma_log ( x ) + gamma_log ( y ) - gamma_log ( x + y ) );

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

void beta_values ( int *n_data, double *x, double *y, double *fxy )

//****************************************************************************80
//
//  Purpose:
//
//    BETA_VALUES returns some values of the Beta function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and
//    N_DATA is set to the index of the test data.  On each subsequent
//    call, N_DATA is incremented and that test data is returned.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, double *X, *Y, the arguments of the function.
//
//    Output, double *FXY, the value of the function.
//
{
# define N_MAX 17

  double b_vec[N_MAX] = {
    5.000000, 2.500000, 1.666667, 1.250000, 
    5.000000, 2.500000, 1.000000, 1.666667E-01, 
    0.333333E-01, 7.142857E-03, 1.587302E-03, 0.238095E-01, 
    5.952381E-03, 1.984127E-03, 7.936508E-04, 3.607504E-04, 
    8.325008E-05 };

  double x_vec[N_MAX] = {
    0.2, 0.4, 0.6, 0.8, 
    1.0, 1.0, 1.0, 2.0, 
    3.0, 4.0, 5.0, 6.0, 
    6.0, 6.0, 6.0, 6.0, 
    7.0 };

  double y_vec[N_MAX] = {
    1.0, 1.0, 1.0, 1.0, 
    0.2, 0.4, 1.0, 2.0, 
    3.0, 4.0, 5.0, 2.0, 
    3.0, 4.0, 5.0, 6.0, 
    7.0 };

  if ( n_data < 0 )
  {
    *n_data = 0;
  }
  
  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *y = 0.0;
    *fxy = 0.0;
  }
  else
  {
    *x = x_vec[*n_data];
    *y = y_vec[*n_data];
    *fxy = b_vec[*n_data];
    *n_data = *n_data + 1;
  }

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

void bpab ( int n, double x, double a, double b, double bern[] )

//****************************************************************************80
//
//  Purpose:
//
//    BPAB evaluates at X the Bernstein polynomials based in [A,B].
//
//  Formula:
//
//    BERN(N,I)(X) = [N!/(I!*(N-I)!)] * (B-X)**(N-I) * (X-A)**I / (B-A)**N
//
//  First values:
//
//    B(0,0)(X) =   1
//
//    B(1,0)(X) = (      B-X                ) / (B-A)
//    B(1,1)(X) = (                 X-A     ) / (B-A)
//
//    B(2,0)(X) = (     (B-X)**2            ) / (B-A)**2
//    B(2,1)(X) = ( 2 * (B-X)    * (X-A)    ) / (B-A)**2
//    B(2,2)(X) = (                (X-A)**2 ) / (B-A)**2
//
//    B(3,0)(X) = (     (B-X)**3            ) / (B-A)**3
//    B(3,1)(X) = ( 3 * (B-X)**2 * (X-A)    ) / (B-A)**3
//    B(3,2)(X) = ( 3 * (B-X)    * (X-A)**2 ) / (B-A)**3
//    B(3,3)(X) = (                (X-A)**3 ) / (B-A)**3
//
//    B(4,0)(X) = (     (B-X)**4            ) / (B-A)**4
//    B(4,1)(X) = ( 4 * (B-X)**3 * (X-A)    ) / (B-A)**4
//    B(4,2)(X) = ( 6 * (B-X)**2 * (X-A)**2 ) / (B-A)**4
//    B(4,3)(X) = ( 4 * (B-X)    * (X-A)**3 ) / (B-A)**4
//    B(4,4)(X) = (                (X-A)**4 ) / (B-A)**4
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the degree of the Bernstein polynomials to be used.
//    For any N, there is a set of N+1 Bernstein polynomials, each of
//    degree N, which form a basis for polynomials on [A,B].
//
//    Input, double X, the point at which the polynomials are to be evaluated.
//
//    Input, double A, B, the endpoints of the interval on which the
//    polynomials are to be based.  A and B should not be equal.
//
//    Output, double BERN[N+1], the values of the N+1 Bernstein polynomials at X.
//
{
  int i;
  int j;

  if ( b == a )
  {
    cout << "\n";
    cout << "BPAB - Fatal error!\n";
    cout << "  A = B = " << a << "\n";
    exit ( 1 );
  }

  if ( n == 0 )
  {
    bern[0] = 1.0;
  }
  else if ( 0 < n );
  {
    bern[0] = ( b - x ) / ( b - a );
    bern[1] = ( x - a ) / ( b - a );

    for ( i = 2; i <= n; i++ )
    {
      bern[i] = ( x - a ) * bern[i-1] / ( b - a );
      for ( j = i-1; 1 <= j; j-- )
      {
        bern[j] = ( ( b - x     ) * bern[j] 
                  + (     x - a ) * bern[j-1] ) 
                  / ( b     - a );
      }
      bern[0] = ( b - x ) * bern[0] / ( b - a );
    }

  }

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

void cardan ( int n, double x, double s, double cx[] )

//****************************************************************************80
//
//  Purpose:
//
//    CARDAN evaluates the Cardan polynomials.
//
//  First terms:
//
//     N  C(N,S,X)
//
//     0  2
//     1  X
//     2  X**2  -  2 S
//     3  X**3  -  3 S X
//     4  X**4  -  4 S X**2 +  2 S**2
//     5  X**5  -  5 S X**3 +  5 S**2 X
//     6  X**6  -  6 S X**4 +  9 S**2 X**2 -  2 S**3
//     7  X**7  -  7 S X**5 + 14 S**2 X**3 -  7 S**3 X
//     8  X**8  -  8 S X**6 + 20 S**2 X**4 - 16 S**3 X**2 +  2 S**4
//     9  X**9  -  9 S X**7 + 27 S**2 X**5 - 30 S**3 X**3 +  9 S**4 X
//    10  X**10 - 10 S X**8 + 35 S**2 X**6 - 50 S**3 X**4 + 25 S**4 X**2 -  2 S**5
//    11  X**11 - 11 S X**9 + 44 S**2 X**7 - 77 S**3 X**5 + 55 S**4 X**3 - 11 S**5 X
//
//  Recursion:
//
//    Writing the N-th polynomial in terms of its coefficients:
//
//      C(N,S,X) = sum ( 0 <= I <= N ) D(N,I) * S**(N-I)/2 * X**I
//
//    then
//
//    D(0,0) = 1
//
//    D(1,1) = 1
//    D(1,0) = 0
//
//    D(N,N) = 1
//    D(N,K) = D(N-1,K-1) - D(N-2,K)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Thomas Osler,
//    Cardan Polynomials and the Reduction of Radicals,
//    Mathematics Magazine,
//    Volume 74, Number 1, February 2001, pages 26-32.
//
//  Parameters:
//
//    Input, int N, the highest polynomial to compute.
//
//    Input, double X, the point at which the polynomials are to be computed.
//
//    Input, double S, the value of the parameter, which must be positive.
//
//    Output, double CX[N+1], the values of the Cardan polynomials at X.
//
{
  double fact;
  int i;
  double s2;
  double x2;

  s2 = sqrt ( s );
  x2 = 0.5 * x / s2;

  cheby_t_polynomial ( n, x2, cx );

  fact = 1.0;

  for ( i = 0; i <= n; i++ )
  {
    cx[i] = 2.0 * fact * cx[i];
    fact = fact * s2;
  }

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

void cardan_poly_coef ( int n, double s, double c[] )

//****************************************************************************80
//
//  Purpose:
//
//    CARDAN_POLY_COEF computes the coefficients of the N-th Cardan polynomial.
//
//  First terms:
//
//    2
//    0       1
//   -2 S     0       1
//    0      -3 S     0       1
//    2 S**2  0      -4 S     0       1
//    0       5 S**2  0      -5 S     0       1
//   -2 S**3  0       9 S**2  0      -6 S     0       1
//    0       7 S**3  0      14 S**2  0      -7 S     0       1
//    2 S**4  0     -16 S**3  0      20 S**2  0      -8 S     0        1
//    0       9 S**4  0     -30 S**3  0      27 S**2  0      -9 S      0     1
//   -2 S**5  0      25 S**4  0     -50 S**3  0      35 S**2  0      -10 S   0   1
//    0     -11 S**5  0      55 S**4  0     -77 S**3  0     +44 S**2   0   -11 S 0 1
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Thomas Osler,
//    Cardan Polynomials and the Reduction of Radicals,
//    Mathematics Magazine,
//    Volume 74, Number 1, February 2001, pages 26-32.
//
//  Parameters:
//
//    Input, int N, the order of the polynomial
//
//    Input, double S, the value of the parameter, which must be positive.
//
//    Output, double C[N+1], the coefficients.  C(0) is the constant term,
//    and C(N) is the coefficient of X**N.
//
{
  double *cm1;
  double *cm2;
  int i;
  int j;

  if ( n < 0 )
  {
    return;
  }

  c[0] = 2.0;
  for ( i = 1; i <= n; i++ )
  {
    c[i] = 0.0;
  }

  if ( n == 0 )
  {
    return;
  }

  cm1 = new double[n+1];
  cm2 = new double[n+1];

  for ( i = 0; i <= n; i++ )
  {
    cm1[i] = c[i];
  }

  c[0] = 0.0;
  c[1] = 1.0;
  for ( i = 2; i <= n; i++ )
  {
    c[i] = 0.0;
  }

  for ( i = 2; i <= n; i++ )
  {

    for ( j = 0; j <= i-2; j++ )
    {
      cm2[j] = cm1[j];
    }

    for ( j = 0; j <= i-1; j++ )
    {
      cm1[j] = c[j];
    }

    c[0] = 0.0;
    for ( j = 1; j <= i; j++ )
    {
      c[j] = cm1[j-1];
    }

    for ( j = 0; j <= i-2; j++ )
    {
      c[j] = c[j] - s * cm2[j];
    }
  }

  delete [] cm1;
  delete [] cm2;

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

void catalan ( int n, int c[] )

//****************************************************************************80
//
//  Purpose:
//
//    CATALAN computes the Catalan numbers, from C(0) to C(N).
//
//  First values:
//
//     C(0)     1
//     C(1)     1
//     C(2)     2
//     C(3)     5
//     C(4)    14
//     C(5)    42
//     C(6)   132
//     C(7)   429
//     C(8)  1430
//     C(9)  4862
//    C(10) 16796
//
//  Formula:
//
//    C(N) = (2*N)! / ( (N+1) * (N!) * (N!) )
//         = 1 / (N+1) * COMB ( 2N, N )
//         = 1 / (2N+1) * COMB ( 2N+1, N+1).
//
//  Recursion:
//
//    C(N) = 2 * (2*N-1) * C(N-1) / (N+1)
//    C(N) = sum ( 1 <= I <= N-1 ) C(I) * C(N-I)
//
//  Discussion:
//
//    The Catalan number C(N) counts:
//
//    1) the number of binary trees on N vertices;
//    2) the number of ordered trees on N+1 vertices;
//    3) the number of full binary trees on 2N+1 vertices;
//    4) the number of well formed sequences of 2N parentheses;
//    5) the number of ways 2N ballots can be counted, in order,
//       with N positive and N negative, so that the running sum
//       is never negative;
//    6) the number of standard tableaus in a 2 by N rectangular Ferrers diagram;
//    7) the number of monotone functions from [1..N} to [1..N} which
//       satisfy f(i) <= i for all i;
//    8) the number of ways to triangulate a polygon with N+2 vertices.
//
//  Example:
//
//    N = 3
//
//    ()()()
//    ()(())
//    (()())
//    (())()
//    ((()))
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Dennis Stanton and Dennis White,
//    Constructive Combinatorics,
//    Springer Verlag, New York, 1986.
//
//  Parameters:
//
//    Input, int N, the number of Catalan numbers desired.
//
//    Output, int C[N+1], the Catalan numbers from C(0) to C(N).
//
{
  int i;

  if ( n < 0 )
  {
    return;
  }

  c[0] = 1;
//
//  The extra parentheses ensure that the integer division is
//  done AFTER the integer multiplication.
//
  for ( i = 1; i <= n; i++ )
  {
    c[i] = ( c[i-1] * 2 * ( 2 * i - 1 ) ) / ( i + 1 );
  }

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

void catalan_row_next ( bool next, int n, int irow[] )

//****************************************************************************80*
//
//  Purpose:
//
//    CATALAN_ROW_NEXT computes row N of Catalan's triangle.
//
//  Example:
//
//    I\J 0   1   2   3   4   5   6
//
//    0   1
//    1   1   1
//    2   1   2   2
//    3   1   3   5   5
//    4   1   4   9  14  14
//    5   1   5  14  28  42  42
//    6   1   6  20  48  90 132 132
//
//  Recursion:
//
//    C(0,0) = 1
//    C(I,0) = 1
//    C(I,J) = 0 for I < J
//    C(I,J) = C(I,J-1) + C(I-1,J)
//    C(I,I) is the I-th Catalan number.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, bool NEXT, indicates whether this is a call for
//    the 'next' row of the triangle.
//    NEXT = FALSE, this is a startup call.  Row N is desired, but
//    presumably this is a first call, or row N-1 was not computed
//    on the previous call.
//    NEXT = TRUE, this is not the first call, and row N-1 was computed
//    on the previous call.  In this case, much work can be saved
//    by using the information from the previous values of IROW
//    to build the next values.
//
//    Input, int N, the index of the row of the triangle desired.
//
//    Input/output, int IROW(0:N), the row of coefficients.
//    If NEXT = FALSE, then IROW is not required to be set on input.
//    If NEXT = TRUE, then IROW must be set on input to the value of
//    row N-1.
//
{
  int i;
  int j;

  if ( n < 0 )
  {
    return;
  }

  if ( !next )
  {
    irow[0] = 1;
    for ( i = 1; i <= n; i++ )
    {
      irow[i] = 0;
    }

    for ( i = 1; i <= n; i++ )
    {
      irow[0] = 1;

      for ( j = 1; j <= i-1; j++ )
      {
        irow[j] = irow[j] + irow[j-1];
      }

      irow[i] = irow[i-1];

    }
  }
  else
  {
    irow[0] = 1;

    for ( j = 1; j <= n-1; j++ )
    {
      irow[j] = irow[j] + irow[j-1];
    }

    if ( 1 <= n )
   {
      irow[n] = irow[n-1];
    }

  }

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

void catalan_values ( int *n_data, int *n, int *c )

//****************************************************************************80
//
//  Purpose:
//
//    CATALAN_VALUES returns some values of the Catalan numbers.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and N_DATA
//    is set to 1.  On each subsequent call, the input value of N_DATA is
//    incremented and that test data item is returned, if available.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, int *N, the order of the Catalan number.
//
//    Output, int *C, the value of the Catalan number.
//
{
# define N_MAX 11

  int c_vec[N_MAX] = { 1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796 };
  int n_vec[N_MAX] = { 0,  1,  2,  3,  4, 5,  6,  7,  8,  9,  10 };

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *n = 0;
    *c = 0;
  }
  else
  {
    *n = n_vec[*n_data];
    *c = c_vec[*n_data];
    *n_data = *n_data + 1;
  }

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

double charlier ( int n, double a, double x )

//****************************************************************************80
//
//  Purpose:
//
//    CHARLIER evaluates a Poisson-Charlier polynomial at a point.
//
//  Discussion:
//
//    p(n,a)(x) = a^(n/2) * (n!)^(-1/2) * sum ( 0 <= j <= n ) (-1)*(n-j)
//      * C(n,j) * j! * a^(-j) * C(x,j)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    17 February 2005
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    J Simoes Pereira,
//    Algorithm 234: Poisson-Charliers Polynomials,
//    Communications of the ACM,
//    Volume 7, Number 7, page 420, July 1964.
//
//    G Szego,
//    Orthogonal Polynomials,
//    American Mathematical Society, 1975.
//
//    Eric Weisstein,
//    CRC Concise Encyclopedia of Mathematics,
//    CRC Press, 2002,
//    Second edition,
//    ISBN: 1584883472,
//    LC: QA5.W45.
//
//  Parameters:
//
//    Input, int N, the order of the polynomial.  N must be at least 0.
//
//    Input, double A, the parameter.  A must be greater than
//    or equal to 0.
//
//    Input, double X, the point at which the polynomial is
//    to be evaluated.
//
//    Output, double CHARLIER, the value of the polynomial.
//
{
  double c;
  int j;
  double s;
  double u;
  double value;

  if ( a < 0.0 )
  {
    value = 0.0;
  }
  else if ( a == 0.0 )
  {
    value = 0.0;
  }
  else if ( n < 0 )
  {
    value = 0.0;
  }
  else if ( n == 0 )
  {
    value = 1.0;
  }
  else
  {
    if ( ( n % 2 ) == 0 )
    {
      u = 1.0;
    }
    else
    {
      u = -1.0;
    }

    s = u;

    for ( j = 0; j <= n - 1; j++ )
    {
      u = - u * ( double ) ( n - j ) * ( x - ( double ) ( j ) )
        / ( a * ( double ) ( j + 1 ) );
      s = s + u;
    }

    c = 1.0;
    for ( j = 1; j <= n; j++ )
    {
      c = c * ( double ) ( j );
    }

    value = s * sqrt ( pow ( a, n ) / c );

  }

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

void cheby_t_polynomial ( int n, double x, double cx[] )

//****************************************************************************80
//
//  Purpose:
//
//    CHEBY_T_POLYNOMIAL evaluates Chebyshev polynomials T(N)(X) of the first kind.
//
//  Discussion:
//
//    Chebyshev polynomials are useful as a basis for representing the
//    approximation of functions since they are well conditioned, in the sense
//    that in the interval [-1,1] they each have maximum absolute value 1.
//    Hence an error in the value of a coefficient of the approximation, of
//    size epsilon, is exactly reflected in an error of size epsilon between
//    the computed approximation and the theoretical approximation.
//
//    Typical usage is as follows, where we assume for the moment
//    that the interval of approximation is [-1,1].  The value
//    of N is chosen, the highest polynomial to be used in the
//    approximation.  Then the function to be approximated is
//    evaluated at the N+1 points XJ which are the zeroes of the N+1-th
//    Chebyshev polynomial.  Let these values be denoted by F(XJ).
//
//    The coefficients of the approximation are now defined by
//
//      C(I) = 2/(N+1) * sum ( 1 <= J <= N+1 ) F(XJ) T(I)(XJ)
//
//    except that C(0) is given a value which is half that assigned
//    to it by the above formula,
//
//    and the representation is
//
//    F(X) approximated by sum ( 0 <= J <= N ) C(J) T(J)(X)
//
//    Now note that, again because of the fact that the Chebyshev polynomials
//    have maximum absolute value 1, if the higher order terms of the
//    coefficients C are small, then we have the option of truncating
//    the approximation by dropping these terms, and we will have an
//    exact value for maximum perturbation to the approximation that
//    this will cause.
//
//    It should be noted that typically the error in approximation
//    is dominated by the first neglected basis function (some multiple of
//    T(N+1)(X) in the example above).  If this term were the exact error,
//    then we would have found the minimax polynomial, the approximating
//    polynomial of smallest maximum deviation from the original function.
//    The minimax polynomial is hard to compute, and another important
//    feature of the Chebyshev approximation is that it tends to behave
//    like the minimax polynomial while being easy to compute.
//
//    To evaluate a sum like
//
//      sum ( 0 <= J <= N ) C(J) T(J)(X),
//
//    Clenshaw's recurrence formula is recommended instead of computing the
//    polynomial values, forming the products and summing.
//
//    Assuming that the coefficients C(J) have been computed
//    for J = 0 to N, then the coefficients of the representation of the
//    indefinite integral of the function may be computed by
//
//      B(I) = ( C(I-1) - C(I+1))/2*(I-1) for I=1 to N+1,
//
//    with
//
//      C(N+1)=0
//      B(0) arbitrary.
//
//    Also, the coefficients of the representation of the derivative of the
//    function may be computed by:
//
//      D(I) = D(I+2)+2*I*C(I) for I=N-1, N-2, ..., 0,
//
//    with
//
//      D(N+1) = D(N)=0.
//
//    Some of the above may have to adjusted because of the irregularity of C(0).
//
//  Differential equation:
//
//    (1-X*X) Y'' - X Y' + N N Y = 0
//
//  Formula:
//
//    T(N)(X) = COS(N*ARCCOS(X))
//
//  First terms:
//
//    T(0)(X) =  1
//    T(1)(X) =  1 X
//    T(2)(X) =  2 X**2 -   1
//    T(3)(X) =  4 X**3 -   3 X
//    T(4)(X) =  8 X**4 -   8 X**2 +  1
//    T(5)(X) = 16 X**5 -  20 X**3 +  5 X
//    T(6)(X) = 32 X**6 -  48 X**4 + 18 X**2 - 1
//    T(7)(X) = 64 X**7 - 112 X**5 + 56 X**3 - 7 X
//
//  Inequality:
//
//    abs ( T(N)(X) ) <= 1 for -1 <= X <= 1
//
//  Orthogonality:
//
//    For integration over [-1,1] with weight
//
//      W(X) = 1 / sqrt(1-X*X),
//
//    if we write the inner product of T(I)(X) and T(J)(X) as
//
//      < T(I)(X), T(J)(X) > = integral ( -1 <= X <= 1 ) W(X) T(I)(X) T(J)(X) dX
//
//    then the result is:
//
//      0 if I /= J
//      PI/2 if I == J /= 0
//      PI if I == J == 0
//
//    A discrete orthogonality relation is also satisfied at each of
//    the N zeroes of T(N)(X):  sum ( 1 <= K <= N ) T(I)(X) * T(J)(X)
//                              = 0 if I /= J
//                              = N/2 if I == J /= 0
//                              = N if I == J == 0
//
//  Recursion:
//
//    T(0)(X) = 1,
//    T(1)(X) = X,
//    T(N)(X) = 2 * X * T(N-1)(X) - T(N-2)(X)
//
//    T'(N)(X) = N * ( -X * T(N)(X) + T(N-1)(X) ) / ( 1 - X**2 )
//
//  Special values:
//
//    T(N)(1) = 1
//    T(N)(-1) = (-1)**N
//    T(2N)(0) = (-1)**N
//    T(2N+1)(0) = 0
//    T(N)(X) = (-1)**N * T(N)(-X)
//
//  Zeroes:
//
//    M-th zero of T(N)(X) is cos((2*M-1)*PI/(2*N)), M = 1 to N
//
//  Extrema:
//
//    M-th extremum of T(N)(X) is cos(PI*M/N), M = 0 to N
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the highest polynomial to compute.
//
//    Input, double X, the point at which the polynomials are to be computed.
//
//    Output, double CX[N+1], the values of the N+1 Chebyshev polynomials.
//
{
  int i;

  if ( n < 0 )
  {
    return;
  }

  cx[0] = 1.0;

  if ( n < 1 )
  {
    return;
  }

  cx[1] = x;

  for ( i = 2; i <= n; i++ )
  {
    cx[i] = 2.0 * x * cx[i-1] - cx[i-2];
  }

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

void cheby_t_polynomial_coef ( int n, double c[] )

//****************************************************************************80
//
//  Purpose:
//
//    CHEBY_T_POLYNOMIAL_COEF evaluates the Chebyshev T polynomial coefficients.
//
//  First terms:
//
//    N/K     0     1      2      3       4     5      6    7      8    9   10
//
//     0      1
//     1      0     1
//     2     -1     0      2
//     3      0    -3      0      4
//     4      1     0     -8      0       8
//     5      0     5      0    -20       0    16
//     6     -1     0     18      0     -48     0     32
//     7      0    -7      0     56       0  -112      0    64
//
//  Recursion:
//
//    T(0)(X) = 1,
//    T(1)(X) = X,
//    T(N)(X) = 2 * X * T(N-1)(X) - T(N-2)(X)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    16 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input, int N, the highest order polynomial to compute.
//    Note that polynomials 0 through N will be computed.
//
//    Output, double C[(N+1)*(N+1)], the coefficients of the Chebyshev T
//    polynomials.
//
{
  int i;
  int j;

  if ( n < 0 )
  {
    return;
  }

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

  c[0+0*(n+1)] = 1.0;

  if ( n == 0 )
  {
    return;
  }

  c[1+1*(n+1)] = 1.0;

  for ( i = 2; i <= n; i++ )
  {
    c[i+0*(n+1)]     =                    - c[i-2+0*(n+1)];
    for ( j = 1; j <= i-2; j++ )
    {
      c[i+j*(n+1)] = 2.0 * c[i-1,j-1] - c[i-2+(j-1)*(n+1)];
    }
    c[i+(i-1)*(n+1)] = 2.0 * c[i-1+(i-2)*(n+1)];
    c[i+ i   *(n+1)] = 2.0 * c[i-1+(i-1)*(n+1)];
  }

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

void cheby_t_polynomial_values ( int *n_data, int *n, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    CHEBY_T_POLYNOMIAL_VALUES returns values of the Chebyshev T polynomial.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    04 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and N_DATA
//    is set to 1.  On each subsequent call, the input value of N_DATA is
//    incremented and that test data item is returned, if available.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, int *N, the order of the function.
//
//    Output, double *X, the point where the function is evaluated.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 13

  double fx_vec[N_MAX] = {
     1.0000000000,  0.8000000000,  0.2800000000, 
    -0.3520000000, -0.8432000000, -0.9971200000, 
    -0.7521920000, -0.2063872000,  0.4219724800, 
     0.8815431680,  0.9884965888,  0.7000513741, 
     0.1315856097 };
  int n_vec[N_MAX] = {
     0,  1,  2, 
     3,  4,  5, 
     6,  7,  8, 
     9, 10, 11, 
    12 };
  double x_vec[N_MAX] = {
    0.8,  0.8,  0.8, 
    0.8,  0.8,  0.8, 
    0.8,  0.8,  0.8, 
    0.8,  0.8,  0.8, 
    0.8 };

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *n = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *n = n_vec[*n_data];
    *x = x_vec[*n_data];
    *fx = fx_vec[*n_data];
    *n_data = *n_data + 1;
  }

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

void cheby_u_polynomial ( int n, double x, double cx[] )

//****************************************************************************80
//
//  Purpose:
//
//    CHEBY_U_POLYNOMIAL evaluates the Chebyshev polynomials of the second kind.
//
//  Differential equation:
//
//    (1-X*X) Y'' - 3 X Y' + N (N+2) Y = 0
//
//  First terms:
//
//    U(0)(X) =   1
//    U(1)(X) =   2 X
//    U(2)(X) =   4 X**2 -   1
//    U(3)(X) =   8 X**3 -   4 X
//    U(4)(X) =  16 X**4 -  12 X**2 +  1
//    U(5)(X) =  32 X**5 -  32 X**3 +  6 X
//    U(6)(X) =  64 X**6 -  80 X**4 + 24 X**2 - 1
//    U(7)(X) = 128 X**7 - 192 X**5 + 80 X**3 - 8X
//
//  Recursion:
//
//    U(0)(X) = 1,
//    U(1)(X) = 2 * X,
//    U(N)(X) = 2 * X * U(N-1)(X) - U(N-2)(X)
//
//  Norm:
//
//    Integral ( -1 <= X <= 1 ) ( 1 - X**2 ) * U(N)(X)**2 dX = PI/2
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the highest polynomial to compute.
//
//    Input, double X, the point at which the polynomials are to be computed.
//
//    Output, double CX[N+1], the values of the N+1 Chebyshev polynomials.
//
{
  int i;

  if ( n < 0 )
  {
    return;
  }

  cx[0] = 1.0;

  if ( n < 1 )
  {
    return;
  }

  cx[1] = 2.0 * x;

  for ( i = 2; i <= n; i++ )
  {
    cx[i] = 2.0 * x * cx[i-1] - cx[i-2];
  }

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

void cheby_u_polynomial_coef ( int n, double c[] )

//****************************************************************************80
//
//  Purpose:
//
//    CHEBY_U_POLYNOMIAL_COEF evaluates the Chebyshev U polynomial coefficients.
//
//  First terms:
//
//    N/K     0     1      2      3       4     5      6    7      8    9   10
//
//     0      1
//     1      0     2
//     2     -1     0      4
//     3      0    -4      0      8
//     4      1     0    -12      0      16
//     5      0     6      0    -32       0    32
//     6     -1     0     24      0     -80     0     64
//     7      0    -8      0     80       0  -192      0   128
//
//  Recursion:
//
//    U(0)(X) = 1,
//    U(1)(X) = 2*X,
//    U(N)(X) = 2 * X * U(N-1)(X) - U(N-2)(X)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    16 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input, int N, the highest order polynomial to compute.
//    Note that polynomials 0 through N will be computed.
//
//    Output, double C[(N+1)*((N+1)], the coefficients of the Chebyshev U
//    polynomials.
//
{
  int i;
  int j;

  if ( n < 0 )
  {
    return;
  }

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

  c[0+0*(n+1)] = 1.0;

  if ( n == 0 )
  {
    return;
  }

  c[1+1*(n+1)] = 2.0;

  for ( i = 2; i <= n; i++ )
  {
    c[i+0*(n+1)]     =                    - c[i-2+0*(n+1)];
    for ( j = 1; j <= i-2; j++ )
    {
      c[i+j*(n+1)] = 2.0 * c[i-1+(j-1)*(n+1)] - c[i-2+j*(n+1)];
    }

    c[i+(i-1)*(n+1)] = 2.0 * c[i-1+(i-2)*(n+1)];
    c[i+ i   *(n+1)] = 2.0 * c[i-1+(i-1)*(n+1)];
  }

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

void cheby_u_polynomial_values ( int *n_data, int *n, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    CHEBY_U_POLYNOMIAL_VALUES returns values of the Chebyshev U polynomial.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    04 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and N_DATA
//    is set to 1.  On each subsequent call, the input value of N_DATA is
//    incremented and that test data item is returned, if available.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, int *N, the order of the function.
//
//    Output, double *X, the point where the function is evaluated.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 13

  double fx_vec[N_MAX] = {
     1.0000000000,  0.4000000000, -0.8400000000, 
    -0.7360000000,  0.5456000000,  0.9542400000, 
    -0.1639040000, -1.0198016000, -0.2440166400, 
     0.9221949440,  0.6128946176, -0.6770370970, 
    -0.8837094564 };
  int n_vec[N_MAX] = {
     0,  1,  2, 
     3,  4,  5, 
     6,  7,  8, 
     9, 10, 11, 
    12 };
  double x_vec[N_MAX] = {
    0.2,  0.2,  0.2, 
    0.2,  0.2,  0.2, 
    0.2,  0.2,  0.2, 
    0.2,  0.2,  0.2, 
    0.2 };

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *n = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *n = n_vec[*n_data];
    *x = x_vec[*n_data];
    *fx = fx_vec[*n_data];
    *n_data = *n_data + 1;
  }

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

int collatz_count ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    COLLATZ_COUNT counts the number of terms in a Collatz sequence.
//
//  Discussion:
//
//    The rules for generation of the Collatz sequence are recursive.
//    If T is the current entry of the sequence, (T is
//    assumed to be a positive integer), then the next
//    entry, U is determined as follows:
//
//      if T is 1 (or less)
//        terminate the sequence;
//      else if T is even
//        U = T/2.
//      else (if T is odd and not 1)
//        U = 3*T+1;
//
//     N  Sequence                                                Length
//
//     1                                                               1
//     2   1                                                           2
//     3  10,  5, 16,  8,  4,  2,  1                                   8
//     4   2   1                                                       3
//     5  16,  8,  4,  2,  1                                           6
//     6   3, 10,  5, 16,  8,  4,  2,  1                               9
//     7  22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1   17
//     8   4,  2,  1                                                   4
//     9  28, 14,  7, ...                                             20
//    10   5, 16,  8,  4,  2,  1                                       7
//    11  34, 17, 52, 26, 13, 40, 20, 10,  5, 16, 8, 4, 2, 1          15
//    12   6,  3, 10,  5, 16,  8,  4,  2,  1                          10
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    09 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Eric Weisstein,
//    CRC Concise Encyclopedia of Mathematics,
//    CRC Press, 2002,
//    Second edition,
//    ISBN: 1584883472,
//    LC: QA5.W45.
//
//  Parameters:
//
//    Input, int N, the first element of the sequence.
//
//    Output, int COLLATZ_COUNT, the number of elements in
//    the Collatz sequence that begins with N.
//
{
  int count;

  count = 1;

  for ( ; ; )
  {
    if ( n <= 1 )
    {
      break;
    }
    else if ( ( n % 2 ) == 0 )
    {
      n = n / 2;
    }
    else
    {
      n = 3 * n + 1;
    }
    count = count + 1;
  }

  return count;
}
//****************************************************************************80

void collatz_count_values ( int *n_data, int *n, int *count )

//****************************************************************************80
//
//  Purpose:
//
//    COLLATZ_COUNT_VALUES returns some values of the Collatz count function.
//
//  Discussion:
//
//    The rules for generation of the Collatz sequence are recursive.
//    If T is the current entry of the sequence, (T is
//    assumed to be a positive integer), then the next
//    entry, U is determined as follows:
//   
//      if T is 1 (or less)
//        terminate the sequence;
//      else if T is even
//        U = T/2.
//      else (if T is odd and not 1)
//        U = 3*T+1;
//
//    The Collatz count is the length of the Collatz sequence for a given
//    starting value.  By convention, we include the initial value in the
//    count, so the minimum value of the count is 1.
//
//     N  Sequence                                                 Count
//
//     1                                                               1
//     2   1                                                           2
//     3  10,  5, 16,  8,  4,  2,  1                                   8
//     4   2   1                                                       3
//     5  16,  8,  4,  2,  1                                           6
//     6   3, 10,  5, 16,  8,  4,  2,  1                               9
//     7  22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1   17
//     8   4,  2,  1                                                   4
//     9  28, 14,  7, ...                                             20
//    10   5, 16,  8,  4,  2,  1                                       7
//    11  34, 17, 52, 26, 13, 40, 20, 10,  5, 16, 8, 4, 2, 1          15
//    12   6,  3, 10,  5, 16,  8,  4,  2,  1                          10
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    07 March 2006
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Eric Weisstein,
//    CRC Concise Encyclopedia of Mathematics,
//    CRC Press, 2002,
//    Second edition,
//    ISBN: 1584883472,
//    LC: QA5.W45.
//
//  Parameters:
//
//    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, int *N, the initial value of a Collatz sequence.
//
//    Output, int *COUNT, the length of the Collatz sequence starting
//    with N.
//
{
# define N_MAX 20

  int count_vec[N_MAX] = {
      1,   2,   8,   3,   6,   9,   17,   4,  20,   7, 
    112,  25,  26,  27,  17,  28,  111,  18,  83,  29 };
  int n_vec[N_MAX] = {
      1,   2,   3,   4,   5,   6,   7,   8,   9,  10, 
     27,  50, 100, 200, 300, 400, 500, 600, 700, 800 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *count = 0;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *count = count_vec[*n_data-1];
  }

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

void comb_row ( bool next, int n, int row[] )

//****************************************************************************80
//
//  Purpose:
//
//    COMB_ROW computes row N of Pascal's triangle.
//
//  Discussion:
//
//    Row N contains the N+1 combinatorial coefficients
//
//      C(N,0), C(N,1), C(N,2), ... C(N,N)
//
//  Discussion:
//
//    The sum of the elements of row N is equal to 2**N.
//
//  Formula:
//
//    C(N,K) = N! / ( K! * (N-K)! )
//
//  First terms:
//
//     N K:0  1   2   3   4   5   6   7  8  9 10
//
//     0   1
//     1   1  1
//     2   1  2   1
//     3   1  3   3   1
//     4   1  4   6   4   1
//     5   1  5  10  10   5   1
//     6   1  6  15  20  15   6   1
//     7   1  7  21  35  35  21   7   1
//     8   1  8  28  56  70  56  28   8  1
//     9   1  9  36  84 126 126  84  36  9  1
//    10   1 10  45 120 210 252 210 120 45 10  1
//
//  Recursion:
//
//    C(N,K) = C(N-1,K-1)+C(N-1,K)
//
//  Special values:
//
//    C(N,0) = C(N,N) = 1
//    C(N,1) = C(N,N-1) = N
//    C(N,N-2) = sum ( 1 <= I <= N ) N
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, bool NEXT, indicates whether this is a call for
//    the 'next' row of the triangle.
//    NEXT = FALSE means this is a startup call.  Row N is desired, but
//    presumably this is a first call, or row N-1 was not computed
//    on the previous call.
//    NEXT = TRUE means this is not the first call, and row N-1 was computed
//    on the previous call.  In this case, much work can be saved
//    by using the information from the previous values of ROW
//    to build the next values.
//
//    Input, int N, the row of the triangle desired.  The triangle
//    begins with row 0.
//
//    Output, int ROW[N+1], the row of coefficients.
//    ROW(I) = C(N,I-1).
//
{
  int i;
  int j;

  if ( n < 0 )
  {
    return;
  }

  if ( next )
  {
    for ( i = n-1; 1 <= i; i-- )
    {
      row[i] = row[i] + row[i-1];
    }

    row[n] = 1;
  }
  else
  {
    row[0] = 1;
    for ( i = 1; i <= n; i++ )
    {
      row[i] = 0;
    }

    for ( j = 1; j <= n; j++ )
    {
      for ( i = j; 1 <= i; i-- )
      {
        row[i] = row[i] + row[i-1];
      }
    }

  }

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

int commul ( int iarray[], int n, int nfact )

//****************************************************************************80
//
//  Purpose:
//
//    COMMUL computes a multinomial combinatorial coefficient.
//
//  Definition:
//
//    The multinomial coefficient is a generalization of the binomial
//    coefficient.  It may be interpreted as the number of combinations of
//    N objects, where IARRAY(1) objects are indistinguishable of type 1,
//    ... and IARRAY(K) are indistinguishable of type NFACT.
//
//  Formula:
//
//    COMMUL = N! / ( IARRAY(1)! IARRAY(2)! ... IARRAY(NFACT)! )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int IARRAY(NFACT).
//    IARRAY contains the NFACT values used in the denominator.
//    Note that the sum of these entries should be N,
//    and that all entries should be nonnegative.
//
//    Input, inte N, determines the numerator.
//
//    Input, int NFACT, the number of factors in the numerator.
//
//    Output, int COMMUL, the value of the multinomial coefficient.
//
{
  double arg;
  double fack;
  double facn;
  int i;
  int isum;
//
  for ( i = 0; i < nfact; i++ )
  {
    if ( iarray[i] < 0 )
    {
      cout << "\n";
      cout << "COMMUL - Fatal error\n";
      cout << "  Entry " << i << " of IARRAY = " << iarray[i] << "\n";
      cout << "  But this value must be nonnegative.\n";
      exit ( 1 );
    }
  }

  isum = 0;
  for ( i = 0; i < nfact; i++ )
  {
    isum = isum + iarray[i];
  }

  if ( isum != n )
  {
    cout << "\n";
    cout << "COMMUL - Fatal error!\n";
    cout << "  The sum of the IARRAY entries is " << isum << "\n";
    cout << "  But it must equal N = " << n << "\n";
    exit ( 1 );
  }

  arg = ( double ) ( n + 1 );
  facn = gamma_log ( arg );

  for ( i = 0; i < nfact; i++ )
  {
    arg = ( double ) ( iarray[i] + 1 );
    fack = gamma_log ( arg );
    facn = facn - fack;
  }

  return ( r8_nint ( exp ( facn ) ) );
}
//****************************************************************************80

double cot ( double angle )

//****************************************************************************80
//
//  Purpose:
//
//    COT returns the cotangent of an angle.
//
//  Definition:
//
//    COT ( THETA ) = COS ( THETA ) / SIN ( THETA )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ANGLE, the angle, in radians.
//
//    Output, double COT, the cotangent of the angle.
//
{
  double value;

  value = cos ( angle ) / sin ( angle );

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

double cotd ( double angle )

//****************************************************************************80
//
//  Purpose:
//
//    COTD returns the cotangent of an angle given in degrees.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double ANGLE, the angle, in degrees.
//
//    Output, double COTD, the cotangent of the angle.
//
{
  double value;

  angle = r8_pi ( ) * angle / 180.0;

  value = cos ( angle ) / sin ( angle );

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

double csc ( double theta )

//****************************************************************************80
//
//  Purpose:
//
//    CSC returns the cosecant of X.
//
//  Definition:
//
//    CSC ( THETA ) = 1.0 / SIN ( THETA )
//
//  Discussion:
//
//    CSC is not a built-in library function, and occasionally it
//    is handier, or more concise, to be able to refer to it directly
//    rather than through its definition in terms of the sine function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double THETA, the angle, in radians, whose cosecant is desired.
//    It must be the case that SIN ( THETA ) is not zero.
//
//    Output, double CSC, the cosecant of THETA.
//
{
  double value;

  value = 1.0 / sin ( theta );

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

double e_constant ( void )

//****************************************************************************80
//
//  Purpose:
//
//    E_CONSTANT returns the value of the base of the natural logarithm system.
//
//  Definition:
//
//    E = Limit ( N -> Infinity ) ( 1 + 1 / N )**N
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double E_CONSTANT, the base of the natural logarithm system.
//
{
  return 2.718281828459045235360287;
}
//****************************************************************************80

void erf_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    ERF_VALUES returns some values of the ERF or "error" function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and
//    N_DATA is set to the index of the test data.  On each subsequent
//    call, N_DATA is incremented and that test data is returned.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, double *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 21

  double bvec[N_MAX] = {
    0.0000000000, 0.1124629160, 0.2227025892, 0.3286267595, 
    0.4283923550, 0.5204998778, 0.6038560908, 0.6778011938, 
    0.7421009647, 0.7969082124, 0.8427007929, 0.8802050696, 
    0.9103139782, 0.9340079449, 0.9522851198, 0.9661051465, 
    0.9763483833, 0.9837904586, 0.9890905016, 0.9927904292, 
    0.9953222650 };
  double xvec[N_MAX] = {
    0.0, 0.1, 0.2, 0.3, 
    0.4, 0.5, 0.6, 0.7, 
    0.8, 0.9, 1.0, 1.1, 
    1.2, 1.3, 1.4, 1.5, 
    1.6, 1.7, 1.8, 1.9, 
    2.0 };

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *x = xvec[*n_data];
    *fx = bvec[*n_data];
    *n_data = *n_data + 1;
  }

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

double error_f ( double x )

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

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

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

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

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

    erfx = ( 0.5 - erfx ) + 0.5;

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

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

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

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

    }

  }

  return erfx;
}
//****************************************************************************80

double euler_constant ( void )

//****************************************************************************80
//
//  Purpose:
//
//    EULER_CONSTANT returns the value of the Euler-Mascheroni constant.
//
//  Discussion:
//
//    The Euler-Mascheroni constant is often denoted by a lower-case
//    Gamma.  Gamma is defined as
//
//      Gamma = limit ( M -> Infinity )
//        ( sum ( 1 <= N <= M ) 1 / N ) - log ( M )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double EULER_CONSTANT, the value of the Euler-Mascheroni constant.
//
{
  return ( 0.577215664901532860606512090082402431042 );
}
//****************************************************************************80

void euler_number ( int n, int e[] )

//****************************************************************************80
//
//  Purpose:
//
//    EULER_NUMBER computes the Euler numbers.
//
//  Discussion:
//
//    The Euler numbers can be evaluated in Mathematica with the call
//
//      EulerE[n]
//
//    These numbers rapidly get too big to store in an ordinary integer!
//
//    The terms of odd index are 0.
//
//    E(N) = -C(N,N-2) * E(N-2) - C(N,N-4) * E(N-4) - ... - C(N,0) * E(0).
//
//  First terms:
//
//    E0  = 1
//    E1  = 0
//    E2  = -1
//    E3  = 0
//    E4  = 5
//    E5  = 0
//    E6  = -61
//    E7  = 0
//    E8  = 1385
//    E9  = 0
//    E10 = -50521
//    E11 = 0
//    E12 = 2702765
//    E13 = 0
//    E14 = -199360981
//    E15 = 0
//    E16 = 19391512145
//    E17 = 0
//    E18 = -2404879675441
//    E19 = 0
//    E20 = 370371188237525
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Parameters:
//
//    Input, int N, the index of the last Euler number to compute.
//
//    Output, int E[N+1], the Euler numbers from index 0 to N.
//
{
  int cnk;
  int i;
  int j;
  int sgn;

  if ( n < 0 )
  {
    return;
  }

  e[0] = 1;

  if ( n == 0 )
  {
    return;
  }

  e[1] = 0;

  if ( n == 1 )
  {
    return;
  }

  e[2] = -1;

  for ( i = 3; i <= n; i++ )
  {
    e[i] = 0;

    if ( ( i % 2 ) == 0 )
    {
      for ( j = 2; j <= i; j = j + 2 )
      {
        e[i] = e[i] - i4_choose ( i, j ) * e[i-j];
      }
    }
  }
  return;
}
//****************************************************************************80

double euler_number2 ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    EULER_NUMBER2 computes the Euler numbers.
//
//  Discussion:
//
//    The Euler numbers can be evaluated in Mathematica with the call
//
//      EulerE[n]
//
//  First terms:
//
//    E0  = 1
//    E1  = 0
//    E2  = -1
//    E3  = 0
//    E4  = 5
//    E5  = 0
//    E6  = -61
//    E7  = 0
//    E8  = 1385
//    E9  = 0
//    E10 = -50521
//    E11 = 0
//    E12 = 2702765
//    E13 = 0
//    E14 = -199360981
//    E15 = 0
//    E16 = 19391512145
//    E17 = 0
//    E18 = -2404879675441
//    E19 = 0
//    E20 = 370371188237525
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Parameters:
//
//    Input, int N, the index of the Euler number to compute.
//
//    Output, double EULER_NUMBER2, the value of E(N).
//
{
  double e[7] = { 1.0, -1.0, 5.0, -61.0, 1385.0, 
       -50521.0, 2702765.0 };
  int i;
  int itmax = 1000;
  double sum1;
  double term;
  double value;

  if ( n < 0 )
  {
    return 0.0;
  }

  if ( n == 0 )
  {
    return e[0];
  }

  if ( ( n % 2 ) == 1 )
  {
    return 0.0;
  }

  if ( n <= 12 )
  {
    return e[n/2];
  }

  sum1 = 0.0;

  for ( i = 1; i <= itmax; i++ )
  {
    term = 1.0 / pow ( ( double ) ( 2 * i - 1 ), n + 1 );

    if ( ( i % 2 ) == 1 )
    {
      sum1 = sum1 + term;
    }
    else
    {
      sum1 = sum1 - term;
    }

    if ( fabs ( term ) < 1.0E-10 )
    {
      break;
    }
    else if ( fabs ( term ) < 1.0E-08 * fabs ( sum1 ) )
    {
      break;
    }

  }

  value = pow ( 2.0, n + 2 ) * sum1 * r8_factorial ( n ) 
    / pow ( r8_pi ( ), n + 1 );

  if ( ( n % 4 ) != 0 )
  {
    value = -value;
  }

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

void euler_number_values ( int *n_data, int *n, int *c )

//****************************************************************************80
//
//  Purpose:
//
//    EULER_NUMBER_VALUES returns some values of the Euler numbers.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and N_DATA
//    is set to 1.  On each subsequent call, the input value of N_DATA is
//    incremented and that test data item is returned, if available.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, int *N, the order of the Euler number.
//
//    Output, int *C, the value of the Euler number.
//
{
# define N_MAX 8
 
  int c_vec[N_MAX] = {
    1, 0, -1, 5, 61, 1385, -50521, 2702765 };

  int n_vec[N_MAX] = {
     0, 1, 2, 4, 6, 8, 10, 12 };

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *n = 0;
    *c = 0;
  }
  else
  {
    *n = n_vec[*n_data];
    *c = c_vec[*n_data];
    *n_data = *n_data + 1;
  }

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

double euler_poly ( int n, double x )

//****************************************************************************80
//
//  Purpose:
//
//    EULER_POLY evaluates the N-th Euler polynomial at X.
//
//  First values:
//
//    E(0,X) = 1
//    E(1,X) = X - 1/2
//    E(2,X) = X**2 - X
//    E(3,X) = X**3 - 3/2 X**2 + 1/4
//    E(4,X) = X**4 - 2*X**3 + X
//    E(5,X) = X**5 - 5/2 X**4 + 5/2 X**2 - 1/2
//    E(6,X) = X**6 - 3 X**5 + 5 X**3 - 3 X
//    E(7,X) = X**7 - 7/2 X**6 + 35/4 X**4 - 21/2 X**2 + 17/8
//    E(8,X) = X**8 - 4 X**7 + 14 X**5 - 28 X**3 + 17 X
//
//  Special values:
//
//    E'(N,X) = N * E(N-1,X)
//
//    E(N,1/2) = E(N) / 2**N, where E(N) is the N-th Euler number.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    05 February 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the order of the Euler polynomial to
//    be evaluated.  N must be 0 or greater.
//
//    Input, double X, the value at which the polynomial is to
//    be evaluated.
//
//    Output, double EULER_POLY, the value of E(N,X).
//
{
  double bx1;
  double bx2;
  double value;

  bx1 = bernoulli_poly2 ( n+1, x );
  bx2 = bernoulli_poly2 ( n+1, 0.5 * x );

  value = 2.0 * ( bx1 - bx2 * pow ( ( double ) 2, n+1 ) ) 
    / ( double ) ( n + 1 );

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

void eulerian ( int n, int e[] )

//****************************************************************************80
//
//  Purpose:
//
//    EULERIAN computes the Eulerian number E(N,K).
//
//  Definition:
//
//    A run in a permutation is a sequence of consecutive ascending values.
//
//    E(N,K) is the number of permutations of N objects which contain
//    exactly K runs.
//
//  Examples:
//
//     N = 7
//
//     1     0     0     0     0     0     0
//     1     1     0     0     0     0     0
//     1     4     1     0     0     0     0
//     1    11    11     1     0     0     0
//     1    26    66    26     1     0     0
//     1    57   302   302    57     1     0
//     1   120  1191  2416  1191   120     1
//
//  Recursion:
//
//    E(N,K) = K * E(N-1,K) + (N-K+1) * E(N-1,K-1).
//
//  Properties:
//
//    E(N,1) = E(N,N) = 1.
//    E(N,K) = 0 if K <= 0 or N < K.
//    sum ( 1 <= K <= N ) E(N,K) = N!.
//    X**N = sum ( 0 <= K <= N ) COMB(X+K-1, N ) E(N,K)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Dennis Stanton and Dennis White,
//    Constructive Combinatorics,
//    Springer Verlag, 1986
//
//  Parameters:
//
//    Input, int N, the number of rows desired.
//
//    Output, int E[N*N], the first N rows of Eulerian numbers.
//
{
  int i;
  int j;

  if ( n < 1 )
  {
    return;
  }
//
//  Construct rows 1, 2, ..., N of the Eulerian triangle.
//
  e[1-1+(1-1)*n] = 1;
  for ( j = 2; j <= n; j++ )
  {
    e[1-1+(j-1)*n] = 0;
  }

  for ( i = 2; i <= n; i++ )
  {
    e[i-1+(1-1)*n] = 1;
    for ( j = 2; j <= n; j++ )
    {
      e[i-1+(j-1)*n] = j * e[i-2+(j-1)*n] + ( i - j + 1 ) * e[i-2+(j-2)*n];
    }
  }

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

int f_hofstadter ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    F_HOFSTADTER computes the Hofstadter F sequence.
//
//  Discussion:
//
//    F(N) = 0                if N = 0
//         = N - F ( N - 1 ), otherwise.
//
//    F(N) is defined for all nonnegative integers, and turns out
//    to be equal to int ( ( N + 1 ) / 2 ).
//
//  Table:
//
//     N  F(N)
//    --  ----
//
//     0     0
//     1     1
//     2     1
//     3     2
//     4     2
//     5     3
//     6     3
//     7     4
//     8     4
//     9     5
//    10     5
//    11     6
//    12     6
//    13     7
//    14     7
//    15     8
//    16     8
//    17     9
//    18     9
//    19    10
//    20    10
//    21    11
//    22    11
//    23    12
//    24    12
//    25    13
//    26    13
//    27    14
//    28    14
//    29    15
//    30    15
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Douglas Hofstadter,
//    Goedel, Escher, Bach,
//    Basic Books, 1979.
//
//  Parameters:
//
//    Input, int N, the argument of the function.
//
//    Output, int F_HOFSTADTER, the value of the function.
//
{
  if ( n <= 0 )
  {
    return 0;
  }
  else
  {
    return ( n - f_hofstadter ( n-1 ) );
  }
}
//****************************************************************************80

int fibonacci_direct ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    FIBONACCI_DIRECT computes the N-th Fibonacci number directly.
//
//  Formula:
//
//      F(N) = ( PHIP**N - PHIM**N ) / sqrt(5)
//
//    where
//
//      PHIP = ( 1 + sqrt(5) ) / 2,
//      PHIM = ( 1 - sqrt(5) ) / 2.
//
//  Example:
//
//     N   F
//    --  --
//     0   0
//     1   1
//     2   1
//     3   2
//     4   3
//     5   5
//     6   8
//     7  13
//     8  21
//     9  34
//    10  55
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the index of the Fibonacci number to compute.
//    N should be nonnegative.
//
//    Output, int FIBONACCI_DIRECT, the value of the N-th Fibonacci number.
//
{
  int f;
  double sqrt5 = 2.236068;
  double phim = ( 1.0 - sqrt5 ) / 2.0;
  double phip = ( 1.0 + sqrt5 ) / 2.0;

  if ( n < 0 )
  {
    f = 0;
  }
  else
  {
    f = r8_nint ( ( pow ( phip, n ) - pow ( phim, n ) ) / sqrt ( 5.0 ) );
  }

  return f;
}
//****************************************************************************80

void fibonacci_floor ( int n, int *f, int *i )

//****************************************************************************80
//
//  Purpose:
//
//    FIBONACCI_FLOOR returns the largest Fibonacci number less or equal to N.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the positive integer whose Fibonacci "floor" is desired.
//
//    Output, int *F, the largest Fibonacci number less than or equal to N.
//
//    Output, int *I, the index of the F.
//
{
  if ( n <= 0 )
  {
    *i = 0;
    *f = 0;
  }
  else
  {
    *i = ( int ) ( 
        log ( 0.5 * ( double ) ( 2 * n + 1 ) * sqrt ( 5.0 ) )
      / log ( 0.5 * ( 1.0 + sqrt ( 5.0 ) ) ) );

    *f = fibonacci_direct ( *i );

    if ( n < *f )
    {
      *i = *i - 1;
      *f = fibonacci_direct ( *i );
    }

  }

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

void fibonacci_recursive ( int n, int f[] )

//****************************************************************************80
//
//  Purpose:
//
//    FIBONACCI_RECURSIVE computes the first N Fibonacci numbers.
//
//  Algebraic equation:
//
//    The 'golden ratio' PHI = (1+sqrt(5))/2 satisfies the equation
//
//      X*X-X-1=0
//
//    which is often written as:
//
//       X        1
//      --- =  ------
//       1      X - 1
//
//    expressing the fact that a rectangle, whose sides are in proportion X:1,
//    is similar to the rotated rectangle after a square of side 1 is removed.
//
//      <----X---->
//
//      +-----*---*
//      |     |   |  1
//      |     |   |
//      +-----*---+
//      <--1-><X-1>
//
//  Formula:
//
//    Let
//
//      PHIP = ( 1 + sqrt(5) ) / 2
//      PHIM = ( 1 - sqrt(5) ) / 2
//
//    Then
//
//      F(N) = ( PHIP**N + PHIM**N ) / sqrt(5)
//
//    Moreover, F(N) can be computed by computing PHIP**N / sqrt(5) and rounding
//    to the nearest whole number.
//
//  First terms:
//
//      1
//      1
//      2
//      3
//      5
//      8
//     13
//     21
//     34
//     55
//     89
//    144
//
//    The 40th number is                  102,334,155.
//    The 50th number is               12,586,269,025.
//    The 100th number is 354,224,848,179,261,915,075.
//
//  Recursion:
//
//    F(1) = 1
//    F(2) = 1
//
//    F(N) = F(N-1) + F(N-2)
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the highest Fibonacci number to compute.
//
//    Output, int F[N], the first N Fibonacci numbers.
//
{
  int i;

  if ( n <= 0 )
  {
    return;
  }

  f[0] = 1;

  if ( n <= 1 )
  {
    return;
  }

  f[1] = 1;

  for ( i = 2; i < n; i++ )
  {
    f[i] = f[i-1] + f[i-2];
  }

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

int g_hofstadter ( int n )

//****************************************************************************80
//
//  Purpose:
//
//    G_HOFSTADTER computes the Hofstadter G sequence.
//
//  Discussion:
//
//    G(N) = 0                      if N = 0
//         = N - G ( G ( N - 1 ) ), otherwise.
//
//    G(N) is defined for all nonnegative integers.
//
//    The value of G(N) turns out to be related to the Zeckendorf
//    representation of N as a sum of non-consecutive Fibonacci numbers.
//    To compute G(N), determine the Zeckendorf representation:
//
//      N = sum ( 1 <= I <= M ) F(I)
//
//    and reduce the index of each Fibonacci number by 1:
//
//      G(N) = sum ( 1 <= I <= M ) F(I-1)
//
//    However, this is NOT how the computation is done in this routine.
//    Instead, a straightforward recursive function call is defined
//    to correspond to the definition of the mathematical function.
//
//  Table:
//
//     N  G(N)  Zeckendorf   Decremented
//    --  ----  ----------   -----------
//
//     1   1    1            1
//     2   1    2            1
//     3   2    3            2
//     4   3    3 + 1        2 + 1
//     5   3    5            3
//     6   4    5 + 1        3 + 1
//     7   4    5 + 2        3 + 1
//     8   5    8            5
//     9   6    8 + 1        5 + 1
//    10   6    8 + 2        5 + 1
//    11   7    8 + 3        5 + 2
//    12   8    8 + 3 + 1    5 + 2 + 1
//    13   8    13           8
//    14   9    13 + 1       8 + 1
//    15   9    13 + 2       8 + 1
//    16  10    13 + 3       8 + 2
//    17  11    13 + 3 + 1   8 + 2 + 1
//    18  11    13 + 5       8 + 3
//    19  12    13 + 5 + 1   8 + 3 + 1
//    20  12    13 + 5 + 2   8 + 3 + 1
//    21  13    21           13
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Douglas Hofstadter,
//    Goedel, Escher, Bach,
//    Basic Books, 1979.
//
//  Parameters:
//
//    Input, int N, the argument of the function.
//
//    Output, int G_HOFSTADTER, the value of the function.
//
{
  if ( n <= 0 )
  {
    return 0;
  }
  else
  {
    return ( n - g_hofstadter ( g_hofstadter ( n-1 ) ) );
  }

}
//****************************************************************************80

double gamma_log ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    GAMMA_LOG calculates the natural logarithm of GAMMA ( X ) for positive X.
//
//  Discussion:
//
//    The program uses rational functions that theoretically approximate
//    LOG(GAMMA(X)) to at least 18 significant decimal digits.  
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    FORTRAN77 original version by William Cody, Laura Stoltz
//    C++ version by John Burkardt
//
//  Reference:
//
//    William Cody and Kenneth Hillstrom,
//    Chebyshev Approximations for the Natural Logarithm of the Gamma Function,
//    Mathematics of Computation,
//    Volume 21, 1967, pages 198-203.
//
//    Kenneth Hillstrom,
//    ANL/AMD Program ANLC366S, DGAMMA/DLGAMA,
//    May 1969.
//
//    Hart, Et. Al.,
//    Computer Approximations,
//    Wiley and sons, New York, 1968.
//
//  Parameters:
//
//    Input, double X, the argument of the Gamma function.  X must be positive.
//
//    Output, double GAMMA_LOG, the logarithm of the Gamma function of X.
//    If X <= 0.0, or if overflow would occur, the program returns the
//    value XINF, the largest representable floating point number.
//
//  Local Parameters:
//
//  BETA   - radix for the floating-point representation.
//
//  MAXEXP - the smallest positive power of BETA that overflows.
//
//  XBIG   - largest argument for which LN(GAMMA(X)) is representable
//           in the machine, i.e., the solution to the equation
//             LN(GAMMA(XBIG)) = BETA**MAXEXP.
//
//  FRTBIG - Rough estimate of the fourth root of XBIG
//
//
//  Approximate values for some important machines are:
//
//                            BETA      MAXEXP         XBIG
//
//  CRAY-1        (S.P.)        2        8191       9.62E+2461
//  Cyber 180/855
//    under NOS   (S.P.)        2        1070       1.72E+319
//  IEEE (IBM/XT,
//    SUN, etc.)  (S.P.)        2         128       4.08E+36
//  IEEE (IBM/XT,
//    SUN, etc.)  (D.P.)        2        1024       2.55D+305
//  IBM 3033      (D.P.)       16          63       4.29D+73
//  VAX D-Format  (D.P.)        2         127       2.05D+36
//  VAX G-Format  (D.P.)        2        1023       1.28D+305
//
//
//                           FRTBIG
//
//  CRAY-1        (S.P.)   3.13E+615
//  Cyber 180/855
//    under NOS   (S.P.)   6.44E+79
//  IEEE (IBM/XT,
//    SUN, etc.)  (S.P.)   1.42E+9
//  IEEE (IBM/XT,
//    SUN, etc.)  (D.P.)   2.25D+76
//  IBM 3033      (D.P.)   2.56D+18
//  VAX D-Format  (D.P.)   1.20D+9
//  VAX G-Format  (D.P.)   1.89D+76
//
{
  double c[7] = {
    -1.910444077728E-03, 
     8.4171387781295E-04, 
    -5.952379913043012E-04, 
     7.93650793500350248E-04, 
    -2.777777777777681622553E-03, 
     8.333333333333333331554247E-02, 
     5.7083835261E-03 };
  double corr;
  double d1 = - 5.772156649015328605195174E-01;
  double d2 =   4.227843350984671393993777E-01;
  double d4 =   1.791759469228055000094023;
  double frtbig = 1.42E+09;
  int i;
  double p1[8] = {
    4.945235359296727046734888, 
    2.018112620856775083915565E+02, 
    2.290838373831346393026739E+03, 
    1.131967205903380828685045E+04, 
    2.855724635671635335736389E+04, 
    3.848496228443793359990269E+04, 
    2.637748787624195437963534E+04, 
    7.225813979700288197698961E+03 };
  double p2[8] = {
    4.974607845568932035012064, 
    5.424138599891070494101986E+02, 
    1.550693864978364947665077E+04, 
    1.847932904445632425417223E+05, 
    1.088204769468828767498470E+06, 
    3.338152967987029735917223E+06, 
    5.106661678927352456275255E+06, 
    3.074109054850539556250927E+06 };
  double p4[8] = {
    1.474502166059939948905062E+04, 
    2.426813369486704502836312E+06, 
    1.214755574045093227939592E+08, 
    2.663432449630976949898078E+09, 
    2.940378956634553899906876E+010,
    1.702665737765398868392998E+011,
    4.926125793377430887588120E+011, 
    5.606251856223951465078242E+011 };
  double pnt68 = 0.6796875;
  double q1[8] = {
    6.748212550303777196073036E+01, 
    1.113332393857199323513008E+03, 
    7.738757056935398733233834E+03, 
    2.763987074403340708898585E+04, 
    5.499310206226157329794414E+04, 
    6.161122180066002127833352E+04, 
    3.635127591501940507276287E+04, 
    8.785536302431013170870835E+03 };
  double q2[8] = {
    1.830328399370592604055942E+02, 
    7.765049321445005871323047E+03, 
    1.331903827966074194402448E+05, 
    1.136705821321969608938755E+06, 
    5.267964117437946917577538E+06, 
    1.346701454311101692290052E+07, 
    1.782736530353274213975932E+07, 
    9.533095591844353613395747E+06 };
  double q4[8] = {
    2.690530175870899333379843E+03, 
    6.393885654300092398984238E+05, 
    4.135599930241388052042842E+07, 
    1.120872109616147941376570E+09, 
    1.488613728678813811542398E+010, 
    1.016803586272438228077304E+011, 
    3.417476345507377132798597E+011, 
    4.463158187419713286462081E+011 };
  double res;
  double sqrtpi = 0.9189385332046727417803297;
  double xbig = 4.08E+36;
  double xden;
  double xm1;
  double xm2;
  double xm4;
  double xnum;
  double xsq;
//
//  Return immediately if the argument is out of range.
//
  if ( x <= 0.0 || xbig < x )
  {
    return r8_huge ( );
  }

  if ( x <= r8_epsilon ( ) )
  {
    res = - log ( x );
  }
  else if ( x <= 1.5 )
  {
    if ( x < pnt68 )
    {
      corr = - log ( x );
      xm1 = x;
    }
    else
    {
      corr = 0.0;
      xm1 = ( x - 0.5 ) - 0.5;
    }

    if ( x <= 0.5 || pnt68 <= x )
    {
      xden = 1.0;
      xnum = 0.0;

      for ( i = 0; i < 8; i++ )
      {
        xnum = xnum * xm1 + p1[i];
        xden = xden * xm1 + q1[i];
      }

      res = corr + ( xm1 * ( d1 + xm1 * ( xnum / xden ) ) );
    }
    else
    {
      xm2 = ( x - 0.5 ) - 0.5;
      xden = 1.0;
      xnum = 0.0;
      for ( i = 0; i < 8; i++ )
      {
        xnum = xnum * xm2 + p2[i];
        xden = xden * xm2 + q2[i];
      }

      res = corr + xm2 * ( d2 + xm2 * ( xnum / xden ) );

    }
  }
  else if ( x <= 4.0 )
  {
    xm2 = x - 2.0;
    xden = 1.0;
    xnum = 0.0;
    for ( i = 0; i < 8; i++ )
    {
      xnum = xnum * xm2 + p2[i];
      xden = xden * xm2 + q2[i];
    }

    res = xm2 * ( d2 + xm2 * ( xnum / xden ) );
  }
  else if ( x <= 12.0 )
  {
    xm4 = x - 4.0;
    xden = - 1.0;
    xnum = 0.0;
    for ( i = 0; i < 8; i++ )
    {
      xnum = xnum * xm4 + p4[i];
      xden = xden * xm4 + q4[i];
    }

    res = d4 + xm4 * ( xnum / xden );
  }
  else
  {
    res = 0.0;

    if ( x <= frtbig )
    {

      res = c[6];
      xsq = x * x;

      for ( i = 0; i < 6; i++ )
      {
        res = res / xsq + c[i];
      }

    }

    res = res / x;
    corr = log ( x );
    res = res + sqrtpi - 0.5 * corr;
    res = res + x * ( corr - 1.0 );

  }

  return res;
}
//****************************************************************************80

void gamma_log_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    GAMMA_LOG_VALUES returns some values of the Log Gamma function.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//  Parameters:
//
//    Input/output, int *N_DATA.
//    On input, if N_DATA is 0, the first test data is returned, and
//    N_DATA is set to the index of the test data.  On each subsequent
//    call, N_DATA isincremented and that test data is returned.  When
//    there is no more test data, N_DATA is set to 0.
//
//    Output, double *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 18

  double bvec[N_MAX] = {
     1.524064183,    0.7966780066,   0.3982337117,  
     0.1520599127,   0.000000000,   -0.04987246543, 
    -0.08537410945, -0.1081747934,  -0.1196128950,  
    -0.1207822040,  -0.1125917658,  -0.09580771625, 
    -0.07108385116, -0.03898428380,  0.000000000,   
    12.80182743,    39.33988571,    71.25704193 };
  double xvec[N_MAX] = {
    0.2,  0.4,  0.6,  0.8, 
    1.0,  1.1,  1.2,  1.3, 
    1.4,  1.5,  1.6,  1.7, 
    1.8,  1.9,  2.0, 10.0, 
   20.0, 30.0 };

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

  if ( N_MAX <= *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *x = xvec[*n_data];
    *fx = bvec[*n_data];
    *n_data = *n_data + 1;
  }

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

void gamma_values ( int *n_data, double *x, double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    GAMMA_VALUES returns some values of the Gamma function.
//
//  Discussion:
//
//    The Gamma function is defined as:
//
//      Gamma(Z) = Integral ( 0 <= T < Infinity) T**(Z-1) exp(-T) dT
//
//    It satisfies the recursion:
//
//      Gamma(X+1) = X * Gamma(X)
//
//    Gamma is undefined for nonpositive integral X.
//    Gamma(0.5) = sqrt(PI)
//    For N a positive integer, Gamma(N+1) = N!, the standard factorial.
//
//    In Mathematica, the function can be evaluated by:
//
//      Gamma[x]
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    20 May 2007
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Parameters:
//
//    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, double *X, the argument of the function.
//
//    Output, double *FX, the value of the function.
//
{
# define N_MAX 25

  double fx_vec[N_MAX] = { 
     -0.3544907701811032E+01,  
     -0.1005871979644108E+03,  
      0.9943258511915060E+02,  
      0.9513507698668732E+01,  
      0.4590843711998803E+01,  
      0.2218159543757688E+01,  
      0.1772453850905516E+01,  
      0.1489192248812817E+01,  
      0.1164229713725303E+01,  
      0.1000000000000000E+01,  
      0.9513507698668732E+00,  
      0.9181687423997606E+00,  
      0.8974706963062772E+00,  
      0.8872638175030753E+00,  
      0.8862269254527580E+00,  
      0.8935153492876903E+00,  
      0.9086387328532904E+00,  
      0.9313837709802427E+00,  
      0.9617658319073874E+00,  
      0.1000000000000000E+01,  
      0.2000000000000000E+01,  
      0.6000000000000000E+01,  
      0.3628800000000000E+06,  
      0.1216451004088320E+18,  
      0.8841761993739702E+31 };

  double x_vec[N_MAX] = { 
     -0.50E+00,  
     -0.01E+00,  
      0.01E+00,  
      0.10E+00,  
      0.20E+00,  
      0.40E+00,  
      0.50E+00,  
      0.60E+00,  
      0.80E+00,  
      1.00E+00,  
      1.10E+00,  
      1.20E+00,  
      1.30E+00,  
      1.40E+00,  
      1.50E+00,  
      1.60E+00,  
      1.70E+00,  
      1.80E+00,  
      1.90E+00,  
      2.00E+00,  
      3.00E+00,  
      4.00E+00,  
     10.00E+00,  
     20.00E+00,  
     30.00E+00 }; 

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

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

void gegenbauer_poly ( int n, double alpha, double x, double cx[] )

//****************************************************************************80
//
//  Purpose:
//
//    GEGENBAUER_POLY computes the Gegenbauer polynomials C(0:N,ALPHA)(X).
//
//  Discussion:
//
//    The Gegenbauer polynomial can be evaluated in Mathematica with
//    the command
//
//      GegenbauerC[n,m,x]
//
//  Differential equation:
//
//    (1-X*X) Y'' - (2 ALPHA + 1) X Y' + N (N + 2 ALPHA) Y = 0
//
//  Recursion:
//
//    C(0,ALPHA,X) = 1,
//    C(1,ALPHA,X) = 2*ALPHA*X
//    C(N,ALPHA,X) = ( 2*(N-1+ALPHA)*X * C(N-1,ALPHA,X) - (N-2+2*ALPHA) * C(N-2,ALPHA,X) )/N
//
//  Restrictions:
//
//    ALPHA must be greater than -0.5.
//
//  Special values:
//
//    If ALPHA = 1, the Gegenbauer polynomials reduce to the Chebyshev
//    polynomials of the second kind.
//
//  Norm:
//
//    Integral ( -1 <= X <= 1 )
//      ( 1 - X**2 )**( ALPHA - 0.5 ) * C(N,ALPHA,X)**2 dX
//
//    = PI * 2**( 1 - 2 * ALPHA ) * Gamma ( N + 2 * ALPHA )
//      / ( N! * ( N + ALPHA ) * ( Gamma ( ALPHA ) )**2 )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Parameters:
//
//    Input, int N, the highest order polynomial to compute.
//    Note that polynomials 0 through N will be computed.
//
//    Input, double ALPHA, a parameter which is part of the definition of
//    the Gegenbauer polynomials.  It must be greater than -0.5.
//
//    Input, double X, the point at which the polynomials are to be evaluated.
//
//    Output, double CX[N+1], the values of the first N+1 Gegenbauer
//    polynomials at the point X.
//
{
  int i;

  if ( alpha <= -0.5 )
  {
    cout << "\n";
    cout << "GEGENBAUER_POLY - Fatal error!\n";
    cout << "  Illegal value of ALPHA = " <<alpha << "\n";
    cout << "  but ALPHA must be greater than -0.5.\n";
    exit ( 1 );
  }

  if ( n < 0 )
  {
    return;
  }

  cx[0] = 1.0;

  if ( n == 0 )
  {
    return;
  }

  cx[1] = 2.0 * alpha * x;

  for ( i = 2; i <= n; i++ )
  {
    cx[i] = ( ( ( double ) ( 2 * i - 2 ) + 2.0 * alpha ) * x * cx[i-1] 
            + ( ( double ) (   - i + 2 ) - 2.0 * alpha )     * cx[i-2] ) 
            /   ( double )       i;
  }

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

void gegenbauer_poly_values ( int *n_data, int *n, double *a, double *x, 
  double *fx )

//****************************************************************************80
//
//  Purpose:
//
//    GEGENBAUER_POLY_VALUES returns some values of the Gegenbauer polynomials.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    06 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Milton Abramowitz, Irene Stegun,
//    Handbook of Mathematical Functions,
//    National Bureau of Standards, 1964,
//    ISBN: 0-486-61272-4,
//    LC: QA47.A34.
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Cambridge University Press, 1999,
//    ISBN: 0-521-64314-7,
//    LC: QA76.95.W65.
//
//  Parameters:
//
//    Input/output, int *N_DATA.  The user sets N_DATA to 0 before the
//    first call.  On each call, the routine increments N_DATA by 1, and
//    returns the corresponding data; when there is no more data, the
//    output value of N_DATA will be 0 again.
//
//    Output, int *N, the order parameter of the function.
//
//    Output, double A, the real parameter of the function.
//
//    Output, double X, the argument of the function.
//
//    Output, double FX, the value of the function.
//
{
# define N_MAX 38

  double a_vec[N_MAX] = {
     0.5,  0.5,  0.5, 
     0.5,  0.5,  0.5, 
     0.5,  0.5,  0.5, 
     0.5,  0.5,  0.0, 
     1.0,  2.0,  3.0, 
     4.0,  5.0,  6.0, 
     7.0,  8.0,  9.0, 
    10.0,  3.0,  3.0, 
     3.0,  3.0,  3.0, 
     3.0,  3.0,  3.0, 
     3.0,  3.0,  3.0, 
     3.0,  3.0,  3.0, 
     3.0,  3.0 };
  double fx_vec[N_MAX] = {
    1.0000000000,   0.2000000000,  -0.4400000000, 
   -0.2800000000,   0.2320000000,   0.3075200000, 
   -0.0805760000,  -0.2935168000,  -0.0395648000, 
    0.2459712000,   0.1290720256,   0.0000000000, 
   -0.3600000000,  -0.0800000000,   0.8400000000, 
    2.4000000000,   4.6000000000,   7.4400000000, 
   10.9200000000,  15.0400000000,  19.8000000000, 
   25.2000000000,  -9.0000000000,  -0.1612800000, 
   -6.6729600000,  -8.3750400000,  -5.5267200000, 
    0.0000000000,   5.5267200000,   8.3750400000, 
    6.6729600000,   0.1612800000,  -9.0000000000, 
  -15.4252800000,  -9.6969600000,  22.4409600000, 
  100.8892800000, 252.0000000000 };
  int n_vec[N_MAX] = {
     0,  1,  2, 
     3,  4,  5, 
     6,  7,  8, 
     9, 10,  2, 
     2,  2,  2, 
     2,  2,  2, 
     2,  2,  2, 
     2,  5,  5, 
     5,  5,  5, 
     5,  5,  5, 
     5,  5,  5, 
     5,  5,  5, 
     5,  5 };
  double x_vec[N_MAX] = {
    0.20,  0.20,  0.20, 
    0.20,  0.20,  0.20, 
    0.20,  0.20,  0.20, 
    0.20,  0.20,  0.40, 
    0.40,  0.40,  0.40, 
    0.40,  0.40,  0.40, 
    0.40,  0.40,  0.40, 
    0.40, -0.50, -0.40, 
   -0.30, -0.20, -0.10, 
    0.00,  0.10,  0.20, 
    0.30,  0.40,  0.50, 
    0.60,  0.70,  0.80, 
    0.90,  1.00 };

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

  *n_data = *n_data + 1;

  if ( N_MAX < *n_data )
  {
    *n_data = 0;
    *n = 0;
    *a = 0.0;
    *x = 0.0;
    *fx = 0.0;
  }
  else
  {
    *n = n_vec[*n_data-1];
    *a = a_vec[*n_data-1];
    *x = x_vec[*n_data-1];
    *fx = fx_vec[*n_data-1];
  }

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

double gud ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    GUD evaluates the Gudermannian function.
//
//  Definition:
//
//    The Gudermannian function relates the hyperbolic and trigonometric
//    functions.  For any argument X, there is a corresponding value
//    GAMMA so that
//
//      sinh(x) = tan(gamma).
//
//    The value GAMMA is called the Gudermannian of X.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    08 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the argument of the Gudermannian.
//
//    Output, double GUD, the value of the Gudermannian.
//
{
  double value;

  value = 2.0 * atan ( tanh ( 0.5 * x ) );

  return value;
}
//****