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

using namespace std;

# include "gegenbauer_polynomial.hpp"

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

bool gegenbauer_alpha_check ( double alpha )

//****************************************************************************80
//
//  Purpose:
//
//    gegenbauer_alpha_check() checks the value of ALPHA.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    30 November 2015
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double ALPHA, a parameter which is part of the definition of
//    the Gegenbauer polynomials.  It must be greater than -0.5.
//
//  Output:
//
//    bool GEGENBAUER_ALPHA_CHECK.
//    TRUE, ALPHA is acceptable.
//    FALSE, ALPHA is not acceptable. 
//
{
  bool check;
  bool squawk;

  squawk = false;

  if ( -0.5 < alpha )
  {
    check = true;
  }
  else
  {
    check = false;
    if ( squawk )
    {
      cerr << "\n";
      cerr << "gegenbauer_alpha_check(): Fatal error!\n";
      cerr << "  Illegal value of ALPHA.\n";
      cerr << "  ALPHA = " << alpha << "\n";
      cerr << "  but ALPHA must be greater than -0.5.\n";
    }
  }

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

void gegenbauer_ek_compute ( int n, double alpha, double x[], double w[] )

//****************************************************************************80
//
//  Purpose:
//
//    gegenbauer_ek_compute() computes a Gauss-Gegenbauer quadrature rule.
//
//  Discussion:
//
//    The integral:
//
//      Integral ( -1 <= X <= 1 ) (1-X^2)^ALPHA * F(X) dX
//
//    The quadrature rule:
//
//      Sum ( 1 <= I <= N ) WEIGHT(I) * F ( XTAB(I) )
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    20 November 2015
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Sylvan Elhay, Jaroslav Kautsky,
//    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of 
//    Interpolatory Quadrature,
//    ACM Transactions on Mathematical Software,
//    Volume 13, Number 4, December 1987, pages 399-415.
//
//  Input:
//
//    int N, the order of the quadrature rule.
//
//    double ALPHA, the exponent of (1-X^2) in the weight.  
//    -1.0 < ALPHA is required.
//
//    double A, B, the left and right endpoints 
//    of the interval.
//
//  Output:
//
//    double X[N], the abscissas.
//
//    double W[N], the weights.
//
{
  double abi;
  double *bj;
  bool check;
  int i;
  double zemu;
//
//  Check N.
//
  if ( n < 1 )
  {
    cerr << "\n";
    cerr << "GEGENBAUER_EK_COMPUTE - Fatal error!\n";
    cerr << "  1 <= N is required.\n";
    exit ( 1 );
  }
//
//  Check ALPHA.
//
  check = gegenbauer_alpha_check ( alpha );
  if ( ! check )
  {
    cerr << "\n";
    cerr << "GEGENBAUER_EK_COMPUTE - Fatal error!\n";
    cerr << "  Illegal value of ALPHA.\n";
    exit ( 1 );
  }
//
//  Define the zero-th moment.
//
  zemu = pow ( 2.0, 2.0 * alpha + 1.0 )
    * tgamma ( alpha + 1.0 )
    * tgamma ( alpha + 1.0 )
    / tgamma ( 2.0 * alpha + 2.0 );
//
//  Define the Jacobi matrix.
//
  for ( i = 0; i < n; i++ )
  {
    x[i] = 0.0;
  }

  bj = new double[n];

  bj[0] = 4.0 * pow ( alpha + 1.0, 2 )
    / ( ( 2.0 * alpha + 3.0 ) * pow ( 2.0 * alpha + 2.0, 2 ) );

  for ( i = 2; i <= n; i++ )
  {
    abi = 2.0 * ( alpha + ( double ) i );
    bj[i-1] = 4.0 * ( double ) ( i ) * pow ( alpha + i, 2 ) * ( 2.0 * alpha + i )
      / ( ( abi - 1.0 ) * ( abi + 1.0 ) * abi * abi );
  }

  for ( i = 0; i < n; i++ )
  {
    bj[i] = sqrt ( bj[i] );
  }

  w[0] = sqrt ( zemu );
  for ( i = 1; i < n; i++ )
  {
    w[i] = 0.0;
  }
//
//  Diagonalize the Jacobi matrix.
//
  imtqlx ( n, x, bj, w );

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

  delete [] bj;

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

double gegenbauer_integral ( int expon, double alpha )

//****************************************************************************80
//
//  Purpose:
//
//    gegenbauer_integral() evaluates the integral of a monomial with Gegenbauer weight.
//
//  Discussion:
//
//    The integral:
//
//      integral ( -1 <= x <= +1 ) x^expon (1-x^2)^alpha dx
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    26 February 2008
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int EXPON, the exponent.
//
//    Input, double ALPHA, the exponent of (1-X^2) in the weight factor.
//
//    Output, double GEGENBAUER_INTEGRAL, the value of the integral.
//
{
  double arg1;
  double arg2;
  double arg3;
  double arg4;
  double c;
  double value;
  double value1;

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

  c = ( double ) ( expon );

  arg1 = - alpha;
  arg2 =   1.0 + c;
  arg3 =   2.0 + alpha + c;
  arg4 = - 1.0;

  value1 = r8_hyper_2f1 ( arg1, arg2, arg3, arg4 );

  value = tgamma ( 1.0 + c ) * 2.0 
    * tgamma ( 1.0 + alpha  ) * value1 
    / tgamma ( 2.0 + alpha  + c );

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

double *gegenbauer_polynomial_value ( int m, int n, double alpha, double x[] )

//****************************************************************************80
//
//  Purpose:
//
//    GEGENBAUER_POLYNOMIAL_VALUE computes the Gegenbauer polynomials C(I,ALPHA)(X).
//
//  Differential equation:
//
//    (1-X*X) Y'' - (2 ALPHA + 1) X Y' + M (M + 2 ALPHA) Y = 0
//
//  Recursion:
//
//    C(0,ALPHA,X) = 1,
//    C(1,ALPHA,X) = 2*ALPHA*X
//    C(M,ALPHA,X) = (  ( 2*M-2+2*ALPHA) * X * C(M-1,ALPHA,X) 
//                    + (  -M+2-2*ALPHA)   *   C(M-2,ALPHA,X) ) / M
//
//  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(M,ALPHA,X)^2 dX
//
//    = PI * 2^( 1 - 2 * ALPHA ) * Gamma ( M + 2 * ALPHA ) 
//      / ( M! * ( M + ALPHA ) * ( Gamma ( ALPHA ) )^2 )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 December 2015
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Stephen Wolfram,
//    The Mathematica Book,
//    Fourth Edition,
//    Wolfram Media / Cambridge University Press, 1999.
//
//  Parameters:
//
//    Input, int M, the highest order polynomial to compute.
//    Note that polynomials 0 through N will be computed.
//
//    Input, int N, the number of evaluation points.
//
//    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[N], the evaluation points.
//
//    Output, double GEGENBAUER_POLYNOMIAL_VALUE(1:M+1,N), the values of 
//    Gegenbauer polynomials 0 through M
//    at the N points X.  
//
{
  double *c;
  bool check;
  int i;
  double i_r8;
  int j;

  check = gegenbauer_alpha_check ( alpha );
  if ( ! check )
  {
    cerr << "\n";
    cerr << "GEGENBAUER_POLYNOMIAL_VALUE - Fatal error!\n";
    cerr << "  Illegal value of ALPHA.\n";
    exit ( 1 );
  }

  c = new double[(m+1)*n];

  if ( m < 0 )
  {
    return c;
  }

  if ( n == 0 )
  {
    return c;
  }

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

  if ( m < 1 )
  {
    return c;
  }

  for ( j = 0; j < n; j++ )
  {
    c[1+j*(m+1)] = 2.0 * alpha * x[j];
  }

  for ( i = 2; i <= m; i++ )
  {
    i_r8 = ( double ) i;
    for ( j = 0; j < n; j++ )
    {
      c[i+j*(m+1)] = (  (     2.0 * i_r8 - 2.0  + 2.0 * alpha ) * x[j] * c[i-1+j*(m+1)]
                     +  (         - i_r8 + 2.0  - 2.0 * alpha ) *        c[i-2+j*(m+1)] )
                     /              i_r8 ;
    }
  }

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

void gegenbauer_polynomial_values ( int &n_data, int &n, double &a, double &x,
  double &fx )

//****************************************************************************80
//
//  Purpose:
//
//    GEGENBAUER_POLYNOMIAL_VALUES returns some values of the Gegenbauer polynomials.
//
//  Discussion:
//
//    The Gegenbauer polynomials are also known as the "spherical
//    polynomials" or "ultraspherical polynomials".
//
//    In Mathematica, the function can be evaluated by:
//
//      GegenbauerC[n,m,x]
//
//  Licensing:
//
//    This code is distributed under the MIT 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

  static double a_vec[N_MAX] = {
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.5E+00,
      0.0E+00,
      1.0E+00,
      2.0E+00,
      3.0E+00,
      4.0E+00,
      5.0E+00,
      6.0E+00,
      7.0E+00,
      8.0E+00,
      9.0E+00,
     10.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00,
      3.0E+00 };

  static double fx_vec[N_MAX] = {
       1.0000000000E+00,
       0.2000000000E+00,
      -0.4400000000E+00,
      -0.2800000000E+00,
       0.2320000000E+00,
       0.3075200000E+00,
      -0.0805760000E+00,
      -0.2935168000E+00,
      -0.0395648000E+00,
       0.2459712000E+00,
       0.1290720256E+00,
       0.0000000000E+00,
      -0.3600000000E+00,
      -0.0800000000E+00,
       0.8400000000E+00,
       2.4000000000E+00,
       4.6000000000E+00,
       7.4400000000E+00,
      10.9200000000E+00,
      15.0400000000E+00,
      19.8000000000E+00,
      25.2000000000E+00,
      -9.0000000000E+00,
      -0.1612800000E+00,
      -6.6729600000E+00,
      -8.3750400000E+00,
      -5.5267200000E+00,
       0.0000000000E+00,
       5.5267200000E+00,
       8.3750400000E+00,
       6.6729600000E+00,
       0.1612800000E+00,
      -9.0000000000E+00,
     -15.4252800000E+00,
      -9.6969600000E+00,
      22.4409600000E+00,
     100.8892800000E+00,
     252.0000000000E+00 };

  static 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 };

  static double x_vec[N_MAX] = {
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.20E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
      0.40E+00,
     -0.50E+00,
     -0.40E+00,
     -0.30E+00,
     -0.20E+00,
     -0.10E+00,
      0.00E+00,
      0.10E+00,
      0.20E+00,
      0.30E+00,
      0.40E+00,
      0.50E+00,
      0.60E+00,
      0.70E+00,
      0.80E+00,
      0.90E+00,
      1.00E+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

void gegenbauer_ss_compute ( int order, double alpha, double xtab[], 
  double weight[] )

//****************************************************************************80
//
//  Purpose:
//
//    GEGENBAUER_SS_COMPUTE computes a Gauss-Gegenbauer quadrature rule.
//
//  Discussion:
//
//    The integral:
//
//      Integral ( -1 <= X <= 1 ) (1-X^2)^ALPHA * F(X) dX
//
//    The quadrature rule:
//
//      Sum ( 1 <= I <= ORDER ) WEIGHT(I) * F ( XTAB(I) )
//
//    Thanks to Janiki Raman for pointing out a problem in an earlier
//    version of the code that occurred when ALPHA was -0.5.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    24 June 2008
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Arthur Stroud, Don Secrest,
//    Gaussian Quadrature Formulas,
//    Prentice Hall, 1966,
//    LC: QA299.4G3S7.
//
//  Parameters:
//
//    Input, int ORDER, the order of the quadrature rule.
//
//    Input, double ALPHA, the exponent of (1-X^2) in the weight.  
//    -1.0 < ALPHA is required.
//
//    Output, double XTAB[ORDER], the abscissas.
//
//    Output, double WEIGHT[ORDER], the weights.
//
{
  double an;
  double *c;
  double cc;
  double delta;
  double dp2;
  int i;
  double p1;
  double prod;
  double r1;
  double r2;
  double r3;
  double temp;
  double x;
//
//  Check ORDER.
//
  if ( order < 1 )
  {
    cerr << "\n";
    cerr << "GEGENBAUER_SS_COMPUTE - Fatal error!\n";
    cerr << "  1 <= ORDER is required.\n";
    exit ( 1 );
  }
  
  c = new double[order];
//
//  Check ALPHA.
//
  if ( alpha <= -1.0 )
  {
    cerr << "\n";
    cerr << "GEGENBAUER_SS_COMPUTE - Fatal error!\n";
    cerr << "  -1.0 < ALPHA is required.\n";
    exit ( 1 );
  }
//
//  Set the recursion coefficients.
//
  c[0] = 0.0;
  if ( 2 <= order )
  {
    c[1] = 1.0 / ( 2.0 * alpha + 3.0 );
  }

  for ( i = 3; i <= order; i++ )
  {
    c[i-1] = ( double ) ( i - 1 ) 
          * ( alpha + alpha + ( double ) ( i - 1 ) ) / 
          ( ( alpha + alpha + ( double ) ( 2 * i - 1 ) ) 
          * ( alpha + alpha + ( double ) ( 2 * i - 3 ) ) );
  }

  delta = tgamma ( alpha         + 1.0 ) 
        * tgamma (         alpha + 1.0 ) 
        / tgamma ( alpha + alpha + 2.0 );

  prod = 1.0;
  for ( i = 2; i <= order; i++ )
  {
    prod = prod * c[i-1];
  }
  cc = delta * pow ( 2.0, alpha + alpha + 1.0 ) * prod;

  for ( i = 1; i <= order; i++ )
  {
    if ( i == 1 )
    {
      an = alpha / ( double ) ( order );

      r1 = ( 1.0 + alpha ) 
        * ( 2.78 / ( 4.0 + ( double ) ( order * order ) ) 
        + 0.768 * an / ( double ) ( order ) );

      r2 = 1.0 + 2.44 * an + 1.282 * an * an;

      x = ( r2 - r1 ) / r2;
    }
    else if ( i == 2 )
    {
      r1 = ( 4.1 + alpha ) / 
        ( ( 1.0 + alpha ) * ( 1.0 + 0.156 * alpha ) );

      r2 = 1.0 + 0.06 * ( ( double ) ( order ) - 8.0 ) * 
        ( 1.0 + 0.12 * alpha ) / ( double ) ( order );

      r3 = 1.0 + 0.012 * alpha * 
        ( 1.0 + 0.25 * fabs ( alpha ) ) / ( double ) ( order );

      x = x - r1 * r2 * r3 * ( 1.0 - x );
    }
    else if ( i == 3 )
    {
      r1 = ( 1.67 + 0.28 * alpha ) / ( 1.0 + 0.37 * alpha );

      r2 = 1.0 + 0.22 * ( ( double ) ( order ) - 8.0 ) 
        / ( double ) ( order );

      r3 = 1.0 + 8.0 * alpha / 
        ( ( 6.28 + alpha ) * ( double ) ( order * order ) );

      x = x - r1 * r2 * r3 * ( xtab[0] - x );
    }
    else if ( i < order - 1 )
    {
      x = 3.0 * xtab[i-2] - 3.0 * xtab[i-3] + xtab[i-4];
    }
    else if ( i == order - 1 )
    {
      r1 = ( 1.0 + 0.235 * alpha ) / ( 0.766 + 0.119 * alpha );

      r2 = 1.0 / ( 1.0 + 0.639 
        * ( ( double ) ( order ) - 4.0 ) 
        / ( 1.0 + 0.71 * ( ( double ) ( order ) - 4.0 ) ) );

      r3 = 1.0 / ( 1.0 + 20.0 * alpha / ( ( 7.5 + alpha ) * 
        ( double ) ( order * order ) ) );

      x = x + r1 * r2 * r3 * ( x - xtab[i-3] );
    }
    else if ( i == order )
    {
      r1 = ( 1.0 + 0.37 * alpha ) / ( 1.67 + 0.28 * alpha );

      r2 = 1.0 / 
        ( 1.0 + 0.22 * ( ( double ) ( order ) - 8.0 ) 
        / ( double ) ( order ) );

      r3 = 1.0 / ( 1.0 + 8.0 * alpha / 
        ( ( 6.28 + alpha ) * ( double ) ( order * order ) ) );

      x = x + r1 * r2 * r3 * ( x - xtab[i-3] );
    }

    gegenbauer_ss_root ( x, order, alpha, dp2, p1, c );

    xtab[i-1] = x;
    weight[i-1] = cc / ( dp2 * p1 );
  }
//
//  Reverse the order of the values.
//
  for ( i = 1; i <= order/2; i++ )
  {
    temp          = xtab[i-1];
    xtab[i-1]     = xtab[order-i];
    xtab[order-i] = temp;
  }

  for ( i = 1; i <=order/2; i++ )
  {
    temp            = weight[i-1];
    weight[i-1]     = weight[order-i];
    weight[order-i] = temp;
  }

  delete [] c;

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

void gegenbauer_ss_recur ( double &p2, double &dp2, double &p1, double x,
  int order, double alpha, double c[] )

//****************************************************************************80
//
//  Purpose:
//
//    GEGENBAUER_SS_RECUR: value and derivative of a Gegenbauer polynomial.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    26 February 2008
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Arthur Stroud, Don Secrest,
//    Gaussian Quadrature Formulas,
//    Prentice Hall, 1966,
//    LC: QA299.4G3S7.
//
//  Parameters:
//
//    Output, double &P2, the value of J(ORDER)(X).
//
//    Output, double &DP2, the value of J'(ORDER)(X).
//
//    Output, double &P1, the value of J(ORDER-1)(X).
//
//    Input, double X, the point at which polynomials are evaluated.
//
//    Input, int ORDER, the order of the polynomial to be computed.
//
//    Input, double ALPHA, the exponents of (1-X^2).
//
//    Input, double C[ORDER], the recursion coefficients.
//
{
  double dp0;
  double dp1;
  int i;
  double p0;

  p1 = 1.0;
  dp1 = 0.0;

  p2 = x;
  dp2 = 1.0;

  for ( i = 2; i <= order; i++ )
  {
    p0 = p1;
    dp0 = dp1;

    p1 = p2;
    dp1 = dp2;

    p2 = x *  p1 - c[i-1] * p0;
    dp2 = x * dp1 + p1 - c[i-1] * dp0;
  }
  return;
}
//****************************************************************************80

void gegenbauer_ss_root ( double &x, int order, double alpha,  double &dp2, 
  double &p1, double c[] )

//****************************************************************************80
//
//  Purpose:
//
//    GEGENBAUER_SS_ROOT improves a root of a Gegenbauer polynomial.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    26 February 2008
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Arthur Stroud, Don Secrest,
//    Gaussian Quadrature Formulas,
//    Prentice Hall, 1966,
//    LC: QA299.4G3S7.
//
//  Parameters:
//
//    Input/output, double &X, the approximate root, which
//    should be improved on output.
//
//    Input, int ORDER, the order of the polynomial to be computed.
//
//    Input, double ALPHA, the exponents of (1-X^2).
//
//    Output, double &DP2, the value of J'(ORDER)(X).
//
//    Output, double &P1, the value of J(ORDER-1)(X).
//
//    Input, double C[ORDER], the recursion coefficients.
//
{
  double d;
  double eps;
  double p2;
  int step;
  int step_max = 10;

  eps = DBL_EPSILON;

  for ( step = 1; step <= step_max; step++ )
  {
    gegenbauer_ss_recur ( p2, dp2, p1, x, order, alpha, c );

    d = p2 / dp2;
    x = x - d;

    if ( fabs ( d ) <= eps * ( fabs ( x ) + 1.0 ) )
    {
      return;
    }
  }
  return;
}
//****************************************************************************80

double *gegenbauer_to_monomial_matrix ( int n, double alpha )

//****************************************************************************80
//
//  Purpose:
//
//    gegenbauer_to_monomial_matrix(): Gegenbauer to monomial conversion matrix.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    01 April 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int n: the order of A.
//
//    double alpha: the parameter.
//
//  Output:
//
//    double A[N,N]: the matrix.
//
{
  double *A;
  double c1;
  double c2;
  int i;
  int j;
  int nn;

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

  if ( n <= 0 )
  {
    return A;
  }

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

  A[0+0*n] = 1.0;

  if ( n == 1 )
  {
    return A;
  }

  A[1+1*n] = 2.0 * alpha;
//
//  Perform convex sum.
//  Translating "(n+1) C(n+1) = 2 (n+alpha) x C(n) - ( n + 2 alpha - 1 ) C(n-1)"
//  drove me nuts, between indexing at 1 rather than 0, and dealing with
//  the interpretation of "n+1", because I now face the rare "off by 2" error!
//
  for ( j = 2; j < n; j++ )
  {
    nn = j - 1;
    c1 = ( 2.0 * nn + 2.0 * alpha       ) / ( nn + 1 );
    c2 = (     - nn - 2.0 * alpha + 1.0 ) / ( nn + 1 );
    for ( i = 1; i <= j; i++ )
    {
      A[i+j*n] = c1 * A[i-1+(j-1)*n];
    }
    for ( i = 0; i <= j - 2; i++ )
    {
      A[i+j*n] = A[i+j*n] + c2 * A[i+(j-2)*n];
    }
  }

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

void hyper_2f1_values ( int &n_data, double &a, double &b, double &c,
  double &x, double &fx )

//****************************************************************************80
//
//  Purpose:
//
//    hyper_2f1_values() returns some values of the hypergeometric function 2F1.
//
//  Discussion:
//
//    In Mathematica, the function can be evaluated by:
//
//      fx = Hypergeometric2F1 [ a, b, c, x ]
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    09 September 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.
//
//    Shanjie Zhang, Jianming Jin,
//    Computation of Special Functions,
//    Wiley, 1996,
//    ISBN: 0-471-11963-6,
//    LC: QA351.C45
//
//    Daniel Zwillinger,
//    CRC Standard Mathematical Tables and Formulae,
//    30th Edition, CRC Press, 1996, pages 651-652.
//
//  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, &C, &X, the parameters of the function.
//
//    Output, double &FX, the value of the function.
//
{
# define N_MAX 24

  static double a_vec[N_MAX] = {
   -2.5,
   -0.5,
    0.5,
    2.5,
   -2.5,
   -0.5,
    0.5,
    2.5,
   -2.5,
   -0.5,
    0.5,
    2.5,
    3.3,
    1.1,
    1.1,
    3.3,
    3.3,
    1.1,
    1.1,
    3.3,
    3.3,
    1.1,
    1.1,
    3.3 };
  static double b_vec[N_MAX] = {
    3.3,
    1.1,
    1.1,
    3.3,
    3.3,
    1.1,
    1.1,
    3.3,
    3.3,
    1.1,
    1.1,
    3.3,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7 };
  static double c_vec[N_MAX] = {
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
    6.7,
   -5.5,
   -0.5,
    0.5,
    4.5,
   -5.5,
   -0.5,
    0.5,
    4.5,
   -5.5,
   -0.5,
    0.5,
    4.5 };
  static double fx_vec[N_MAX] = {
    0.72356129348997784913,
    0.97911109345277961340,
    1.0216578140088564160,
    1.4051563200112126405,
    0.46961431639821611095,
    0.95296194977446325454,
    1.0512814213947987916,
    2.3999062904777858999,
    0.29106095928414718320,
    0.92536967910373175753,
    1.0865504094806997287,
    5.7381565526189046578,
    15090.669748704606754,
   -104.31170067364349677,
    21.175050707768812938,
    4.1946915819031922850,
    1.0170777974048815592E+10,
   -24708.635322489155868,
    1372.2304548384989560,
    58.092728706394652211,
    5.8682087615124176162E+18,
   -4.4635010147295996680E+08,
    5.3835057561295731310E+06,
    20396.913776019659426 };
  static double x_vec[N_MAX] = {
    0.25,
    0.25,
    0.25,
    0.25,
    0.55,
    0.55,
    0.55,
    0.55,
    0.85,
    0.85,
    0.85,
    0.85,
    0.25,
    0.25,
    0.25,
    0.25,
    0.55,
    0.55,
    0.55,
    0.55,
    0.85,
    0.85,
    0.85,
    0.85 };

  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;
    c = 0.0;
    x = 0.0;
    fx = 0.0;
  }
  else
  {
    a = a_vec[n_data-1];
    b = b_vec[n_data-1];
    c = c_vec[n_data-1];
    x = x_vec[n_data-1];
    fx = fx_vec[n_data-1];
  }

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

void imtqlx ( int n, double d[], double e[], double z[] )

//****************************************************************************80
//
//  Purpose:
//
//    IMTQLX diagonalizes a symmetric tridiagonal matrix.
//
//  Discussion:
//
//    This routine is a slightly modified version of the EISPACK routine to 
//    perform the implicit QL algorithm on a symmetric tridiagonal matrix. 
//
//    The authors thank the authors of EISPACK for permission to use this
//    routine. 
//
//    It has been modified to produce the product Q' * Z, where Z is an input 
//    vector and Q is the orthogonal matrix diagonalizing the input matrix.  
//    The changes consist (essentially) of applying the orthogonal transformations
//    directly to Z as they are generated.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    08 January 2010
//
//  Author:
//
//    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
//    This version by John Burkardt.
//
//  Reference:
//
//    Sylvan Elhay, Jaroslav Kautsky,
//    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of 
//    Interpolatory Quadrature,
//    ACM Transactions on Mathematical Software,
//    Volume 13, Number 4, December 1987, pages 399-415.
//
//    Roger Martin, James Wilkinson,
//    The Implicit QL Algorithm,
//    Numerische Mathematik,
//    Volume 12, Number 5, December 1968, pages 377-383.
//
//  Parameters:
//
//    Input, int N, the order of the matrix.
//
//    Input/output, double D(N), the diagonal entries of the matrix.
//    On output, the information in D has been overwritten.
//
//    Input/output, double E(N), the subdiagonal entries of the 
//    matrix, in entries E(1) through E(N-1).  On output, the information in
//    E has been overwritten.
//
//    Input/output, double Z(N).  On input, a vector.  On output,
//    the value of Q' * Z, where Q is the matrix that diagonalizes the
//    input symmetric tridiagonal matrix.
//
{
  double b;
  double c;
  double f;
  double g;
  int i;
  int ii;
  int itn = 30;
  int j;
  int k;
  int l;
  int m;
  int mml;
  double p;
  double prec;
  double r;
  double s;

  prec = DBL_EPSILON;

  if ( n == 1 )
  {
    return;
  }

  e[n-1] = 0.0;

  for ( l = 1; l <= n; l++ )
  {
    j = 0;
    for ( ; ; )
    {
      for ( m = l; m <= n; m++ )
      {
        if ( m == n )
        {
          break;
        }

        if ( fabs ( e[m-1] ) <= prec * ( fabs ( d[m-1] ) + fabs ( d[m] ) ) )
        {
          break;
        }
      }
      p = d[l-1];
      if ( m == l )
      {
        break;
      }
      if ( itn <= j )
      {
        cerr << "\n";
        cerr << "IMTQLX - Fatal error!\n";
        cerr << "  Iteration limit exceeded\n";
        exit ( 1 );
      }
      j = j + 1;
      g = ( d[l] - p ) / ( 2.0 * e[l-1] );
      r = sqrt ( g * g + 1.0 );
      g = d[m-1] - p + e[l-1] / ( g + fabs ( r ) * r8_sign ( g ) );
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ( ii = 1; ii <= mml; ii++ )
      {
        i = m - ii;
        f = s * e[i-1];
        b = c * e[i-1];

        if ( fabs ( g ) <= fabs ( f ) )
        {
          c = g / f;
          r = sqrt ( c * c + 1.0 );
          e[i] = f * r;
          s = 1.0 / r;
          c = c * s;
        }
        else
        {
          s = f / g;
          r = sqrt ( s * s + 1.0 );
          e[i] = g * r;
          c = 1.0 / r;
          s = s * c;
        }
        g = d[i] - p;
        r = ( d[i-1] - g ) * s + 2.0 * c * b;
        p = s * r;
        d[i] = g + p;
        g = c * r - b;
        f = z[i];
        z[i] = s * z[i-1] + c * f;
        z[i-1] = c * z[i-1] - s * f;
      }
      d[l-1] = d[l-1] - p;
      e[l-1] = g;
      e[m-1] = 0.0;
    }
  }
//
//  Sorting.
//
  for ( ii = 2; ii <= m; ii++ )
  {
    i = ii - 1;
    k = i;
    p = d[i-1];

    for ( j = ii; j <= n; j++ )
    {
      if ( d[j-1] < p )
      {
         k = j;
         p = d[j-1];
      }
    }

    if ( k != i )
    {
      d[k-1] = d[i-1];
      d[i-1] = p;
      p = z[i-1];
      z[i-1] = z[k-1];
      z[k-1] = p;
    }
  }
  return;
}
//****************************************************************************80

double r8_hyper_2f1 ( double a, double b, double c, double x )

//****************************************************************************80
//
//  Purpose:
//
//    R8_HYPER_2F1 evaluates the hypergeometric function 2F1(A,B,C,X).
//
//  Discussion:
//
//    A minor bug was corrected.  The HW variable, used in several places as
//    the "old" value of a quantity being iteratively improved, was not
//    being initialized.  JVB, 11 February 2008.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    05 July 2009
//
//  Author:
//
//    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
//    This version by John Burkardt.
//
//    The original FORTRAN77 version of this routine is copyrighted by
//    Shanjie Zhang and Jianming Jin.  However, they give permission to
//    incorporate this routine into a user program provided that the copyright
//    is acknowledged.
//
//  Reference:
//
//    Shanjie Zhang, Jianming Jin,
//    Computation of Special Functions,
//    Wiley, 1996,
//    ISBN: 0-471-11963-6,
//    LC: QA351.C45
//
//  Parameters:
//
//    Input, double A, B, C, X, the arguments of the function.
//    C must not be equal to a nonpositive integer.
//    X < 1.
//
//    Output, double R8_HYPER_2F1, the value of the function.
//
{
  double a0;
  double aa;
  double bb;
  double c0;
  double c1;
  double el = 0.5772156649015329;
  double eps;
  double f0;
  double f1;
  double g0;
  double g1;
  double g2;
  double g3;
  double ga;
  double gabc;
  double gam;
  double gb;
  double gbm;
  double gc;
  double gca;
  double gcab;
  double gcb;
  double gm;
  double hf;
  double hw;
  int j;
  int k;
  bool l0;
  bool l1;
  bool l2;
  bool l3;
  bool l4;
  bool l5;
  int m;
  int nm;
  double pa;
  double pb;
  const double r8_pi = 3.141592653589793;
  double r;
  double r0;
  double r1;
  double rm;
  double rp;
  double sm;
  double sp;
  double sp0;
  double x1;

  l0 = ( c == ( int ) ( c ) ) && ( c < 0.0 );
  l1 = ( 1.0 - x < 1.0E-15 ) && ( c - a - b <= 0.0 );
  l2 = ( a == ( int ) ( a ) ) && ( a < 0.0 );
  l3 = ( b == ( int ) ( b ) ) && ( b < 0.0 );
  l4 = ( c - a == ( int ) ( c - a ) ) && ( c - a <= 0.0 );
  l5 = ( c - b == ( int ) ( c - b ) ) && ( c - b <= 0.0 );

  if ( l0 )
  {
    cerr << "\n";
    cerr << "R8_HYPER_2F1 - Fatal error!\n";
    cerr << "  The hypergeometric series is divergent.\n";
    cerr << "  C is integral and negative.\n";
    cerr << "  C = " << c << "\n";
    exit ( 1 );
  }

  if ( l1 )
  {
    cerr << "\n";
    cerr << "R8_HYPER_2F1 - Fatal error!\n";
    cerr << "  The hypergeometric series is divergent.\n";
    cerr << "  1 - X < 0, C - A - B <= 0\n";
    cerr << "  A = " << a << "\n";
    cerr << "  B = " << b << "\n";
    cerr << "  C = " << c << "\n";
    cerr << "  X = " << x << "\n";
    exit ( 1 );
  }

  if ( 0.95 < x )
  {
    eps = 1.0E-08;
  }
  else
  {
    eps = 1.0E-15;
  }

  if ( x == 0.0 || a == 0.0 || b == 0.0 )
  {
    hf = 1.0;
    return hf;
  }
  else if ( 1.0 - x == eps && 0.0 < c - a - b )
  {
    gc = tgamma ( c );
    gcab = tgamma ( c - a - b );
    gca = tgamma ( c - a );
    gcb = tgamma ( c - b );
    hf = gc * gcab / ( gca * gcb );
    return hf;
  }
  else if ( 1.0 + x <= eps && fabs ( c - a + b - 1.0 ) <= eps )
  {
    g0 = sqrt ( r8_pi ) * pow ( 2.0, - a );
    g1 = tgamma ( c );
    g2 = tgamma ( 1.0 + a / 2.0 - b );
    g3 = tgamma ( 0.5 + 0.5 * a );
    hf = g0 * g1 / ( g2 * g3 );
    return hf;
  }
  else if ( l2 || l3 )
  {
    if ( l2 )
    {
      nm = ( int ) ( fabs ( a ) );
    }

    if ( l3 )
    {
      nm = ( int ) ( fabs ( b ) );
    }

    hf = 1.0;
    r = 1.0;

    for ( k = 1; k <= nm; k++ )
    {
      r = r * ( a + k - 1.0 ) * ( b + k - 1.0 ) 
        / ( k * ( c + k - 1.0 ) ) * x;
      hf = hf + r;
    }

    return hf;
  }
  else if ( l4 || l5 )
  {
    if ( l4 )
    {
      nm = ( int ) ( fabs ( c - a ) );
    }

    if ( l5 )
    {
      nm = ( int ) ( fabs ( c - b ) );
    }

    hf = 1.0;
    r  = 1.0;
    for ( k = 1; k <= nm; k++ )
    {
      r = r * ( c - a + k - 1.0 ) * ( c - b + k - 1.0 ) 
        / ( k * ( c + k - 1.0 ) ) * x;
      hf = hf + r;
    }
    hf = pow ( 1.0 - x, c - a - b ) * hf;
    return hf;
  }

  aa = a;
  bb = b;
  x1 = x;

  if ( x < 0.0 )
  {
    x = x / ( x - 1.0 );
    if ( a < c && b < a && 0.0 < b )
    {
      a = bb;
      b = aa;
    }
    b = c - b;
  }

  if ( 0.75 <= x )
  {
    gm = 0.0;

    if ( fabs ( c - a - b - ( int ) ( c - a - b ) ) < 1.0E-15 )
    {
      m = ( int ) ( c - a - b );
      ga = tgamma ( a );
      gb = tgamma ( b );
      gc = tgamma ( c );
      gam = tgamma ( a + m );
      gbm = tgamma ( b + m );

      pa = r8_psi ( a );
      pb = r8_psi ( b );

      if ( m != 0 )
      {
        gm = 1.0;
      }

      for ( j = 1; j <= abs ( m ) - 1; j++ )
      {
        gm = gm * j;
      }

      rm = 1.0;
      for ( j = 1; j <= abs ( m ); j++ )
      {
        rm = rm * j;
      }

      f0 = 1.0;
      r0 = 1.0;;
      r1 = 1.0;
      sp0 = 0.0;;
      sp = 0.0;

      if ( 0 <= m )
      {
        c0 = gm * gc / ( gam * gbm );
        c1 = - gc * pow ( x - 1.0, m ) / ( ga * gb * rm );

        for ( k = 1; k <= m - 1; k++ )
        {
          r0 = r0 * ( a + k - 1.0 ) * ( b + k - 1.0 ) 
            / ( k * ( k - m ) ) * ( 1.0 - x );
          f0 = f0 + r0;
        }

        for ( k = 1; k <= m; k++ )
        {
          sp0 = sp0 + 1.0 / ( a + k - 1.0 ) + 1.0 / ( b + k - 1.0 ) 
          - 1.0 / ( double ) ( k );
        }

        f1 = pa + pb + sp0 + 2.0 * el + log ( 1.0 - x );
        hw = f1;

        for ( k = 1; k <= 250; k++ )
        {
          sp = sp + ( 1.0 - a ) / ( k * ( a + k - 1.0 ) ) 
            + ( 1.0 - b ) / ( k * ( b + k - 1.0 ) );

          sm = 0.0;
          for ( j = 1; j <= m; j++ )
          {
            sm = sm + ( 1.0 - a ) 
              / ( ( j + k ) * ( a + j + k - 1.0 ) ) 
              + 1.0 / ( b + j + k - 1.0 );
          }

          rp = pa + pb + 2.0 * el + sp + sm + log ( 1.0 - x );

          r1 = r1 * ( a + m + k - 1.0 ) * ( b + m + k - 1.0 ) 
            / ( k * ( m + k ) ) * ( 1.0 - x );

          f1 = f1 + r1 * rp;

          if ( fabs ( f1 - hw ) < fabs ( f1 ) * eps )
          {
            break;
          }
          hw = f1;
        }
        hf = f0 * c0 + f1 * c1;
      }
      else if ( m < 0 )
      {
        m = - m;
        c0 = gm * gc / ( ga * gb * pow ( 1.0 - x, m ) );
        c1 = - pow ( - 1.0, m ) * gc / ( gam * gbm * rm );

        for ( k = 1; k <= m - 1; k++ )
        {
          r0 = r0 * ( a - m + k - 1.0 ) * ( b - m + k - 1.0 ) 
            / ( k * ( k - m ) ) * ( 1.0 - x );
          f0 = f0 + r0;
        }

        for ( k = 1; k <= m; k++ )
        {
          sp0 = sp0 + 1.0 / ( double ) ( k );
        }

        f1 = pa + pb - sp0 + 2.0 * el + log ( 1.0 - x );
        hw = f1;

        for ( k = 1; k <= 250; k++ )
        {
          sp = sp + ( 1.0 - a ) 
            / ( k * ( a + k - 1.0 ) ) 
            + ( 1.0 - b ) / ( k * ( b + k - 1.0 ) );

          sm = 0.0;
          for ( j = 1; j <= m; j++ )
          {
            sm = sm + 1.0 / ( double ) ( j + k );
          }

          rp = pa + pb + 2.0 * el + sp - sm + log ( 1.0 - x );

          r1 = r1 * ( a + k - 1.0 ) * ( b + k - 1.0 ) 
            / ( k * ( m + k ) ) * ( 1.0 - x );

          f1 = f1 + r1 * rp;

          if ( fabs ( f1 - hw ) < fabs ( f1 ) * eps )
          {
            break;
          }

          hw = f1;
        }

        hf = f0 * c0 + f1 * c1;
      }
    }
    else
    {
      ga = tgamma ( a );
      gb = tgamma ( b );
      gc = tgamma ( c );
      gca = tgamma ( c - a );
      gcb = tgamma ( c - b );
      gcab = tgamma ( c - a - b );
      gabc = tgamma ( a + b - c );
      c0 = gc * gcab / ( gca * gcb );
      c1 = gc * gabc / ( ga * gb ) * pow ( 1.0 - x, c - a - b );
      hf = 0.0;
      hw = hf;
      r0 = c0;
      r1 = c1;

      for ( k = 1; k <= 250; k++ )
      {
        r0 = r0 * ( a + k - 1.0 ) * ( b + k - 1.0 ) 
          / ( k * ( a + b - c + k ) ) * ( 1.0 - x );

        r1 = r1 * ( c - a + k - 1.0 ) * ( c - b + k - 1.0 ) 
          / ( k * ( c - a - b + k ) ) * ( 1.0 - x );

        hf = hf + r0 + r1;

        if ( fabs ( hf - hw ) < fabs ( hf ) * eps )
        {
          break;
        }
        hw = hf;
      }
      hf = hf + c0 + c1;
    }
  }
  else
  {
    a0 = 1.0;

    if ( a < c && c < 2.0 * a && b < c && c < 2.0 * b )
    {
      a0 = pow ( 1.0 - x, c - a - b );
      a = c - a;
      b = c - b;
    }

    hf = 1.0;
    hw = hf;
    r = 1.0;

    for ( k = 1; k <= 250; k++ )
    {
      r = r * ( a + k - 1.0 ) * ( b + k - 1.0 ) 
        / ( k * ( c + k - 1.0 ) ) * x;

      hf = hf + r;

      if ( fabs ( hf - hw ) <= fabs ( hf ) * eps )
      {
        break;
      }

      hw = hf;
    }
    hf = a0 * hf;
  }

  if ( x1 < 0.0 )
  {
    x = x1;
    c0 = 1.0 / pow ( 1.0 - x, aa );
    hf = c0 * hf;
  }

  a = aa;
  b = bb;

  if ( 120 < k )
  {
    cout << "\n";
    cout << "R8_HYPER_2F1 - Warning!\n";
    cout << "  A large number of iterations were needed.\n";
    cout << "  The accuracy of the results should be checked.\n";
  }

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

double r8_psi ( double xx )

//****************************************************************************80
//
//  Purpose:
//
//    R8_PSI evaluates the function Psi(X).
//
//  Discussion:
//
//    This routine evaluates the logarithmic derivative of the
//    Gamma function,
//
//      PSI(X) = d/dX ( GAMMA(X) ) / GAMMA(X)
//             = d/dX LN ( GAMMA(X) )
//
//    for real X, where either
//
//      - XMAX1 < X < - XMIN, and X is not a negative integer,
//
//    or
//
//      XMIN < X.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    09 February 2008
//
//  Author:
//
//    Original FORTRAN77 version by William Cody.
//    This version by John Burkardt.
//
//  Reference:
//
//    William Cody, Anthony Strecok, Henry Thacher,
//    Chebyshev Approximations for the Psi Function,
//    Mathematics of Computation,
//    Volume 27, Number 121, January 1973, pages 123-127.
//
//  Parameters:
//
//    Input, double XX, the argument of the function.
//
//    Output, double R8_PSI, the value of the function.
//
{
  double aug;
  double den;
  double four = 4.0;
  double fourth = 0.25;
  double half = 0.5;
  int i;
  int n;
  int nq;
  double one = 1.0;
  double p1[9] = { 
   4.5104681245762934160E-03, 
   5.4932855833000385356, 
   3.7646693175929276856E+02, 
   7.9525490849151998065E+03, 
   7.1451595818951933210E+04, 
   3.0655976301987365674E+05, 
   6.3606997788964458797E+05, 
   5.8041312783537569993E+05, 
   1.6585695029761022321E+05 };
  double p2[7] = { 
  -2.7103228277757834192, 
  -1.5166271776896121383E+01, 
  -1.9784554148719218667E+01, 
  -8.8100958828312219821, 
  -1.4479614616899842986, 
  -7.3689600332394549911E-02, 
  -6.5135387732718171306E-21 };
  double piov4 = 0.78539816339744830962;
  double q1[8] = { 
   9.6141654774222358525E+01, 
   2.6287715790581193330E+03, 
   2.9862497022250277920E+04, 
   1.6206566091533671639E+05, 
   4.3487880712768329037E+05, 
   5.4256384537269993733E+05, 
   2.4242185002017985252E+05, 
   6.4155223783576225996E-08 };
  double q2[6] = { 
   4.4992760373789365846E+01, 
   2.0240955312679931159E+02, 
   2.4736979003315290057E+02, 
   1.0742543875702278326E+02, 
   1.7463965060678569906E+01, 
   8.8427520398873480342E-01 };
  double sgn;
  double three = 3.0;
  double upper;
  double value;
  double w;
  double x;
  double x01 = 187.0;
  double x01d = 128.0;
  double x02 = 6.9464496836234126266E-04;
  double xinf = 1.70E+38;
  double xlarge = 2.04E+15;
  double xmax1 = 3.60E+16;
  double xmin1 = 5.89E-39;
  double xsmall = 2.05E-09;
  double z;
  double zero = 0.0;

  x = xx;
  w = fabs ( x );
  aug = zero;
//
//  Check for valid arguments, then branch to appropriate algorithm.
//
  if ( xmax1 <= - x || w < xmin1 )
  {
    if ( zero < x )
    {
      value = - xinf;
    }
    else
    {
      value = xinf;
    }
    return value;
  }

  if ( x < half )
  {
//
//  X < 0.5, use reflection formula: psi(1-x) = psi(x) + pi * cot(pi*x)
//  Use 1/X for PI*COTAN(PI*X)  when  XMIN1 < |X| <= XSMALL.
//
    if ( w <= xsmall )
    {
      aug = - one / x;
    }
//
//  Argument reduction for cotangent.
//
    else
    {
      if ( x < zero )
      {
        sgn = piov4;
      }
      else
      {
        sgn = - piov4;
      }

      w = w - ( double ) ( ( int ) ( w ) );
      nq = int ( w * four );
      w = four * ( w - ( double ) ( nq ) * fourth );
//
//  W is now related to the fractional part of 4.0 * X.
//  Adjust argument to correspond to values in the first
//  quadrant and determine the sign.
//
      n = nq / 2;

      if ( n + n != nq )
      {
        w = one - w;
      }

      z = piov4 * w;

      if ( ( n % 2 ) != 0 )
      {
        sgn = - sgn;
      }
//
//  Determine the final value for  -pi * cotan(pi*x).
//
      n = ( nq + 1 ) / 2;
      if ( ( n % 2 ) == 0 )
      {
//
//  Check for singularity.
//
        if ( z == zero )
        {
          if ( zero < x )
          {
            value = -xinf;
          }
          else
          {
            value = xinf;
          }
          return value;
        }
        aug = sgn * ( four / tan ( z ) );
      }
      else
      {
        aug = sgn * ( four * tan ( z ) );
      }
    }
    x = one - x;
  }
//
//  0.5 <= X <= 3.0.
//
  if ( x <= three )
  {
    den = x;
    upper = p1[0] * x;
    for ( i = 1; i <= 7; i++ )
    {
      den = ( den + q1[i-1] ) * x;
      upper = ( upper + p1[i]) * x;
    }
    den = ( upper + p1[8] ) / ( den + q1[7] );
    x = ( x - x01 / x01d ) - x02;
    value = den * x + aug;
    return value;
  }
//
//  3.0 < X.
//
  if ( x < xlarge )
  {
    w = one / ( x * x );
    den = w;
    upper = p2[0] * w;
    for ( i = 1; i <= 5; i++ )
    {
      den = ( den + q2[i-1] ) * w;
      upper = ( upper + p2[i] ) * w;
    }
    aug = ( upper + p2[6] ) / ( den + q2[5] ) - half / x + aug;
  }

  value = aug + log ( x );

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

double r8_sign ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    R8_SIGN returns the sign of an R8.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    18 October 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the number whose sign is desired.
//
//    Output, double R8_SIGN, the sign of X.
//
{
  double value;

  if ( x < 0.0 )
  {
    value = -1.0;
  } 
  else
  {
    value = 1.0;
  }
  return value;
}
//****************************************************************************80

double r8_uniform_ab ( double a, double b, int &seed )

//****************************************************************************80
//
//  Purpose:
//
//    R8_UNIFORM_AB returns a scaled pseudorandom R8.
//
//  Discussion:
//
//    The pseudorandom number should be uniformly distributed
//    between A and B.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    09 April 2012
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double A, B, the limits of the interval.
//
//    Input/output, int &SEED, the "seed" value, which should NOT be 0.
//    On output, SEED has been updated.
//
//    Output, double R8_UNIFORM_AB, a number strictly between A and B.
//
{
  const int i4_huge = 2147483647;
  int k;
  double value;

  if ( seed == 0 )
  {
    cerr << "\n";
    cerr << "R8_UNIFORM_AB - Fatal error!\n";
    cerr << "  Input value of SEED = 0.\n";
    exit ( 1 );
  }

  k = seed / 127773;

  seed = 16807 * ( seed - k * 127773 ) - k * 2836;

  if ( seed < 0 )
  {
    seed = seed + i4_huge;
  }

  value = ( double ) ( seed ) * 4.656612875E-10;

  value = a + ( b - a ) * value;

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

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

//****************************************************************************80
//
//  Purpose:
//
//    R8VEC_PRINT prints an R8VEC.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    16 August 2004
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N, the number of components of the vector.
//
//    double A[N], the vector to be printed.
//
//    string TITLE, a title.
//
{
  int i;

  cout << "\n";
  cout << title << "\n";
  cout << "\n";
  for ( i = 0; i < n; i++ )
  {
    cout << "  " << setw(8)  << i
         << ": " << setw(14) << a[i]  << "\n";
  }

  return;
}