# 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 GNU LGPL license. Modified: 23 November 2015 Author: John Burkardt Parameters: 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 GNU LGPL 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 ) * r8_gamma ( alpha + 1.0 ) * r8_gamma ( alpha + 1.0 ) / r8_gamma ( 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 GNU LGPL 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 = r8_gamma ( 1.0 + c ) * 2.0 * r8_gamma ( 1.0 + alpha ) * value1 / r8_gamma ( 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 GNU LGPL 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 GNU LGPL license. Modified: 06 August 2004 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Parameters: Input/output, int *N_DATA. The user sets N_DATA to 0 before the first call. On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. Output, int *N, the order parameter of the function. Output, double *A, the real parameter of the function. Output, double *X, the argument of the function. Output, double *FX, the value of the function. */ { # define N_MAX 38 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 GNU LGPL 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 = r8_gamma ( alpha + 1.0 ) * r8_gamma ( alpha + 1.0 ) / r8_gamma ( 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 GNU LGPL 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 GNU LGPL 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 = r8_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; } /******************************************************************************/ 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 GNU LGPL 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 GNU LGPL license. Modified: 11 January 2010 Author: Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky. C version by John Burkardt. Reference: Sylvan Elhay, Jaroslav Kautsky, Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of Interpolatory Quadrature, ACM Transactions on Mathematical Software, Volume 13, Number 4, December 1987, pages 399-415. 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 = r8_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_epsilon ( ) /******************************************************************************/ /* Purpose: R8_EPSILON returns the R8 round off unit. Discussion: R8_EPSILON is a number R which is a power of 2 with the property that, to the precision of the computer's arithmetic, 1 < 1 + R but 1 = ( 1 + R / 2 ) Licensing: This code is distributed under the GNU LGPL license. Modified: 01 September 2012 Author: John Burkardt Parameters: Output, double R8_EPSILON, the R8 round-off unit. */ { const double value = 2.220446049250313E-016; return value; } /******************************************************************************/ double r8_gamma ( double x ) /******************************************************************************/ /* Purpose: R8_GAMMA evaluates Gamma(X) for a real argument. Discussion: This routine calculates the gamma function for a real argument X. Computation is based on an algorithm outlined in reference 1. The program uses rational functions that approximate the gamma function to at least 20 significant decimal digits. Coefficients for the approximation over the interval (1,2) are unpublished. Those for the approximation for 12 <= X are from reference 2. Licensing: This code is distributed under the GNU LGPL license. Modified: 18 January 2008 Author: Original FORTRAN77 version by William Cody, Laura Stoltz. C version by John Burkardt. Reference: William Cody, An Overview of Software Development for Special Functions, in Numerical Analysis Dundee, 1975, edited by GA Watson, Lecture Notes in Mathematics 506, Springer, 1976. John Hart, Ward Cheney, Charles Lawson, Hans Maehly, Charles Mesztenyi, John Rice, Henry Thatcher, Christoph Witzgall, Computer Approximations, Wiley, 1968, LC: QA297.C64. Parameters: Input, double X, the argument of the function. Output, double R8_GAMMA, the value of the function. */ { double c[7] = { -1.910444077728E-03, 8.4171387781295E-04, -5.952379913043012E-04, 7.93650793500350248E-04, -2.777777777777681622553E-03, 8.333333333333333331554247E-02, 5.7083835261E-03 }; double eps = 2.22E-16; double fact; double half = 0.5; int i; int n; double one = 1.0; double p[8] = { -1.71618513886549492533811E+00, 2.47656508055759199108314E+01, -3.79804256470945635097577E+02, 6.29331155312818442661052E+02, 8.66966202790413211295064E+02, -3.14512729688483675254357E+04, -3.61444134186911729807069E+04, 6.64561438202405440627855E+04 }; int parity; const double r8_pi = 3.1415926535897932384626434; double q[8] = { -3.08402300119738975254353E+01, 3.15350626979604161529144E+02, -1.01515636749021914166146E+03, -3.10777167157231109440444E+03, 2.25381184209801510330112E+04, 4.75584627752788110767815E+03, -1.34659959864969306392456E+05, -1.15132259675553483497211E+05 }; double res; double sqrtpi = 0.9189385332046727417803297; double sum; double twelve = 12.0; double two = 2.0; double value; double xbig = 171.624; double xden; double xinf = 1.79E+308; double xminin = 2.23E-308; double xnum; double y; double y1; double ysq; double z; double zero = 0.0;; parity = 0; fact = 1.0; n = 0; y = x; /* Argument is negative. */ if ( y <= zero ) { y = - x; y1 = ( double ) ( int ) ( y ); res = y - y1; if ( res != zero ) { if ( y1 != ( double ) ( int ) ( y1 * half ) * two ) { parity = 1; } fact = - r8_pi / sin ( r8_pi * res ); y = y + one; } else { res = xinf; value = res; return value; } } /* Argument is positive. */ if ( y < eps ) { /* Argument < EPS. */ if ( xminin <= y ) { res = one / y; } else { res = xinf; value = res; return value; } } else if ( y < twelve ) { y1 = y; /* 0.0 < argument < 1.0. */ if ( y < one ) { z = y; y = y + one; } /* 1.0 < argument < 12.0. Reduce argument if necessary. */ else { n = ( int ) ( y ) - 1; y = y - ( double ) ( n ); z = y - one; } /* Evaluate approximation for 1.0 < argument < 2.0. */ xnum = zero; xden = one; for ( i = 0; i < 8; i++ ) { xnum = ( xnum + p[i] ) * z; xden = xden * z + q[i]; } res = xnum / xden + one; /* Adjust result for case 0.0 < argument < 1.0. */ if ( y1 < y ) { res = res / y1; } /* Adjust result for case 2.0 < argument < 12.0. */ else if ( y < y1 ) { for ( i = 1; i <= n; i++ ) { res = res * y; y = y + one; } } } else { /* Evaluate for 12.0 <= argument. */ if ( y <= xbig ) { ysq = y * y; sum = c[6]; for ( i = 0; i < 6; i++ ) { sum = sum / ysq + c[i]; } sum = sum / y - y + sqrtpi; sum = sum + ( y - half ) * log ( y ); res = exp ( sum ); } else { res = xinf; value = res; return value; } } /* Final adjustments and return. */ if ( parity ) { res = - res; } if ( fact != one ) { res = fact / res; } value = res; return value; } /******************************************************************************/ 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 GNU LGPL license. Modified: 05 July 2009 Author: Original FORTRAN77 version by Shanjie Zhang, Jianming Jin. C version by John Burkardt. The original FORTRAN77 version of this routine is copyrighted by Shanjie Zhang and Jianming Jin. However, they give permission to incorporate this routine into a user program provided that the copyright is acknowledged. Reference: Shanjie Zhang, Jianming Jin, Computation of Special Functions, Wiley, 1996, ISBN: 0-471-11963-6, LC: QA351.C45 Parameters: Input, double A, B, C, X, the arguments of the function. C must not be equal to a nonpositive integer. X < 1. Output, double R8_HYPER_2F1, the value of the function. */ { double a0; double aa; double bb; double c0; double c1; double el = 0.5772156649015329; double eps; double f0; double f1; double g0; double g1; double g2; double g3; double ga; double gabc; double gam; double gb; double gbm; double gc; double gca; double gcab; double gcb; double gm; double hf; double hw; int j; int k; 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 = r8_gamma ( c ); gcab = r8_gamma ( c - a - b ); gca = r8_gamma ( c - a ); gcb = r8_gamma ( 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 = r8_gamma ( c ); g2 = r8_gamma ( 1.0 + a / 2.0 - b ); g3 = r8_gamma ( 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 = r8_gamma ( a ); gb = r8_gamma ( b ); gc = r8_gamma ( c ); gam = r8_gamma ( a + m ); gbm = r8_gamma ( 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 = r8_gamma ( a ); gb = r8_gamma ( b ); gc = r8_gamma ( c ); gca = r8_gamma ( c - a ); gcb = r8_gamma ( c - b ); gcab = r8_gamma ( c - a - b ); gabc = r8_gamma ( 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 GNU LGPL license. Modified: 09 February 2008 Author: Original FORTRAN77 version by William Cody. C 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 GNU LGPL 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 GNU LGPL 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 GNU LGPL license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, 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; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the GNU LGPL license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE }