# include # include # include # include # include # include "gegenbauer_polynomial.h" /******************************************************************************/ bool gegenbauer_alpha_check ( double alpha ) /******************************************************************************/ /* Purpose: gegenbauer_alpha_check() checks the value of ALPHA. Licensing: This code is distributed under the MIT license. Modified: 23 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "gegenbauer_alpha_check(): Fatal error!\n" ); fprintf ( stderr, " Illegal value of ALPHA.\n" ); fprintf ( stderr, " ALPHA = %g\n", alpha ); fprintf ( stderr, " but ALPHA must be greater than -0.5.\n" ); } } return check; } /******************************************************************************/ void gegenbauer_ek_compute ( int n, double alpha, double x[], double w[] ) /******************************************************************************/ /* Purpose: gegenbauer_ek_compute() computes a Gauss-Gegenbauer quadrature rule. Discussion: The integral: integral ( -1 <= x <= +1 ) (1-x*x)^alpha * f(x) dx The quadrature rule: sum ( 1 <= i <= n ) w(i) * f ( x(i) ) Licensing: This code is distributed under the MIT license. Modified: 01 December 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. Parameters: Input, int N, the order. Input, double ALPHA, the exponent of (1-X*X) in the weight. -1.0 < ALPHA. Output, double X[N], the abscissas. Output, double W[N], the weights. */ { double abi; double *bj; bool check; int i; double zemu; /* Check N. */ if ( n < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GEGENBAUER_EK_COMPUTE - Fatal error!\n" ); fprintf ( stderr, " 1 <= N is required.\n" ); exit ( 1 ); } /* Check ALPHA. */ check = gegenbauer_alpha_check ( alpha ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GEGENBAUER_EK_COMPUTE - Fatal error!\n" ); fprintf ( stderr, " 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 = ( double * ) malloc ( n * sizeof ( double ) ); 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 ); } free ( bj ); return; } /******************************************************************************/ double gegenbauer_integral ( int expon, double alpha ) /******************************************************************************/ /* Purpose: gegenbauer_integral(): 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; bool check; double value; double value1; check = gegenbauer_alpha_check ( alpha ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GEGENBAUER_INTEGRAL - Fatal error!\n" ); fprintf ( stderr, " Illegal value of ALPHA.\n" ); exit ( 1 ); } 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; } /******************************************************************************/ double *gegenbauer_polynomial_value ( int m, int n, double alpha, double x[] ) /******************************************************************************/ /* 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GEGENBAUER_POLYNOMIAL_VALUE - Fatal error!\n" ); fprintf ( stderr, " Illegal value of ALPHA.\n" ); exit ( 1 ); } c = ( double * ) malloc ( ( m + 1 ) * n * sizeof ( double ) ); 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; } /******************************************************************************/ void gegenbauer_polynomial_values ( int *n_data, int *n, double *a, double *x, double *fx ) /******************************************************************************/ /* 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 } /******************************************************************************/ void gegenbauer_ss_compute ( int order, double alpha, double xtab[], double weight[] ) /******************************************************************************/ /* 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. 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; bool check; 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GEGENBAUER_SS_COMPUTE - Fatal error!\n" ); fprintf ( stderr, " 1 <= ORDER is required.\n" ); exit ( 1 ); } /* Check ALPHA. */ check = gegenbauer_alpha_check ( alpha ); if ( ! check ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GEGENBAUER_SS_COMPUTE - Fatal error!\n" ); fprintf ( stderr, " Illegal value of ALPHA.\n" ); exit ( 1 ); } /* Set the recursion coefficients. */ c = ( double * ) malloc ( order * sizeof ( double ) ); 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; } free ( c ); return; } /******************************************************************************/ void gegenbauer_ss_recur ( double *p2, double *dp2, double *p1, double x, int order, double alpha, double c[] ) /******************************************************************************/ /* 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. 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; } /******************************************************************************/ void gegenbauer_ss_root ( double *x, int order, double alpha, double *dp2, double *p1, double c[] ) /******************************************************************************/ /* Purpose: gegenbauer_ss_root() improves an approximate 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. 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; } /******************************************************************************/ double *gegenbauer_to_monomial_matrix ( int n, double alpha ) /******************************************************************************/ /* 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; } /******************************************************************************/ void hyper_2f1_values ( int *n_data, double *a, double *b, double *c, double *x, double *fx ) /******************************************************************************/ /* 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 } /******************************************************************************/ void imtqlx ( int n, double d[], double e[], double z[] ) /******************************************************************************/ /* 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: 11 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "imtqlx(): Fatal error!\n" ); fprintf ( stderr, " 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; } /******************************************************************************/ double r8_hyper_2f1 ( double a, double b, double c, double x ) /******************************************************************************/ /* 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. 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; int l0; int l1; int l2; int l3; int l4; int 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 || l1 ) { hf = 0.0; fprintf ( stderr, "\n" ); fprintf ( stderr, "r8_hyper_2f1():Fatal error!\n" ); fprintf ( stderr, " The hypergeometric series is divergent.\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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "r8_hyper_2f1(): Warning!\n" ); fprintf ( stderr, " A large number of iterations were needed.\n" ); fprintf ( stderr, " The accuracy of the results should be checked.\n" ); } return hf; } /******************************************************************************/ double r8_psi ( double xx ) /******************************************************************************/ /* 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; int i; int n; int nq; 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 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; x = xx; w = fabs ( x ); aug = 0.0; /* Check for valid arguments, then branch to appropriate algorithm. */ if ( xmax1 <= - x || w < xmin1 ) { if ( 0.0 < x ) { value = - xinf; } else { value = xinf; } return value; } if ( x < 0.5 ) { /* 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 = - 1.0 / x; } /* Argument reduction for cotangent. */ else { if ( x < 0.0 ) { sgn = piov4; } else { sgn = - piov4; } w = w - ( double ) ( ( int ) ( w ) ); nq = ( int ) ( w * 4.0 ); w = 4.0 * ( w - ( double ) ( nq ) * 0.25 ); /* 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 = 1.0 - 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 == 0.0 ) { if ( 0.0 < x ) { value = -xinf; } else { value = xinf; } return value; } aug = sgn * ( 4.0 / tan ( z ) ); } else { aug = sgn * ( 4.0 * tan ( z ) ); } } x = 1.0 - x; } /* 0.5 <= X <= 3.0. */ if ( x <= 3.0 ) { 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 = 1.0 / ( 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] ) - 0.5 / x + aug; } value = aug + log ( x ); return value; } /******************************************************************************/ double r8_sign ( double x ) /******************************************************************************/ /* Purpose: r8_sign() returns the sign of an R8. Licensing: This code is distributed under the MIT license. Modified: 08 May 2006 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; } /******************************************************************************/ double r8_uniform_ab ( double a, double b, int *seed ) /******************************************************************************/ /* Purpose: r8_uniform_ab() returns a pseudorandom R8 scaled to [A,B]. Licensing: This code is distributed under the MIT license. Modified: 21 November 2004 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 r; double value; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; value = a + ( b - a ) * r; return value; } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8vec_print() prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Input: int N, the number of components of the vector. double A[N], the vector to be printed. char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14g\n", i, a[i] ); } return; }