subroutine gegenbauer_alpha_check ( alpha, check )

c*********************************************************************72
c
cc gegenbauer_alpha_check() checks the value of ALPHA.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    23 November 2015
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double precision ALPHA, a parameter which is part of the 
c    definition of the Gegenbauer polynomials.  It must be greater than -0.5.
c
c  Output:
c
c    logical CHECK.
c    TRUE, ALPHA is acceptable.
c    FALSE, ALPHA is not acceptable. 
c
      implicit none

      double precision alpha
      logical check
      logical squawk

      squawk = .false.

      if ( - 0.5D+00 < alpha ) then

        check = .true.

      else

        check = .false.

        if ( squawk ) then
          write ( *, '(a)' ) ''
          write ( *, '(a)' ) 
     &      'gegenbauer_polynomial_value(): Fatal error!';
          write ( *, '(a)' ) '  Illegal value of ALPHA.'
          write ( *, '(a,g14.6)' ) '  ALPHA = ', alpha 
          write ( *, '(a)' ) '  but ALPHA must be greater than -0.5.'
        end if

      end if

      return
      end
      subroutine gegenbauer_ek_compute ( n, alpha, x, w )

c*********************************************************************72
c
cc gegenbauer_ek_compute() computes a Gauss-Gegenbauer quadrature rule.
c
c  Discussion:
c
c    The integral:
c
c      integral ( -1 <= x <= 1 ) (1-x^2)^(alpha-1/2) * f(x) dx
c
c    The quadrature rule:
c
c      sum ( 1 <= i <= n ) w(i) * f ( x(i) )
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    23 November 2015
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Sylvan Elhay, Jaroslav Kautsky,
c    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of 
c    Interpolatory Quadrature,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 4, December 1987, pages 399-415.
c
c  Input:
c
c    integer N, the order.
c
c    double precision ALPHA, the exponent of (1-X*X) in the weight.
c    -1.0 < ALPHA.
c
c  Output:
c
c    double precision X(N), the abscissas.
c
c    double precision W(N), the weights.
c
      implicit none

      integer n

      double precision abi
      double precision alpha
      double precision bj(n)
      integer i
      double precision i_r8
      double precision w(n)
      double precision x(n)
      double precision zemu
c
c  Check N.
c
      if ( n < 1 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'gegenbauer_ek_compute(): Fatal error!'
        write ( *, '(a)' ) '  1 <= N is required.'
        stop 1
      end if
c
c  Check ALPHA.
c
      if ( alpha <= -1.0D+00 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'gegenbauer_ek_compute(): Fatal error!'
        write ( *, '(a)' ) '  -1.0 < ALPHA is required.'
        stop 1
      end if
c
c  Define the zero-th moment.
c
      zemu = 2.0 ** ( 2.0 * alpha + 1.0 ) 
     &  * gamma ( alpha + 1.0 ) 
     &  * gamma ( alpha + 1.0 ) 
     &  / gamma ( 2.0 * alpha + 2.0 )
c
c  Define the Jacobi matrix.
c
      x(1:n) = 0.0D+00

      bj(1) = 4.0D+00 * ( alpha + 1.0D+00 ) ** 2 
     &  / ( ( 2.0D+00 * alpha + 3.0D+00 ) 
     &  * ( 2.0D+00 * alpha + 2.0D+00 ) ** 2 )

      do i = 2, n
        i_r8 = dble ( i )
        abi = 2.0D+00 * ( alpha + i_r8 );
        bj(i) = 4.0D+00 * i_r8 * ( alpha + i_r8 ) ** 2 
     &    * ( 2.0D+00 * alpha + i_r8 ) 
     &    / ( ( abi - 1.0D+00 ) * ( abi + 1.0D+00 ) * abi * abi )
      end do

      bj(1:n) = sqrt ( bj(1:n) )

      w(1) = sqrt ( zemu )
      w(2:n) = 0.0D+00
c
c  Diagonalize the Jacobi matrix.
c
      call imtqlx ( n, x, bj, w )

      w(1:n) = w(1:n) ** 2

      return
      end
      subroutine gegenbauer_integral ( expon, alpha, value )

c*********************************************************************72
c
cc gegenbauer_integral(): integral of a monomial with Gegenbauer weight.
c
c  Discussion:
c
c    The integral:
c
c      integral ( -1 <= x <= +1 ) x^expon (1-x*x)^(alpha-1/2) dx
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    25 February 2008
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer EXPON, the exponent.
c
c    double precision ALPHA, part of the exponent of (1-X^2).
c    -0.5 < ALPHA.
c
c  Output:
c
c    double precision VALUE, the value of the integral.
c
      implicit none

      double precision alpha
      double precision arg1
      double precision arg2
      double precision arg3
      double precision arg4
      double precision c
      integer expon
      double precision value
      double precision value1

      if ( mod ( expon, 2 ) == 1 ) then
        value = 0.0D+00
        return
      end if

      c = dble ( expon )

      arg1 = - alpha
      arg2 =   1.0D+00 + c
      arg3 =   2.0D+00 + alpha + c
      arg4 = - 1.0D+00

      call r8_hyper_2f1 ( arg1, arg2, arg3, arg4, value1 )

      value = 2.0D+00 * gamma ( 1.0D+00 + c ) 
     &  * gamma ( 1.0D+00 + alpha ) 
     &  * value1 / gamma ( 2.0D+00 + alpha + c )

      return
      end
      subroutine gegenbauer_polynomial_value ( m, n, alpha, x, c )

c*********************************************************************72
c
cc gegenbauer_polynomial_value() computes the Gegenbauer polynomials C(I,ALPHA,X).
c
c  Discussion:
c
c    The Gegenbauer polynomial can be evaluated in Mathematica with
c    the command 
c
c      GegenbauerC[n,m,x]
c
c    ALPHA must be greater than -0.5.
c
c    If ALPHA = 1, the Gegenbauer polynomials reduce to the Chebyshev
c    polynomials of the second kind.
c
c  Differential equation:
c
c    (1-X*X) Y'' - (2 ALPHA + 1) X Y' + N (N + 2 ALPHA) Y = 0
c
c  Recursion:
c
c    C(0,ALPHA,X) = 1,
c    C(1,ALPHA,X) = 2*ALPHA*X
c    C(N,ALPHA,X) = ( (2*N-2+2*ALPHA) * X * C(N-1,ALPHA,X) 
c                   + ( -N+2-2*ALPHA)     * C(N-2,ALPHA,X) ) / N
c
c  Norm:
c
c    Integral ( -1 <= X <= 1 ) 
c      ( 1 - X^2 )^( ALPHA - 0.5 ) * C(N,ALPHA,X)^2 dX
c
c    = PI * 2^( 1 - 2 * ALPHA ) * Gamma ( N + 2 * ALPHA ) 
c      / ( Nc * ( N + ALPHA ) * ( Gamma ( ALPHA ) )^2 )
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    23 November 2015
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Stephen Wolfram,
c    The Mathematica Book,
c    Fourth Edition,
c    Cambridge University Press, 1999,
c    ISBN: 0-521-64314-7,
c    LC: QA76.95.W65.
c
c  Input:
c
c    integer M, the highest order polynomial to compute.
c    Note that polynomials 0 through M will be computed.
c
c    integer N, the number of evaluation points.
c
c    double precision ALPHA, a parameter which is part of the 
c    definition of the Gegenbauer polynomials.  It must be greater than -0.5.
c
c    double precision X(N), the evalulation points.
c
c  Output:
c
c    double precision C(0:M,N), the values of the first N+1 Gegenbauer
c    polynomials at the point X.  
c
      implicit none

      integer m
      integer n

      double precision alpha
      double precision c(0:m,n)
      integer i
      double precision x(n)

      if ( alpha <= -0.5D+00 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'gegenbauer_polynomial_value(): Fatal errorc'
        write ( *, '(a,g14.6)' ) '  Illegal value of ALPHA = ', alpha
        write ( *, '(a)' ) '  but ALPHA must be greater than -0.5.'
        stop 1
      end if

      if ( m < 0 ) then
        return
      end if

      if ( n < 1 ) then
        return
      end if

      c(0,1:n) = 1.0D+00

      if ( m < 1 ) then
        return
      end if

      c(1,1:n) = 2.0D+00 * alpha * x(1:n)

      do i = 2, m
        c(i,1:n) = 
     &    ( ( dble ( 2 * i - 2 ) + 2.0D+00 * alpha ) 
     &    * x(1:n) * c(i-1,1:n)   
     &    + ( dble (   - i + 2 ) - 2.0D+00 * alpha )
     &             * c(i-2,1:n) ) 
     &    /   dble (     i      )
      end do
     
      return
      end
      subroutine gegenbauer_polynomial_values ( n_data, n, a, x, fx )

c*********************************************************************72
c
cc gegenbauer_polynomial_values() returns some values of the Gegenbauer polynomials.
c
c  Discussion:
c
c    The Gegenbauer polynomials are also known as the "spherical
c    polynomials" or "ultraspherical polynomials".
c
c    In Mathematica, the function can be evaluated by:
c
c      GegenbauerC[n,m,x]
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    24 March 2007
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Milton Abramowitz, Irene Stegun,
c    Handbook of Mathematical Functions,
c    National Bureau of Standards, 1964,
c    ISBN: 0-486-61272-4,
c    LC: QA47.A34.
c
c    Stephen Wolfram,
c    The Mathematica Book,
c    Fourth Edition,
c    Cambridge University Press, 1999,
c    ISBN: 0-521-64314-7,
c    LC: QA76.95.W65.
c
c  Parameters:
c
c    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
c    first call.  On each call, the routine increments N_DATA by 1, and
c    returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    Output, integer N, the order parameter of the function.
c
c    Output, double precision A, the real parameter of the function.
c
c    Output, double precision X, the argument of the function.
c
c    Output, double precision FX, the value of the function.
c
      implicit none

      integer n_max
      parameter ( n_max = 38 )

      double precision a
      double precision a_vec(n_max)
      double precision fx
      double precision fx_vec(n_max)
      integer n
      integer n_data
      integer n_vec(n_max)
      double precision x
      double precision x_vec(n_max)

      save a_vec
      save fx_vec
      save n_vec
      save x_vec

      data a_vec /
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.5D+00,
     &   0.0D+00,
     &   1.0D+00,
     &   2.0D+00,
     &   3.0D+00,
     &   4.0D+00,
     &   5.0D+00,
     &   6.0D+00,
     &   7.0D+00,
     &   8.0D+00,
     &   9.0D+00,
     &  10.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00,
     &   3.0D+00 /
      data fx_vec /
     &    1.0000000000D+00,
     &    0.2000000000D+00,
     &   -0.4400000000D+00,
     &   -0.2800000000D+00,
     &    0.2320000000D+00,
     &    0.3075200000D+00,
     &   -0.0805760000D+00,
     &   -0.2935168000D+00,
     &   -0.0395648000D+00,
     &    0.2459712000D+00,
     &    0.1290720256D+00,
     &    0.0000000000D+00,
     &   -0.3600000000D+00,
     &   -0.0800000000D+00,
     &    0.8400000000D+00,
     &    2.4000000000D+00,
     &    4.6000000000D+00,
     &    7.4400000000D+00,
     &   10.9200000000D+00,
     &   15.0400000000D+00,
     &   19.8000000000D+00,
     &   25.2000000000D+00,
     &   -9.0000000000D+00,
     &   -0.1612800000D+00,
     &   -6.6729600000D+00,
     &   -8.3750400000D+00,
     &   -5.5267200000D+00,
     &    0.0000000000D+00,
     &    5.5267200000D+00,
     &    8.3750400000D+00,
     &    6.6729600000D+00,
     &    0.1612800000D+00,
     &   -9.0000000000D+00,
     &  -15.4252800000D+00,
     &   -9.6969600000D+00,
     &   22.4409600000D+00,
     &  100.8892800000D+00,
     &  252.0000000000D+00 /
      data n_vec /
     &   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 /
      data x_vec /
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.20D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &   0.40D+00,
     &  -0.50D+00,
     &  -0.40D+00,
     &  -0.30D+00,
     &  -0.20D+00,
     &  -0.10D+00,
     &   0.00D+00,
     &   0.10D+00,
     &   0.20D+00,
     &   0.30D+00,
     &   0.40D+00,
     &   0.50D+00,
     &   0.60D+00,
     &   0.70D+00,
     &   0.80D+00,
     &   0.90D+00,
     &   1.00D+00 /

      if ( n_data .lt. 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max .lt. n_data ) then
        n_data = 0
        n = 0
        a = 0.0D+00
        x = 0.0D+00
        fx = 0.0D+00
      else
        n = n_vec(n_data)
        a = a_vec(n_data)
        x = x_vec(n_data)
        fx = fx_vec(n_data)
      end if

      return
      end
      subroutine gegenbauer_ss_compute ( n, alpha, x, w )

c*********************************************************************72
c
cc gegenbauer_ss_compute() computes a Gauss-Gegenbauer quadrature rule.
c
c  Discussion:
c
c    The integral:
c
c      integral ( -1 <= x <= 1 ) (1-x^2)^(alpha-1/2) * f(x) dx
c
c    The quadrature rule:
c
c      sum ( 1 <= i <= n ) w(i) * f ( x(i) )
c
c    Thanks to Janiki Raman for pointing out a problem in an earlier
c    version of the code that occurred when ALPHA was -0.5.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    24 June 2008
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Arthur Stroud, Don Secrest,
c    Gaussian Quadrature Formulas,
c    Prentice Hall, 1966,
c    LC: QA299.4G3S7.
c
c  Input:
c
c    integer N, the order.
c
c    double precision ALPHA, the exponent of (1-X*X) in the weight.
c    -1.0 < ALPHA.
c
c  Output:
c
c    double precision X(N), the abscissas.
c
c    double precision W(N), the weights.
c
      implicit none

      integer n

      double precision alpha
      double precision an
      double precision c(n)
      double precision cc
      double precision delta
      double precision dp2
      integer i
      double precision p1
      double precision r1
      double precision r2
      double precision r3
      double precision w(n)
      double precision x(n)
      double precision xval
c
c  Check N.
c
      if ( n < 1 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'gegenbauer_ss_compute(): Fatal errorc'
        write ( *, '(a)' ) '  1 <= N is required.'
        stop 1
      end if
c
c  Check ALPHA.
c
      if ( alpha <= -1.0D+00 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'gegenbauer_ss_compute(): Fatal errorc'
        write ( *, '(a)' ) '  -1.0 < ALPHA is required.'
        stop 1
      end if
c
c  Set the recursion coefficients.
c
      c(1) = 0.0D+00

      if ( 2 <= n ) then
        c(2) = 1.0D+00 / ( 2.0D+00 * alpha + 3.0D+00 )
      end if

      do i = 3, n

        c(i) = dble ( i - 1 ) 
     &        * ( alpha + alpha + dble ( i - 1     ) ) / 
     &        ( ( alpha + alpha + dble ( 2 * i - 1 ) ) 
     &        * ( alpha + alpha + dble ( 2 * i - 3 ) ) )

      end do

      delta = gamma ( alpha         + 1.0D+00 ) 
     &      * gamma (         alpha + 1.0D+00 ) 
     &      / gamma ( alpha + alpha + 2.0D+00 )

      cc = delta * 2.0D+00 ** ( 2.0D+00 * alpha + 1.0D+00 ) 
     &  * product ( c(2:n) )

      do i = 1, n

        if ( i == 1 ) then

          an = alpha / dble ( n )

          r1 = ( 1.0D+00 + alpha ) 
     &      * ( 2.78D+00 / ( 4.0D+00 + dble ( n * n ) ) 
     &      + 0.768D+00 * an / dble ( n ) )

          r2 = 1.0D+00 + 2.44D+00 * an + 1.282D+00 * an * an

          xval = ( r2 - r1 ) / r2

        else if ( i == 2 ) then

          r1 = ( 4.1D+00 + alpha ) / 
     &      ( ( 1.0D+00 + alpha ) * ( 1.0D+00 + 0.156D+00 * alpha ) )

          r2 = 1.0D+00 + 0.06D+00 * ( dble ( n ) - 8.0D+00 ) * 
     &      ( 1.0D+00 + 0.12D+00 * alpha ) / dble ( n )

          r3 = 1.0D+00 + 0.012D+00 * alpha * 
     &      ( 1.0D+00 + 0.25D+00 * abs ( alpha ) ) / dble ( n )

          xval = xval - r1 * r2 * r3 * ( 1.0D+00 - xval )

        else if ( i == 3 ) then

          r1 = ( 1.67D+00 + 0.28D+00 * alpha ) 
     &      / ( 1.0D+00 + 0.37D+00 * alpha )

          r2 = 1.0D+00 + 0.22D+00 * ( dble ( n ) - 8.0D+00 ) 
     &      / dble ( n )

          r3 = 1.0D+00 + 8.0D+00 * alpha / 
     &      ( ( 6.28D+00 + alpha ) * dble ( n * n ) )

          xval = xval - r1 * r2 * r3 * ( x(1) - xval )

        else if ( i < n - 1 ) then

          xval = 3.0D+00 * x(i-1) - 3.0D+00 * x(i-2) + x(i-3)

        else if ( i == n - 1 ) then

          r1 = ( 1.0D+00 + 0.235D+00 * alpha ) 
     &      / ( 0.766D+00 + 0.119D+00 * alpha )

          r2 = 1.0D+00 / ( 1.0D+00 + 0.639D+00 
     &      * ( dble ( n ) - 4.0D+00 ) 
     &      / ( 1.0D+00 + 0.71D+00 * ( dble ( n ) - 4.0D+00 ) ) )

          r3 = 1.0D+00 / ( 1.0D+00 + 20.0D+00 * alpha 
     &      / ( ( 7.5D+00 + alpha ) * dble ( n * n ) ) )

          xval = xval + r1 * r2 * r3 * ( xval - x(i-2) )

        else if ( i == n ) then

          r1 = ( 1.0D+00 + 0.37D+00 * alpha ) 
     &      / ( 1.67D+00 + 0.28D+00 * alpha )

          r2 = 1.0D+00 / 
     &      ( 1.0D+00 + 0.22D+00 * ( dble ( n ) - 8.0D+00 ) 
     &      / dble ( n ) )

          r3 = 1.0D+00 / ( 1.0D+00 + 8.0D+00 * alpha / 
     &      ( ( 6.28D+00 + alpha ) * dble ( n * n ) ) )

          xval = xval + r1 * r2 * r3 * ( xval - x(i-2) )

        end if

        call gegenbauer_ss_root ( xval, n, dp2, p1, c )

        x(i) = xval
        w(i) = cc / ( dp2 * p1 )

      end do
c
c  Reverse the data.
c
      x(1:n) = x(n:1:-1)
      w(1:n) = w(n:1:-1)

      return
      end
      subroutine gegenbauer_ss_recur ( p2, dp2, p1, x, n, c )

c*********************************************************************72
c
cc gegenbauer_ss_recur(): value and derivative of a Gegenbauer polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    05 April 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Arthur Stroud, Don Secrest,
c    Gaussian Quadrature Formulas,
c    Prentice Hall, 1966,
c    LC: QA299.4G3S7.
c
c  Input:
c
c    double precision X, the point at which polynomials are evaluated.
c
c    integer N, the order of the polynomial.
c
c    double precision C(N), the recursion coefficients.
c
c  Output:
c
c    double precision P2, the value of J(N)(X).
c
c    double precision DP2, the value of J'(N)(X).
c
c    double precision P1, the value of J(N-1)(X).
c
      implicit none

      integer n

      double precision c(n)
      double precision dp0
      double precision dp1
      double precision dp2
      integer i
      double precision p0
      double precision p1
      double precision p2
      double precision x

      p1 = 1.0D+00
      dp1 = 0.0D+00

      p2 = x
      dp2 = 1.0D+00

      do i = 2, n

        p0 = p1
        dp0 = dp1

        p1 = p2
        dp1 = dp2

        p2 = x * p1 - c(i) * p0
        dp2 = x * dp1 + p1 - c(i) * dp0

      end do

      return
      end
      subroutine gegenbauer_ss_root ( x, n, dp2, p1, c )

c*********************************************************************72
c
cc gegenbauer_ss_root() improves an approximate root of a Gegenbauer polynomial.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    05 April 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Arthur Stroud, Don Secrest,
c    Gaussian Quadrature Formulas,
c    Prentice Hall, 1966,
c    LC: QA299.4G3S7.
c
c  Input:
c
c    double precision X, the approximate root.
c
c    integer N, the order of the polynomial.
c
c    double precision C(N), the recursion coefficients.
c
c  Output:
c
c    double precision X, the improved root.
c
c    double precision DP2, the value of J'(N)(X).
c
c    double precision P1, the value of J(N-1)(X).
c
      implicit none

      integer n

      double precision c(n)
      double precision d
      double precision dp2
      double precision eps
      double precision p1
      double precision p2
      integer step
      integer, parameter :: step_max = 10
      double precision x

      eps = epsilon ( eps )

      do step = 1, step_max

        call gegenbauer_ss_recur ( p2, dp2, p1, x, n, c )

        d = p2 / dp2
        x = x - d

        if ( abs ( d ) <= eps * ( abs ( x ) + 1.0D+00 ) ) then
          return
        end if

      end do

      return
      end
      subroutine gegenbauer_to_monomial_matrix ( n, alpha, A )

c*********************************************************************72
c
cc gegenbauer_to_monomial_matrix(): Gegenbauer to monomial conversion matrix.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    05 April 2024
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N: the order of A.
c
c    real alpha: the parameter.
c
c  Output:
c
c    double precision A(N,N): the matrix.
c
      implicit none

      integer n

      double precision A(n,n)
      double precision alpha
      double precision c1
      double precision c2
      integer i
      integer j
      integer nn

      if ( n <= 0 ) then
        return
      end if

      do j = 1, n
        do i = 1, n
          A(i,j) = 0.0D+00
        end do
      end do

      A(1,1) = 1.0D+00

      if ( n == 1 ) then
        return
      end if

      A(2,2) = 2.0D+00 * alpha
c
c  Perform convex sum.
c  Translating "(n+1) C(n+1) = 2 (n+alpha) x C(n) - ( n + 2 alpha - 1 ) C(n-1)"
c  drove me nuts, between indexing at 1 rather than 0, and dealing with
c  the interpretation of "n+1", because I now face the rare "off by 2" errorc
c
      do j = 3, n
        nn = j - 2
        c1 = ( 2.0 * nn + 2.0 * alpha       ) / dble ( nn + 1 )
        c2 = (     - nn - 2.0 * alpha + 1.0 ) / dble ( nn + 1 )
        A(2:j,j)   =              c1 * A(1:j-1,j-1)
        A(1:j-2,j) = A(1:j-2,j) + c2 * A(1:j-2,j-2)
      end do

      return
      end
      subroutine hyper_2f1_values ( n_data, a, b, c, x, fx )

c*********************************************************************72
c
cc hyper_2f1_values() returns some values of the hypergeometric function 2F1.
c
c  Discussion:
c
c    In Mathematica, the function can be evaluated by:
c
c      fx = Hypergeometric2F1 [ a, b, c, x ]
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    09 September 2007
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Milton Abramowitz, Irene Stegun,
c    Handbook of Mathematical Functions,
c    National Bureau of Standards, 1964,
c    ISBN: 0-486-61272-4,
c    LC: QA47.A34.
c
c    Stephen Wolfram,
c    The Mathematica Book,
c    Fourth Edition,
c    Cambridge University Press, 1999,
c    ISBN: 0-521-64314-7,
c    LC: QA76.95.W65.
c
c    Shanjie Zhang, Jianming Jin,
c    Computation of Special Functions,
c    Wiley, 1996,
c    ISBN: 0-471-11963-6,
c    LC: QA351.C45
c
c    Daniel Zwillinger, editor,
c    CRC Standard Mathematical Tables and Formulae,
c    30th Edition,
c    CRC Press, 1996,
c    ISBN: 0-8493-2479-3,
c    LC: QA47.M315.
c
c  Input:
c
c    integer N_DATA.  The user sets N_DATA to 0 before the first call.
c
c  Output:
c
c    integer N_DATA.  The routine increments N_DATA by 1,
c    and returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    double precision A, B, C: the parameters.
c
c    double precision X: the argument.
c
c    double precision FX, the value of the function.
c
      implicit none

      integer n_max
      parameter ( n_max = 24 )

      double precision a
      double precision a_vec(n_max)
      double precision b
      double precision b_vec(n_max)
      double precision c
      double precision c_vec(n_max)
      double precision fx
      double precision fx_vec(n_max)
      integer n_data
      double precision x
      double precision x_vec(n_max)

      save a_vec
      save b_vec
      save c_vec
      save fx_vec
      save x_vec

      data a_vec /
     &   -2.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    2.5D+00,
     &   -2.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    2.5D+00,
     &   -2.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    2.5D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00 /
      data b_vec /
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    3.3D+00,
     &    1.1D+00,
     &    1.1D+00,
     &    3.3D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00 /
      data c_vec /
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &    6.7D+00,
     &   -5.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    4.5D+00,
     &   -5.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    4.5D+00,
     &   -5.5D+00,
     &   -0.5D+00,
     &    0.5D+00,
     &    4.5D+00 /
      data fx_vec /
     &    0.72356129348997784913D+00,
     &    0.97911109345277961340D+00,
     &    1.0216578140088564160D+00,
     &    1.4051563200112126405D+00,
     &    0.46961431639821611095D+00,
     &    0.95296194977446325454D+00,
     &    1.0512814213947987916D+00,
     &    2.3999062904777858999D+00,
     &    0.29106095928414718320D+00,
     &    0.92536967910373175753D+00,
     &    1.0865504094806997287D+00,
     &    5.7381565526189046578D+00,
     &    15090.669748704606754D+00,
     &   -104.31170067364349677D+00,
     &    21.175050707768812938D+00,
     &    4.1946915819031922850D+00,
     &    1.0170777974048815592D+10,
     &   -24708.635322489155868D+00,
     &    1372.2304548384989560D+00,
     &    58.092728706394652211D+00,
     &    5.8682087615124176162D+18,
     &   -4.4635010147295996680D+08,
     &    5.3835057561295731310D+06,
     &    20396.913776019659426D+00 /
      data x_vec /
     &    0.25D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.25D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.55D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.85D+00,
     &    0.85D+00 /

      if ( n_data .lt. 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max .lt. n_data ) then
        n_data = 0
        a = 0.0D+00
        b = 0.0D+00
        c = 0.0D+00
        x = 0.0D+00
        fx = 0.0D+00
      else
        a = a_vec(n_data)
        b = b_vec(n_data)
        c = c_vec(n_data)
        x = x_vec(n_data)
        fx = fx_vec(n_data)
      end if

      return
      end
      subroutine imtqlx ( n, d, e, z )

c*********************************************************************72
c
cc imtqlx() diagonalizes a symmetric tridiagonal matrix.
c
c  Discussion:
c
c    This routine is a slightly modified version of the EISPACK routine to
c    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
c
c    The authors thank the authors of EISPACK for permission to use this
c    routine.
c
c    It has been modified to produce the product Q' * Z, where Z is an input
c    vector and Q is the orthogonal matrix diagonalizing the input matrix.
c    The changes consist (essentially) of applying the orthogonal
c    transformations directly to Z as they are generated.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    27 December 2009
c
c  Author:
c
c    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
c    This version by John Burkardt.
c
c  Reference:
c
c    Sylvan Elhay, Jaroslav Kautsky,
c    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
c    Interpolatory Quadrature,
c    ACM Transactions on Mathematical Software,
c    Volume 13, Number 4, December 1987, pages 399-415.
c
c    Roger Martin, James Wilkinson,
c    The Implicit QL Algorithm,
c    Numerische Mathematik,
c    Volume 12, Number 5, December 1968, pages 377-383.
c
c  Parameters:
c
c    Input, integer N, the order of the matrix.
c
c    Input/output, double precision D(N), the diagonal entries of the matrix.
c    On output, the information in D has been overwritten.
c
c    Input/output, double precision E(N), the subdiagonal entries of the
c    matrix, in entries E(1) through E(N-1).  On output, the information in
c    E has been overwritten.
c
c    Input/output, double precision Z(N).  On input, a vector.  On output,
c    the value of Q' * Z, where Q is the matrix that diagonalizes the
c    input symmetric tridiagonal matrix.
c
      implicit none

      integer n

      double precision b
      double precision c
      double precision d(n)
      double precision e(n)
      double precision f
      double precision g
      integer i
      integer ii
      integer, parameter :: itn = 30
      integer j
      integer k
      integer l
      integer m
      integer mml
      double precision p
      double precision prec
      double precision r
      double precision s
      double precision z(n)

      prec = epsilon ( prec )

      if ( n == 1 ) then
        return
      end if

      e(n) = 0.0D+00

      do l = 1, n

        j = 0

        do

          do m = l, n

            if ( m == n ) then
              exit
            end if

            if ( abs ( e(m) ) <= 
     &        prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) ) then
              exit
            end if

          end do

          p = d(l)

          if ( m == l ) then
            exit
          end if

          if ( itn <= j ) then
            write ( *, '(a)' ) ' '
            write ( *, '(a)' ) 'imtqlx(): Fatal error.'
            write ( *, '(a)' ) '  Iteration limit exceeded.'
            write ( *, '(a,i8)' ) '  J = ', j
            write ( *, '(a,i8)' ) '  L = ', l
            write ( *, '(a,i8)' ) '  M = ', m
            write ( *, '(a,i8)' ) '  N = ', n
            stop 1
          end if

          j = j + 1
          g = ( d(l+1) - p ) / ( 2.0D+00 * e(l) )
          r =  sqrt ( g * g + 1.0D+00 )
          g = d(m) - p + e(l) / ( g + sign ( r, g ) )
          s = 1.0D+00
          c = 1.0D+00
          p = 0.0D+00
          mml = m - l

          do ii = 1, mml

            i = m - ii
            f = s * e(i)
            b = c * e(i)

            if ( abs ( g ) <= abs ( f ) ) then
              c = g / f
              r =  sqrt ( c * c + 1.0D+00 )
              e(i+1) = f * r
              s = 1.0D+00 / r
              c = c * s
            else
              s = f / g
              r =  sqrt ( s * s + 1.0D+00 )
              e(i+1) = g * r
              c = 1.0D+00 / r
              s = s * c
            end if

            g = d(i+1) - p
            r = ( d(i) - g ) * s + 2.0D+00 * c * b
            p = s * r
            d(i+1) = g + p
            g = c * r - b
            f = z(i+1)
            z(i+1) = s * z(i) + c * f
            z(i) = c * z(i) - s * f

          end do

          d(l) = d(l) - p
          e(l) = g
          e(m) = 0.0D+00

        end do

      end do
c
c  Sorting.
c
      do ii = 2, n

        i = ii - 1
        k = i
        p = d(i)

        do j = ii, n
          if ( d(j) < p ) then
            k = j
            p = d(j)
          end if
        end do

        if ( k /= i ) then
          d(k) = d(i)
          d(i) = p
          p = z(i)
          z(i) = z(k)
          z(k) = p
        end if

      end do

      return
      end
      subroutine monomial_to_gegenbauer_matrix ( n, alpha, A )

c*********************************************************************72
c
cc monomial_to_gegenbauer_matrix(): monomial to Gegenbauer conversion matrix.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    24 April 2024
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer n: the order of A.
c
c    double precision alpha: the parameter.
c
c  Output:
c
c    double precision A[N,N]: the matrix.
c
      implicit none

      integer n

      double precision A(0:n-1,0:n-1)
      double precision alpha
      double precision bot
      integer i
      integer ilo
      integer j
      double precision r8_factorial
      double precision top

      A(0:n-1,0:n-1) = 0.0D+00

      if ( n <= 0 ) then
        return
      end if

      do j = 0, n - 1

        ilo = mod ( j, 2 )

        do i = ilo, j, 2

          top = ( i + alpha ) * r8_factorial ( j ) 
     &      * gamma ( alpha )

          bot = 2.0D+00 ** j * r8_factorial ( ( j - i ) / 2 ) 
     &      * gamma ( ( j + i ) / 2 + alpha + 1.0D+00 )

          A(i,j) = top / bot

        end do
      end do

      return
      end
      subroutine psi_values ( n_data, x, fx )

c*********************************************************************72
c
cc psi_values() returns some values of the Psi or Digamma function for testing.
c
c  Discussion:
c
c    PSI(X) = d LN ( GAMMA ( X ) ) / d X = GAMMA'(X) / GAMMA(X)
c
c    PSI(1) = - Euler's constant.
c
c    PSI(X+1) = PSI(X) + 1 / X.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    04 June 2013
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Milton Abramowitz, Irene Stegun,
c    Handbook of Mathematical Functions,
c    National Bureau of Standards, 1964,
c    ISBN: 0-486-61272-4,
c    LC: QA47.A34.
c
c  Input:
c
c    integer N_DATA.  The user sets N_DATA to 0 before the first call.  
c
c  Output:
c
c    integer N_DATA.  The routine increments N_DATA by 1, and
c    returns the corresponding data; when there is no more data, the
c    output value of N_DATA will be 0 again.
c
c    double precision X, the argument of the function.
c
c    double precision FX, the value of the function.
c
      implicit none

      integer n_max
      parameter ( n_max = 20 )

      double precision fx
      double precision fxvec ( n_max )
      integer n_data
      double precision x
      double precision xvec ( n_max )

      data fxvec /
     & -10.42375494041108D+00, 
     &  -5.289039896592188D+00, 
     &  -3.502524222200133D+00, 
     &  -2.561384544585116D+00, 
     &  -1.963510026021423D+00, 
     &  -1.540619213893190D+00, 
     &  -1.220023553697935D+00, 
     &  -0.9650085667061385D+00, 
     &  -0.7549269499470514D+00, 
     &  -0.5772156649015329D+00,
     &  -0.4237549404110768D+00,
     &  -0.2890398965921883D+00,
     &  -0.1691908888667997D+00,
     &  -0.6138454458511615D-01,
     &   0.3648997397857652D-01,
     &   0.1260474527734763D+00,
     &   0.2085478748734940D+00,
     &   0.2849914332938615D+00,
     &   0.3561841611640597D+00,
     &   0.4227843350984671D+00 /

      data xvec /
     &  0.1D+00,
     &  0.2D+00,
     &  0.3D+00,
     &  0.4D+00,
     &  0.5D+00,
     &  0.6D+00,
     &  0.7D+00,
     &  0.8D+00,
     &  0.9D+00,
     &  1.0D+00,
     &  1.1D+00,
     &  1.2D+00,
     &  1.3D+00,
     &  1.4D+00,
     &  1.5D+00,
     &  1.6D+00,
     &  1.7D+00,
     &  1.8D+00,
     &  1.9D+00,
     &  2.0D+00 /

      if ( n_data .lt. 0 ) then
        n_data = 0
      end if

      n_data = n_data + 1

      if ( n_max .lt. n_data ) then
        n_data = 0
        x = 0.0D+00
        fx = 0.0D+00
      else
        x = xvec(n_data)
        fx = fxvec(n_data)
      end if

      return
      end
      function r8_factorial ( n )

c*********************************************************************72
c
cc r8_factorial() computes the factorial of N.
c
c  Discussion:
c
c    factorial ( N ) = product ( 1 <= I <= N ) I
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    16 January 1999
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N: the argument of the factorial function.
c    If N is less than 1, the function value is returned as 1.
c
c  Output:
c
c    double precision R8_FACTORIAL: the factorial of N.
c
      implicit none

      double precision r8_factorial
      integer i
      integer n
      double precision value

      value = 1.0D+00

      do i = 1, n
        value = value * dble ( i )
      end do

      r8_factorial = value

      return
      end
      subroutine r8_hyper_2f1 ( a_input, b_input, c_input, x_input, hf )

c*********************************************************************72
c
cc r8_hyper_2f1() evaluates the hypergeometric function F(A,B,C,X).
c
c  Discussion:
c
c    A minor bug was corrected.  The HW variable, used in several places as
c    the "old" value of a quantity being iteratively improved, was not
c    being initialized.  JVB, 11 February 2008.
c
c    The original version of this program allowed the input arguments to
c    be modified, although they were restored to their input values before exit.
c    This is unacceptable if the input arguments are allowed to be constants.
c    The code has been modified so that the input arguments are never modified.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    11 October 2008
c
c  Author:
c
c    Original FORTRAN77 version by Shanjie Zhang, Jianming Jin.
c    This version by John Burkardt.
c
c    The original FORTRAN77 version of this routine is copyrighted by
c    Shanjie Zhang and Jianming Jin.  However, they give permission to
c    incorporate this routine into a user program provided that the copyright
c    is acknowledged.
c
c  Reference:
c
c    Shanjie Zhang, Jianming Jin,
c    Computation of Special Functions,
c    Wiley, 1996,
c    ISBN: 0-471-11963-6,
c    LC: QA351.C45
c
c  Input:
c
c    double precision A_INPUT, B_INPUT, C_INPUT, X_INPUT,
c    the arguments of the function.  The user is allowed to pass these
c    values as constants or variables.
c    C_INPUT must not be equal to a nonpositive integer.
c    X_INPUT < 1.
c
c  Output:
c
c    double precision HF, the value of the function.
c
      implicit none

      double precision a
      double precision a_input
      double precision a0
      double precision aa
      double precision b
      double precision b_input
      double precision bb
      double precision c
      double precision c_input
      double precision c0
      double precision c1
      double precision, parameter :: el = 0.5772156649015329D+00
      double precision eps
      double precision f0
      double precision f1
      double precision g0
      double precision g1
      double precision g2
      double precision g3
      double precision ga
      double precision gabc
      double precision gam
      double precision gb
      double precision gbm
      double precision gc
      double precision gca
      double precision gcab
      double precision gcb
      double precision gm
      double precision hf
      double precision hw
      integer j
      integer k
      logical l0
      logical l1
      logical l2
      logical l3
      logical l4
      logical l5
      integer m
      integer nm
      double precision pa
      double precision pb
      double precision r
      double precision r0
      double precision r1
      double precision, parameter :: r8_pi = 3.141592653589793D+00
      double precision r8_psi
      double precision rm
      double precision rp
      double precision sm
      double precision sp
      double precision sp0
      double precision x
      double precision x_input
      double precision x1
c
c  Immediately copy the input argumentsc
c
      a = a_input
      b = b_input
      c = c_input
      x = x_input

      l0 = ( c == aint ( c ) ) .and. ( c < 0.0D+00 )
      l1 = ( 1.0D+00 - x < 1.0D-15 ) .and. ( c - a - b <= 0.0D+00 )
      l2 = ( a == aint ( a ) ) .and. ( a < 0.0D+00 )
      l3 = ( b == aint ( b ) ) .and. ( b < 0.0D+00 )
      l4 = ( c - a == aint ( c - a ) ) .and. ( c - a <= 0.0D+00 )
      l5 = ( c - b == aint ( c - b ) ) .and. ( c - b <= 0.0D+00 )

      if ( l0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'r8_hyper_2f1(): Fatal error!'
        write ( *, '(a)' ) '  Integral C < 0.'
        write ( *, '(a)' ) '  The hypergeometric series is divergent.'
        stop 1
      end if

      if ( l1 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'r8_hyper_2f1(): Fatal error!'
        write ( *, '(a)' ) '  The hypergeometric series is divergent.'
        write ( *, '(a)' ) '  1 - X < 0, C - A - B < 0.'
        stop 1
      end if

      if ( 0.95D+00 < x ) then
        eps = 1.0D-08
      else
        eps = 1.0D-15
      end if

      if ( x == 0.0D+00 .or. a == 0.0D+00 .or. b == 0.0D+00 ) then

        hf = 1.0D+00
        return

      else if ( 1.0D+00 - x == eps .and. 0.0D+00 < c - a - b ) then

        gc = gamma ( c )
        gcab = gamma ( c - a - b )
        gca = gamma ( c - a )
        gcb = gamma ( c - b )
        hf = gc * gcab / ( gca * gcb )
        return

      else if ( 1.0D+00 + x <= eps .and. 
     &  abs ( c - a + b - 1.0D+00 ) <= eps ) then

        g0 = sqrt ( r8_pi ) * 2.0D+00 ** ( - a )
        g1 = gamma ( c )
        g2 = gamma ( 1.0D+00 + a / 2.0D+00 - b )
        g3 = gamma ( 0.5D+00 + 0.5D+00 * a )
        hf = g0 * g1 / ( g2 * g3 )
        return

      else if ( l2 .or. l3 ) then

        if ( l2 ) then
          nm = int ( abs ( a ) )
        end if

        if ( l3 ) then
          nm = int ( abs ( b ) )
        end if

        hf = 1.0D+00
        r = 1.0D+00

        do k = 1, nm
          r = r * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) 
     &      / ( k * ( c + k - 1.0D+00 ) ) * x
          hf = hf + r
        end do

        return

      else if ( l4 .or. l5 ) then

        if ( l4 ) then
          nm = int ( abs ( c - a ) )
        end if

        if ( l5 ) then
          nm = int ( abs ( c - b ) )
        end if

        hf = 1.0D+00
        r  = 1.0D+00
        do k = 1, nm
          r = r * ( c - a + k - 1.0D+00 ) * ( c - b + k - 1.0D+00 ) 
     &      / ( k * ( c + k - 1.0D+00 ) ) * x
          hf = hf + r
        end do
        hf = ( 1.0D+00 - x )**( c - a - b ) * hf
        return

      end if

      aa = a
      bb = b
      x1 = x

      if ( x < 0.0D+00 ) then
        x = x / ( x - 1.0D+00 )
        if ( a < c .and. b < a .and. 0.0D+00 < b ) then
          a = bb
          b = aa
        end if
        b = c - b
      end if

      if ( 0.75D+00 <= x ) then

        gm = 0.0D+00

        if ( abs ( c - a - b - aint ( c - a - b ) ) < 1.0D-15 ) then

          m = int ( c - a - b )
          ga = gamma ( a )
          gb = gamma ( b )
          gc = gamma ( c )
          gam = gamma ( a + m )
          gbm = gamma ( b + m )

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

          if ( m /= 0 ) then
            gm = 1.0D+00
          end if

          do j = 1, abs ( m ) - 1
            gm = gm * j
          end do

          rm = 1.0D+00
          do j = 1, abs ( m )
            rm = rm * j
          end do

          f0 = 1.0D+00
          r0 = 1.0D+00
          r1 = 1.0D+00
          sp0 = 0.0D+00
          sp = 0.0D+00

          if ( 0 <= m ) then

            c0 = gm * gc / ( gam * gbm )
            c1 = - gc * ( x - 1.0D+00 )**m / ( ga * gb * rm )

            do k = 1, m - 1
              r0 = r0 * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) 
     &          / ( k * ( k - m ) ) * ( 1.0D+00 - x )
              f0 = f0 + r0
            end do

            do k = 1, m
              sp0 = sp0 + 1.0D+00 / ( a + k - 1.0D+00 ) 
     &          + 1.0D+00 / ( b + k - 1.0D+00 ) - 1.0D+00 
     &          / dble ( k )
            end do

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

            do k = 1, 250

              sp = sp + ( 1.0D+00 - a ) / ( k * ( a + k - 1.0D+00 ) ) 
     &          + ( 1.0D+00 - b ) / ( k * ( b + k - 1.0D+00 ) )

              sm = 0.0D+00
              do j = 1, m
                sm = sm + ( 1.0D+00 - a ) 
     &            / ( ( j + k ) * ( a + j + k - 1.0D+00 ) ) 
     &            + 1.0D+00 / ( b + j + k - 1.0D+00 )
              end do

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

              r1 = r1 * ( a + m + k - 1.0D+00 ) 
     &          * ( b + m + k - 1.0D+00 ) 
     &          / ( k * ( m + k ) ) * ( 1.0D+00 - x )

              f1 = f1 + r1 * rp

              if ( abs ( f1 - hw ) < abs ( f1 ) * eps ) then
                exit
              end if

              hw = f1

            end do

            hf = f0 * c0 + f1 * c1

          else if ( m < 0 ) then

            m = - m
            c0 = gm * gc / ( ga * gb * ( 1.0D+00 - x )**m )
            c1 = - ( - 1 )**m * gc / ( gam * gbm * rm )

            do k = 1, m - 1
              r0 = r0 * ( a - m + k - 1.0D+00 ) 
     &          * ( b - m + k - 1.0D+00 ) 
     &          / ( k * ( k - m ) ) * ( 1.0D+00 - x )
              f0 = f0 + r0
            end do

            do k = 1, m
              sp0 = sp0 + 1.0D+00 / dble ( k )
            end do

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

            do k = 1, 250

              sp = sp + ( 1.0D+00 - a ) 
     &          / ( k * ( a + k - 1.0D+00 ) ) 
     &          + ( 1.0D+00 - b ) / ( k * ( b + k - 1.0D+00 ) )

              sm = 0.0D+00
              do j = 1, m
                sm = sm + 1.0D+00 / dble ( j + k )
              end do

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

              r1 = r1 * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) 
     &          / ( k * ( m + k ) ) * ( 1.0D+00 - x )

              f1 = f1 + r1 * rp

              if ( abs ( f1 - hw ) < abs ( f1 ) * eps ) then
                exit
              end if

              hw = f1

            end do

            hf = f0 * c0 + f1 * c1

          end if

        else

          ga = gamma ( a )
          gb = gamma ( b )
          gc = gamma ( c )
          gca = gamma ( c - a )
          gcb = gamma ( c - b )
          gcab = gamma ( c - a - b )
          gabc = gamma ( a + b - c )
          c0 = gc * gcab / ( gca * gcb )
          c1 = gc * gabc / ( ga * gb ) 
     &      * ( 1.0D+00 - x )**( c - a - b )
          hf = 0.0D+00
          hw = hf
          r0 = c0
          r1 = c1

          do k = 1, 250

            r0 = r0 * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) 
     &        / ( k * ( a + b - c + k ) ) * ( 1.0D+00 - x )

            r1 = r1 * ( c - a + k - 1.0D+00 ) * ( c - b + k - 1.0D+00 )
     &        / ( k * ( c - a - b + k ) ) * ( 1.0D+00 - x )

            hf = hf + r0 + r1

            if ( abs ( hf - hw ) < abs ( hf ) * eps ) then
              exit
            end if

            hw = hf

          end do

          hf = hf + c0 + c1

        end if

      else

        a0 = 1.0D+00

        if ( a < c .and. c < 2.0D+00 * a .and. b < c 
     &    .and. c < 2.0D+00 * b ) then

          a0 = ( 1.0D+00 - x )**( c - a - b )
          a = c - a
          b = c - b

        end if

        hf = 1.0D+00
        hw = hf
        r = 1.0D+00

        do k = 1, 250

          r = r * ( a + k - 1.0D+00 ) * ( b + k - 1.0D+00 ) 
     &      / ( k * ( c + k - 1.0D+00 ) ) * x

          hf = hf + r

          if ( abs ( hf - hw ) <= abs ( hf ) * eps ) then
            exit
          end if

          hw = hf

        end do

        hf = a0 * hf

      end if

      if ( x1 < 0.0D+00 ) then
        x = x1
        c0 = 1.0D+00 / ( 1.0D+00 - x ) ** aa
        hf = c0 * hf
      end if

      a = aa
      b = bb

      if ( 120 < k ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'r8_hyper_2f1(): Warning!'
        write ( *, '(a)' ) '  A large number of iterations were needed.'
        write ( *, '(a)' ) '  The results should be checked.'
      end if

      return
      end
      function r8_psi ( xx )

c*********************************************************************72
c
cc r8_psi() evaluates the function Psi(X).
c
c  Discussion:
c
c    This routine evaluates the logarithmic derivative of the
c    GAMMA function,
c
c      PSI(X) = d/dX (GAMMA(X)) / GAMMA(X) 
c             = d/dX LN ( GAMMA(X) )
c
c    for real X, where either
c
c      -XMAX1 < X < -XMIN  and X is not a negative integer), 
c
c    or
c
c      XMIN < X.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    23 January 2008
c
c  Author:
c
c    William Cody
c
c  Reference:
c
c    William Cody, Anthony Strecok, Henry Thacher,
c    Chebyshev Approximations for the Psi Function,
c    Mathematics of Computation,
c    Volume 27, Number 121, January 1973, pages 123-127.
c
c  Input:
c
c    double precision XX, the argument of the function.
c
c  Output:
c
c    double precision R8_PSI, the value of the function.
c
      implicit none

      double precision aug
      double precision den
      double precision four
      double precision fourth
      double precision half
      integer i
      integer n
      integer nq
      double precision one
      double precision p1(9)
      double precision p2(7)
      double precision piov4
      double precision q1(8)
      double precision q2(6)
      double precision r8_psi
      double precision sgn
      double precision three
      double precision xlarge
      double precision upper
      double precision w
      double precision x
      double precision xinf
      double precision xmax1
      double precision xmin1
      double precision xsmall
      double precision x01
      double precision x01d
      double precision x02
      double precision xx
      double precision z
      double precision zero
c
c  Mathematical constants.  PIOV4 = pi / 4
c
      data zero /0.0d0/
      data fourth / 0.25d0/ 
      data half / 0.5d0 /
      data one / 1.0d0 /
      data three /3.0d0/
      data four /4.0d0/
      data piov4 /7.8539816339744830962d-01/
c
c  Machine-dependent constants
c
      data xinf /1.70d+38/
      data xmin1 /5.89d-39/
      data xmax1 /3.60d+16/
      data xsmall /2.05d-09/
      data xlarge /2.04d+15/
c
c  Zero of psi(x)
c
      data x01 /187.0d0/
      data x01d /128.0d0/
      data x02 /6.9464496836234126266d-04/
c
c  Coefficients for approximation to  psi(x)/(x-x0)  over [0.5, 3.0]
c
      data p1/4.5104681245762934160d-03,5.4932855833000385356d+00,
     &        3.7646693175929276856d+02,7.9525490849151998065d+03,
     &        7.1451595818951933210d+04,3.0655976301987365674d+05,
     &        6.3606997788964458797d+05,5.8041312783537569993d+05,
     &        1.6585695029761022321d+05/
      data q1/9.6141654774222358525d+01,2.6287715790581193330d+03,
     &        2.9862497022250277920d+04,1.6206566091533671639d+05,
     &        4.3487880712768329037d+05,5.4256384537269993733d+05,
     &        2.4242185002017985252d+05,6.4155223783576225996d-08/
c
c  Coefficients for approximation to  psi(x) - ln(x) + 1/(2x)
c  for 3.0 < x.
c
      data p2/-2.7103228277757834192d+00,-1.5166271776896121383d+01,
     &        -1.9784554148719218667d+01,-8.8100958828312219821d+00,
     &        -1.4479614616899842986d+00,-7.3689600332394549911d-02,
     &        -6.5135387732718171306d-21/
      data q2/ 4.4992760373789365846d+01, 2.0240955312679931159d+02,
     &         2.4736979003315290057d+02, 1.0742543875702278326d+02,
     &         1.7463965060678569906d+01, 8.8427520398873480342d-01/

      x = xx
      w = abs ( x )
      aug = zero
c
c  Check for valid arguments, then branch to appropriate algorithm.
c
      if ( -x .ge. xmax1 .or. w .lt. xmin1 ) then
        r8_psi = xinf
        if ( zero .lt. x ) then
          r8_psi = -xinf
        end if
        return
      end if

      if ( x .ge. half ) then
        go to 200
c
c  X < 0.5, use reflection formula: psi(1-x) = psi(x) + pi * cot(pi*x)
c  Use 1/X for PI*COTAN(PI*X)  when  XMIN1 < |X| <= XSMALL.
c
      else if ( w .le. xsmall ) then
        aug = -one / x
        go to 150
      end if
c
c  Argument reduction for cotangent.
c
c 100 continue

      if ( x .lt. zero ) then
        sgn = piov4
      else
        sgn = - piov4
      end if

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

      if ( n + n .ne. nq ) then
        w = one - w
      end if

      z = piov4 * w

      if ( mod ( n, 2 ) .ne. 0 ) then
        sgn = - sgn
      end if
c
c  Determine the final value for  -pi * cotan(pi*x).
c
      n = ( nq + 1 ) / 2
      if ( mod ( n, 2 ) .eq. 0 ) then
c
c  Check for singularity.
c
        if ( z .eq. zero ) then
          r8_psi = xinf
          if ( zero .lt. x ) then
            r8_psi = -xinf
          end if
          return
        end if

        aug = sgn * ( four / tan ( z ) )

      else
        aug = sgn * ( four * tan ( z ) )
      end if

  150 continue

      x = one - x

  200 continue
c
c  0.5 <= X <= 3.0.
c
      if ( x .le. three ) then

        den = x
        upper = p1(1) * x
        do i = 1, 7
          den = ( den + q1(i) ) * x
          upper = ( upper + p1(i+1) ) * x
        end do
        den = ( upper + p1(9) ) / ( den + q1(8) )
        x = ( x - x01 / x01d ) - x02
        r8_psi = den * x + aug
        return

      end if
c
c  3.0 < X.
c
      if ( x .lt. xlarge ) then
        w = one / ( x * x )
        den = w
        upper = p2(1) * w
        do i = 1, 5
          den = ( den + q2(i) ) * w
          upper = ( upper + p2(i+1) ) * w
        end do
        aug = ( upper + p2(7) ) / ( den + q2(6) ) - half / x + aug
      end if

      r8_psi = aug + log ( x )

      return
      end
      subroutine r8vec_print ( n, a, title )

c*********************************************************************72
c
cc r8vec_print() prints an R8VEC.
c
c  Discussion:
c
c    An R8VEC is a vector of R8's.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    22 August 2000
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    integer N, the number of components of the vector.
c
c    double precision A(N), the vector to be printed.
c
c    character ( len = * ) TITLE, a title.
c
      implicit none

      integer n

      double precision a(n)
      integer i
      character ( len = * ) title

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) trim ( title )
      write ( *, '(a)' ) ' '

      do i = 1, n
        write ( *, '(2x,i8,a,1x,g16.8)' ) i, ':', a(i)
      end do

      return
      end