# include # include # include # include # include # include using namespace std; # include "test_values.hpp" //****************************************************************************80 void abram0_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ABRAM0_VALUES returns some values of the Abramowitz0 function. // // Discussion: // // The function is defined by: // // ABRAM0(X) = integral ( 0 <= T < +oo ) exp ( -T * T - X / T ) dT // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 21 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. // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.87377726306985360531E+00, 0.84721859650456925922E+00, 0.77288934483988301615E+00, 0.59684345853450151603E+00, 0.29871735283675888392E+00, 0.15004596450516388138E+00, 0.11114662419157955096E+00, 0.83909567153151897766E-01, 0.56552321717943417515E-01, 0.49876496603033790206E-01, 0.44100889219762791328E-01, 0.19738535180254062496E-01, 0.86193088287161479900E-02, 0.40224788162540127227E-02, 0.19718658458164884826E-02, 0.10045868340133538505E-02, 0.15726917263304498649E-03, 0.10352666912350263437E-04, 0.91229759190956745069E-06, 0.25628287737952698742E-09 }; static double x_vec[N_MAX] = { 0.0019531250E+00, 0.0078125000E+00, 0.0312500000E+00, 0.1250000000E+00, 0.5000000000E+00, 1.0000000000E+00, 1.2500000000E+00, 1.5000000000E+00, 1.8750000000E+00, 2.0000000000E+00, 2.1250000000E+00, 3.0000000000E+00, 4.0000000000E+00, 5.0000000000E+00, 6.0000000000E+00, 7.0000000000E+00, 10.0000000000E+00, 15.0000000000E+00, 20.0000000000E+00, 40.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void abram1_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ABRAM1_VALUES returns some values of the Abramowitz1 function. // // Discussion: // // The function is defined by: // // ABRAM1(x) = integral ( 0 <= t < oo ) t * exp ( -t^2 - x / t ) dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 21 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. // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.49828219848799921792E+00, 0.49324391773047288556E+00, 0.47431612784691234649E+00, 0.41095983258760410149E+00, 0.25317617388227035867E+00, 0.14656338138597777543E+00, 0.11421547056018366587E+00, 0.90026307383483764795E-01, 0.64088214170742303375E-01, 0.57446614314166191085E-01, 0.51581624564800730959E-01, 0.25263719555776416016E-01, 0.11930803330196594536E-01, 0.59270542280915272465E-02, 0.30609215358017829567E-02, 0.16307382136979552833E-02, 0.28371851916959455295E-03, 0.21122150121323238154E-04, 0.20344578892601627337E-05, 0.71116517236209642290E-09 }; static double x_vec[N_MAX] = { 0.0019531250E+00, 0.0078125000E+00, 0.0312500000E+00, 0.1250000000E+00, 0.5000000000E+00, 1.0000000000E+00, 1.2500000000E+00, 1.5000000000E+00, 1.8750000000E+00, 2.0000000000E+00, 2.1250000000E+00, 3.0000000000E+00, 4.0000000000E+00, 5.0000000000E+00, 6.0000000000E+00, 7.0000000000E+00, 10.0000000000E+00, 15.0000000000E+00, 20.0000000000E+00, 40.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void abram2_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ABRAM2_VALUES returns some values of the Abramowitz2 function. // // Discussion: // // The function is defined by: // // ABRAM2(x) = Integral ( 0 <= t < +oo ) t^2 * exp( -t^2 - x / t ) dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 22 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. // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.44213858162107913430E+00, 0.43923379545684026308E+00, 0.42789857297092602234E+00, 0.38652825661854504406E+00, 0.26538204413231368110E+00, 0.16848734838334595000E+00, 0.13609200032513227112E+00, 0.11070330027727917352E+00, 0.82126019995530382267E-01, 0.74538781999594581763E-01, 0.67732034377612811390E-01, 0.35641808698811851022E-01, 0.17956589956618269083E-01, 0.94058737143575370625E-02, 0.50809356204299213556E-02, 0.28149565414209719359E-02, 0.53808696422559303431E-03, 0.44821756380146327259E-04, 0.46890678427324100410E-05, 0.20161544850996420504E-08 }; static double x_vec[N_MAX] = { 0.0019531250E+00, 0.0078125000E+00, 0.0312500000E+00, 0.1250000000E+00, 0.5000000000E+00, 1.0000000000E+00, 1.2500000000E+00, 1.5000000000E+00, 1.8750000000E+00, 2.0000000000E+00, 2.1250000000E+00, 3.0000000000E+00, 4.0000000000E+00, 5.0000000000E+00, 6.0000000000E+00, 7.0000000000E+00, 10.0000000000E+00, 15.0000000000E+00, 20.0000000000E+00, 40.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void agm_values ( int &n_data, double &a, double &b, double &fx ) //****************************************************************************80 // // Purpose: // // AGM_VALUES returns some values of the AGM. // // Discussion: // // The AGM is defined for nonnegative A and B. // // The AGM of numbers A and B is defined by setting // // A(0) = A, // B(0) = B // // A(N+1) = ( A(N) + B(N) ) / 2 // B(N+1) = sqrt ( A(N) * B(N) ) // // The two sequences both converge to AGM(A,B). // // In Mathematica, the AGM can be evaluated by // // ArithmeticGeometricMean [ a, b ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 09 February 2008 // // Author: // // John Burkardt // // Reference: // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &A, &B, the argument ofs the function. // // Output, double &FX, the value of the function. // { # define N_MAX 15 static double a_vec[N_MAX] = { 22.0, 83.0, 42.0, 26.0, 4.0, 6.0, 40.0, 80.0, 90.0, 9.0, 53.0, 1.0, 1.0, 1.0, 1.5 }; static double b_vec[N_MAX] = { 96.0, 56.0, 7.0, 11.0, 63.0, 45.0, 75.0, 0.0, 35.0, 1.0, 53.0, 2.0, 4.0, 8.0, 8.0 }; static double fx_vec[N_MAX] = { 52.274641198704240049, 68.836530059858524345, 20.659301196734009322, 17.696854873743648823, 23.867049721753300163, 20.717015982805991662, 56.127842255616681863, 0.000000000000000000, 59.269565081229636528, 3.9362355036495554780, 53.000000000000000000, 1.4567910310469068692, 2.2430285802876025701, 3.6157561775973627487, 4.0816924080221632670 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; a = 0.0; b = 0.0; fx = 0.0; } else { a = a_vec[n_data-1]; b = b_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_ai_values ( int &n_data, double &x, double &ai ) //****************************************************************************80 // // Purpose: // // AIRY_AI_VALUES returns some values of the Airy Ai(x) function. // // Discussion: // // The Airy functions Ai(X) and Bi(X) are a pair of linearly independent // solutions of the differential equation: // // W'' - X * W = 0; // // In Mathematica, the function can be evaluated by: // // AiryAi[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 11 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, double &X, the argument of the function. // // Output, double &AI, the value of the Airy AI function. // { # define N_MAX 11 static double ai_vec[N_MAX] = { 0.3550280538878172E+00, 0.3292031299435381E+00, 0.3037031542863820E+00, 0.2788064819550049E+00, 0.2547423542956763E+00, 0.2316936064808335E+00, 0.2098000616663795E+00, 0.1891624003981501E+00, 0.1698463174443649E+00, 0.1518868036405444E+00, 0.1352924163128814E+00 }; static double x_vec[N_MAX] = { 0.0E+00, 0.1E+00, 0.2E+00, 0.3E+00, 0.4E+00, 0.5E+00, 0.6E+00, 0.7E+00, 0.8E+00, 0.9E+00, 1.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; ai = 0.0; } else { x = x_vec[n_data-1]; ai = ai_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_ai_int_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // AIRY_AI_INT_VALUES returns some values of the integral of the Airy function. // // Discussion: // // The function is defined by: // // AIRY_AI_INT(x) = Integral ( 0 <= t <= x ) Ai(t) dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 22 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. // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { -0.75228838916610124300E+00, -0.57348350185854889466E+00, -0.76569840313421291743E+00, -0.65181015505382467421E+00, -0.55881974894471876922E+00, -0.56902352870716815309E+00, -0.47800749642926168100E+00, -0.46567398346706861416E+00, -0.96783140945618013679E-01, -0.34683049857035607494E-03, 0.34658366917927930790E-03, 0.27657581846051227124E-02, 0.14595330491185717833E+00, 0.23631734191710977960E+00, 0.33289264538612212697E+00, 0.33318759129779422976E+00, 0.33332945170523851439E+00, 0.33333331724248357420E+00, 0.33333333329916901594E+00, 0.33333333333329380187E+00 }; static double x_vec[N_MAX] = { -12.0000000000E+00, -11.0000000000E+00, -10.0000000000E+00, -9.5000000000E+00, -9.0000000000E+00, -6.5000000000E+00, -4.0000000000E+00, -1.0000000000E+00, -0.2500000000E+00, -0.0009765625E+00, 0.0009765625E+00, 0.0078125000E+00, 0.5000000000E+00, 1.0000000000E+00, 4.0000000000E+00, 4.5000000000E+00, 6.0000000000E+00, 8.0000000000E+00, 10.0000000000E+00, 12.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_ai_prime_values ( int &n_data, double &x, double &aip ) //****************************************************************************80 // // Purpose: // // AIRY_AI_PRIME_VALUES returns some values of the Airy function Ai'(x). // // Discussion: // // The Airy functions Ai(X) and Bi(X) are a pair of linearly independent // solutions of the differential equation: // // W'' - X * W = 0; // // In Mathematica, the function can be evaluated by: // // AiryAiPrime[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 11 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, double &X, the argument of the function. // // Output, double &AIP, the derivative of the Airy AI function. // { # define N_MAX 11 static double aip_vec[N_MAX] = { -0.2588194037928068E+00, -0.2571304219075862E+00, -0.2524054702856195E+00, -0.2451463642190548E+00, -0.2358320344192082E+00, -0.2249105326646839E+00, -0.2127932593891585E+00, -0.1998511915822805E+00, -0.1864128638072717E+00, -0.1727638434616347E+00, -0.1591474412967932E+00 }; static double x_vec[N_MAX] = { 0.0E+00, 0.1E+00, 0.2E+00, 0.3E+00, 0.4E+00, 0.5E+00, 0.6E+00, 0.7E+00, 0.8E+00, 0.9E+00, 1.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; aip = 0.0; } else { x = x_vec[n_data-1]; aip = aip_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_bi_values ( int &n_data, double &x, double &bi ) //****************************************************************************80 // // Purpose: // // AIRY_BI_VALUES returns some values of the Airy Bi(x) function. // // Discussion: // // The Airy functions Ai(X) and Bi(X) are a pair of linearly independent // solutions of the differential equation: // // W'' - X * W = 0; // // In Mathematica, the function can be evaluated by: // // AiryBi[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 11 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, double &X, the argument of the function. // // Output, double &BI, the value of the Airy BI function. // { # define N_MAX 11 static double bi_vec[N_MAX] = { 0.6149266274460007E+00, 0.6598616901941892E+00, 0.7054642029186612E+00, 0.7524855850873156E+00, 0.8017730000135972E+00, 0.8542770431031555E+00, 0.9110633416949405E+00, 0.9733286558781659E+00, 0.1042422171231561E+01, 0.1119872813134447E+01, 0.1207423594952871E+01 }; static double x_vec[N_MAX] = { 0.0E+00, 0.1E+00, 0.2E+00, 0.3E+00, 0.4E+00, 0.5E+00, 0.6E+00, 0.7E+00, 0.8E+00, 0.9E+00, 1.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; bi = 0.0; } else { x = x_vec[n_data-1]; bi = bi_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_bi_int_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // AIRY_BI_INT_VALUES returns some values of the integral of the Airy function. // // Discussion: // // The function is defined by: // // AIRY_BI_INT(x) = Integral ( 0 <= t <= x ) Bi(t) dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 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. // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.17660819031554631869E-01, -0.15040424806140020451E-01, 0.14756446293227661920E-01, -0.11847304264848446271E+00, -0.64916741266165856037E-01, 0.97260832464381044540E-01, 0.50760058495287539119E-01, -0.37300500963429492179E+00, -0.13962988442666578531E+00, -0.12001735266723296160E-02, 0.12018836117890354598E-02, 0.36533846550952011043E+00, 0.87276911673800812196E+00, 0.48219475263803429675E+02, 0.44006525804904178439E+06, 0.17608153976228301458E+07, 0.73779211705220007228E+07, 0.14780980310740671617E+09, 0.97037614223613433849E+11, 0.11632737638809878460E+15 }; static double x_vec[N_MAX] = { -12.0000000000E+00, -10.0000000000E+00, -8.0000000000E+00, -7.5000000000E+00, -7.0000000000E+00, -6.5000000000E+00, -4.0000000000E+00, -1.0000000000E+00, -0.2500000000E+00, -0.0019531250E+00, 0.0019531250E+00, 0.5000000000E+00, 1.0000000000E+00, 4.0000000000E+00, 8.0000000000E+00, 8.5000000000E+00, 9.0000000000E+00, 10.0000000000E+00, 12.0000000000E+00, 14.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_bi_prime_values ( int &n_data, double &x, double &bip ) //****************************************************************************80 // // Purpose: // // AIRY_BI_PRIME_VALUES returns some values of the Airy function Bi'(x). // // Discussion: // // The Airy functions Ai(X) and Bi(X) are a pair of linearly independent // solutions of the differential equation: // // W'' - X * W = 0; // // In Mathematica, the function can be evaluated by: // // AiryBiPrime[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 11 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, double &X, the argument of the function. // // Output, double &BIP, the derivative of the Airy BI function. // { # define N_MAX 11 static double bip_vec[N_MAX] = { 0.4482883573538264E+00, 0.4515126311496465E+00, 0.4617892843621509E+00, 0.4800490287524480E+00, 0.5072816760506224E+00, 0.5445725641405923E+00, 0.5931444786342857E+00, 0.6544059191721400E+00, 0.7300069016152518E+00, 0.8219038903072090E+00, 0.9324359333927756E+00 }; static double x_vec[N_MAX] = { 0.0E+00, 0.1E+00, 0.2E+00, 0.3E+00, 0.4E+00, 0.5E+00, 0.6E+00, 0.7E+00, 0.8E+00, 0.9E+00, 1.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; bip = 0.0; } else { x = x_vec[n_data-1]; bip = bip_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_cai_values ( int &n_data, complex &x, complex &cai ) //****************************************************************************80 // // Purpose: // // AIRY_CAI_VALUES returns some values of the Airy Ai(x) for complex argument. // // Discussion: // // The Airy functions Ai(X) and Bi(X) are a pair of linearly independent // solutions of the differential equation: // // W'' - X * W = 0; // // In Mathematica, the function can be evaluated by: // // AiryAi[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 13 April 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, complex &X, the argument of the function. // // Output, complex &CAI, the value of the Airy AI function. // { # define N_MAX 10 static complex cai_vec[N_MAX] = { complex ( 0.1352924163128814, + 0.0000000000000000 ), complex ( 0.1433824486882056, - 0.1092193342707378 ), complex ( 0.2215404472324631, - 0.2588711788891803 ), complex ( 0.4763929771766866, - 0.3036484220291284 ), complex ( 0.5983692170633874, - 0.08154602160771214 ), complex ( 0.5355608832923521, + 0.00000000000000000 ), complex ( 0.5983692170633874, + 0.08154602160771214 ), complex ( 0.4763929771766866, + 0.3036484220291284 ), complex ( 0.2215404472324631, + 0.2588711788891803 ), complex ( 0.1433824486882056, + 0.1092193342707378 ) }; static complex x_vec[N_MAX] = { complex ( 1.0000000000000000, + 0.0000000000000000 ), complex ( 0.8090169943749474, + 0.5877852522924731 ), complex ( 0.3090169943749474, + 0.9510565162951536 ), complex ( -0.3090169943749474, + 0.9510565162951536 ), complex ( -0.8090169943749474, + 0.5877852522924731 ), complex ( -1.0000000000000000, + 0.0000000000000000 ), complex ( -0.8090169943749474, - 0.5877852522924731 ), complex ( -0.3090169943749474, - 0.9510565162951536 ), complex ( 0.3090169943749474, - 0.9510565162951536 ), complex ( 0.8090169943749474, - 0.5877852522924731 ) }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; cai = 0.0; } else { x = x_vec[n_data-1]; cai = cai_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_cbi_values ( int &n_data, complex &x, complex &cbi ) //****************************************************************************80 // // Purpose: // // AIRY_CBI_VALUES returns some values of the Airy Bi(x) for complex argument. // // Discussion: // // The Airy functions Ai(X) and Bi(X) are a pair of linearly independent // solutions of the differential equation: // // W'' - X * W = 0; // // In Mathematica, the function can be evaluated by: // // AiryAi[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 13 April 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, complex &X, the argument of the function. // // Output, complex &CBI, the value of the Airy BI function. // { # define N_MAX 10 static complex cbi_vec[N_MAX] = { complex ( 1.207423594952871, + 0.0000000000000000 ), complex ( 0.9127160108293936, + 0.3800456133135556 ), complex ( 0.6824453575635721, + 0.3343047153635002 ), complex ( 0.5726265660086474, + 0.3988641086982559 ), complex ( 0.2511841251049547, + 0.3401447690712719 ), complex ( 0.1039973894969446, + 0.0000000000000000 ), complex ( 0.2511841251049547, - 0.3401447690712719 ), complex ( 0.5726265660086474, - 0.3988641086982559 ), complex ( 0.6824453575635721, - 0.3343047153635002 ), complex ( 0.9127160108293936, - 0.3800456133135556 ) }; static complex x_vec[N_MAX] = { complex ( 1.0000000000000000, + 0.0000000000000000 ), complex ( 0.8090169943749474, + 0.5877852522924731 ), complex ( 0.3090169943749474, + 0.9510565162951536 ), complex ( -0.3090169943749474, + 0.9510565162951536 ), complex ( -0.8090169943749474, + 0.5877852522924731 ), complex ( -1.0000000000000000, + 0.0000000000000000 ), complex ( -0.8090169943749474, - 0.5877852522924731 ), complex ( -0.3090169943749474, - 0.9510565162951536 ), complex ( 0.3090169943749474, - 0.9510565162951536 ), complex ( 0.8090169943749474, - 0.5877852522924731 ) }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; cbi = 0.0; } else { x = x_vec[n_data-1]; cbi = cbi_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_gi_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // AIRY_GI_VALUES returns some values of the Airy Gi function. // // Discussion: // // The function is defined by: // // AIRY_GI(x) = Integral ( 0 <= t < +oo ) sin ( x*t+t^3/3) dt / pi // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 24 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. // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.20468308070040542435E+00, 0.18374662832557904078E+00, -0.11667221729601528265E+00, 0.31466934902729557596E+00, -0.37089040722426257729E+00, -0.25293059772424019694E+00, 0.28967410658692701936E+00, -0.34644836492634090590E+00, 0.28076035913873049496E+00, 0.21814994508094865815E+00, 0.20526679000810503329E+00, 0.22123695363784773258E+00, 0.23521843981043793760E+00, 0.82834303363768729338E-01, 0.45757385490989281893E-01, 0.44150012014605159922E-01, 0.39951133719508907541E-01, 0.35467706833949671483E-01, 0.31896005100679587981E-01, 0.26556892713512410405E-01 }; static double x_vec[N_MAX] = { -0.0019531250E+00, -0.1250000000E+00, -1.0000000000E+00, -4.0000000000E+00, -8.0000000000E+00, -8.2500000000E+00, -9.0000000000E+00, -10.0000000000E+00, -11.0000000000E+00, -13.0000000000E+00, 0.0019531250E+00, 0.1250000000E+00, 1.0000000000E+00, 4.0000000000E+00, 7.0000000000E+00, 7.2500000000E+00, 8.0000000000E+00, 9.0000000000E+00, 10.0000000000E+00, 12.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void airy_hi_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // AIRY_HI_VALUES returns some values of the Airy Hi function. // // Discussion: // // The function is defined by: // // AIRY_HI(x) = Integral ( 0 <= t < +oo ) exp(x*t-t^3/3) dt / pi // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 24 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. // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.40936798278458884024E+00, 0.37495291608048868619E+00, 0.22066960679295989454E+00, 0.77565356679703713590E-01, 0.39638826473124717315E-01, 0.38450072575004151871E-01, 0.35273216868317898556E-01, 0.31768535282502272742E-01, 0.28894408288051391369E-01, 0.24463284011678541180E-01, 0.41053540139998941517E+00, 0.44993502381204990817E+00, 0.97220515514243332184E+00, 0.83764237105104371193E+02, 0.80327744952044756016E+05, 0.15514138847749108298E+06, 0.11995859641733262114E+07, 0.21472868855967642259E+08, 0.45564115351632913590E+09, 0.32980722582904761929E+12 }; static double x_vec[N_MAX] = { -0.0019531250E+00, -0.1250000000E+00, -1.0000000000E+00, -4.0000000000E+00, -8.0000000000E+00, -8.2500000000E+00, -9.0000000000E+00, -10.0000000000E+00, -11.0000000000E+00, -13.0000000000E+00, 0.0019531250E+00, 0.1250000000E+00, 1.0000000000E+00, 4.0000000000E+00, 7.0000000000E+00, 7.2500000000E+00, 8.0000000000E+00, 9.0000000000E+00, 10.0000000000E+00, 12.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void arccos_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ARCCOS_VALUES returns some values of the arc cosine function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // ArcCos[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 12 June 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 12 static double fx_vec[N_MAX] = { 1.6709637479564564156, 1.5707963267948966192, 1.4706289056333368229, 1.3694384060045658278, 1.2661036727794991113, 1.1592794807274085998, 1.0471975511965977462, 0.92729521800161223243, 0.79539883018414355549, 0.64350110879328438680, 0.45102681179626243254, 0.00000000000000000000 }; static double x_vec[N_MAX] = { -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void arccosh_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ARCCOSH_VALUES returns some values of the hyperbolic arc cosine function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // ArcCosh[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 June 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 15 static double fx_vec[N_MAX] = { 0.0000000000000000000, 0.14130376948564857735, 0.44356825438511518913, 0.62236250371477866781, 0.75643291085695958624, 0.86701472649056510395, 0.96242365011920689500, 1.3169578969248167086, 1.7627471740390860505, 1.8115262724608531070, 2.0634370688955605467, 2.2924316695611776878, 2.9932228461263808979, 5.2982923656104845907, 7.6009022095419886114 }; static double x_vec[N_MAX] = { 1.0, 1.01, 1.1, 1.2, 1.3, 1.4, 1.5, 2.0, 3.0, 3.1415926535897932385, 4.0, 5.0, 10.0, 100.0, 1000.0 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void arcsin_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ARCSIN_VALUES returns some values of the arc sine function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // ArcSin[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 12 June 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 12 static double fx_vec[N_MAX] = { -0.10016742116155979635, 0.00000000000000000000, 0.10016742116155979635, 0.20135792079033079146, 0.30469265401539750797, 0.41151684606748801938, 0.52359877559829887308, 0.64350110879328438680, 0.77539749661075306374, 0.92729521800161223243, 1.1197695149986341867, 1.5707963267948966192 }; static double x_vec[N_MAX] = { -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void arcsinh_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ARCSINH_VALUES returns some values of the hyperbolic arc sine function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // ArcSinh[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 June 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { -2.3124383412727526203, -0.88137358701954302523, 0.00000000000000000000, 0.099834078899207563327, 0.19869011034924140647, 0.29567304756342243910, 0.39003531977071527608, 0.48121182505960344750, 0.56882489873224753010, 0.65266656608235578681, 0.73266825604541086415, 0.80886693565278246251, 0.88137358701954302523, 1.4436354751788103425, 1.8184464592320668235, 2.0947125472611012942, 2.3124383412727526203, 2.9982229502979697388, 5.2983423656105887574, 7.6009027095419886115 }; static double x_vec[N_MAX] = { -5.0, -1.0, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 100.0, 1000.0 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void arctan_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ARCTAN_VALUES returns some values of the arc tangent function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // ArcTan[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 12 June 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 11 static double f_vec[N_MAX] = { 0.00000000000000000000, 0.24497866312686415417, 0.32175055439664219340, 0.46364760900080611621, 0.78539816339744830962, 1.1071487177940905030, 1.2490457723982544258, 1.3258176636680324651, 1.3734007669450158609, 1.4711276743037345919, 1.5208379310729538578 }; static double x_vec[N_MAX] = { 0.00000000000000000000, 0.25000000000000000000, 0.33333333333333333333, 0.50000000000000000000, 1.0000000000000000000, 2.0000000000000000000, 3.0000000000000000000, 4.0000000000000000000, 5.0000000000000000000, 10.000000000000000000, 20.000000000000000000 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = f_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void arctan_int_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ARCTAN_INT_VALUES returns some values of the inverse tangent integral. // // Discussion: // // The function is defined by: // // ARCTAN_INT(x) = Integral ( 0 <= t <= x ) arctan ( t ) / t dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 25 August 2004 // // Author: // // John Burkardt // // Reference: // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double f_vec[N_MAX] = { 0.19531241721588483191E-02, -0.39062433772980711281E-02, 0.78124470192576499535E-02, 0.15624576181996527280E-01, -0.31246610349485401551E-01, 0.62472911335014397321E-01, 0.12478419717389654039E+00, -0.24830175098230686908E+00, 0.48722235829452235711E+00, 0.91596559417721901505E+00, 0.12749694484943800618E+01, -0.15760154034463234224E+01, 0.24258878412859089996E+01, 0.33911633326292997361E+01, 0.44176450919422186583E+01, -0.47556713749547247774E+01, 0.50961912150934111303E+01, 0.53759175735714876256E+01, -0.61649904785027487422E+01, 0.72437843013083534973E+01 }; static double x_vec[N_MAX] = { 0.0019531250E+00, -0.0039062500E+00, 0.0078125000E+00, 0.0156250000E+00, -0.0312500000E+00, 0.0625000000E+00, 0.1250000000E+00, -0.2500000000E+00, 0.5000000000E+00, 1.0000000000E+00, 1.5000000000E+00, -2.0000000000E+00, 4.0000000000E+00, 8.0000000000E+00, 16.0000000000E+00, -20.0000000000E+00, 25.0000000000E+00, 30.0000000000E+00, -50.0000000000E+00, 100.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = f_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void arctan2_values ( int &n_data, double &x, double &y, double &fxy ) //****************************************************************************80 // // Purpose: // // ARCTAN2_VALUES: arc tangent function of two arguments. // // Discussion: // // In Mathematica, the function can be evaluated by: // // ArcTan[x,y] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 29 March 2010 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, &Y, the arguments of the function. // // Output, double &FXY, the value of the function. // { # define N_MAX 19 static double f_vec[N_MAX] = { -1.5707963267948966192, -1.0471975511965977462, -0.52359877559829887308, 0.00000000000000000000, 0.52359877559829887308, 1.0471975511965977462, 1.5707963267948966192, 2.0943951023931954923, 2.6179938779914943654, 3.1415926535897932385, -2.6179938779914943654, -2.0943951023931954923, -1.5707963267948966192, -1.0471975511965977462, -0.52359877559829887308, 0.00000000000000000000, 0.52359877559829887308, 1.0471975511965977462, 1.5707963267948966192 }; static double x_vec[N_MAX] = { 0.00000000000000000000, 0.50000000000000000000, 0.86602540378443864676, 1.00000000000000000000, 0.86602540378443864676, 0.50000000000000000000, 0.00000000000000000000, -0.50000000000000000000, -0.86602540378443864676, -1.00000000000000000000, -0.86602540378443864676, -0.50000000000000000000, 0.00000000000000000000, 0.50000000000000000000, 0.86602540378443864676, 1.00000000000000000000, 0.86602540378443864676, 0.50000000000000000000, 0.00000000000000000000 }; static double y_vec[N_MAX] = { -1.00000000000000000000, -0.86602540378443864676, -0.50000000000000000000, 0.00000000000000000000, 0.50000000000000000000, 0.86602540378443864676, 1.00000000000000000000, 0.86602540378443864676, 0.50000000000000000000, 0.00000000000000000000, -0.50000000000000000000, -0.86602540378443864676, -1.00000000000000000000, -0.86602540378443864676, -0.50000000000000000000, 0.00000000000000000000, 0.50000000000000000000, 0.86602540378443864676, 1.00000000000000000000 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; y = 0.0; fxy = 0.0; } else { x = x_vec[n_data-1]; y = y_vec[n_data-1]; fxy = f_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void arctanh_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // ARCTANH_VALUES returns some values of the hyperbolic arc tangent function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // ArcTanh[x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 June 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 15 static double fx_vec[N_MAX] = { -0.54930614433405484570, 0.00000000000000000000, 0.0010000003333335333335, 0.10033534773107558064, 0.20273255405408219099, 0.30951960420311171547, 0.42364893019360180686, 0.54930614433405484570, 0.69314718055994530942, 0.86730052769405319443, 1.0986122886681096914, 1.4722194895832202300, 2.6466524123622461977, 3.8002011672502000318, 7.2543286192620472067 }; static double x_vec[N_MAX] = { -0.5, 0.0, 0.001, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.99, 0.999, 0.999999 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bei0_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BEI0_VALUES returns some values of the Kelvin BEI function of order NU = 0. // // Discussion: // // The function is defined by: // // BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4)) // // where J(NU,X) is the J Bessel function. // // In Mathematica, BEI(NU,X) can be defined by: // // Im [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 28 June 2006 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 11 static double fx_vec[N_MAX] = { 0.0000000000000000, 0.06249321838219946, 0.2495660400366597, 0.5575600623030867, 0.9722916273066612, 1.457182044159804, 1.937586785266043, 2.283249966853915, 2.292690322699300, 1.686017203632139, 0.1160343815502004 }; static double x_vec[N_MAX] = { 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bei1_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BEI1_VALUES returns some values of the Kelvin BEI function of order NU = 1. // // Discussion: // // The function is defined by: // // BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4)) // // where J(NU,X) is the J Bessel function. // // In Mathematica, BEI(NU,X) can be defined by: // // Im [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 28 June 2006 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 11 static double fx_vec[N_MAX] = { 0.0000000000000000, 0.1711951797170153, 0.3075566313755366, 0.3678649890020899, 0.2997754370020335, 0.03866844396595048, -0.4874541770160708, -1.344042373111174, -2.563821688561078, -4.105685408400878, -5.797907901792625 }; static double x_vec[N_MAX] = { 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bell_values ( int &n_data, int &n, int &c ) //****************************************************************************80 // // Purpose: // // BELL_VALUES returns some values of the Bell numbers. // // Discussion: // // The Bell number B(N) is the number of restricted growth functions on N. // // Note that the Stirling numbers of the second kind, S^m_n, count the // number of partitions of N objects into M classes, and so it is // true that // // B(N) = S^1_N + S^2_N + ... + S^N_N. // // The Bell numbers were named for Eric Temple Bell. // // In Mathematica, the function can be evaluated by // // Sum[StirlingS2[n,m],{m,1,n}] // // Definition: // // The Bell number B(N) is defined as the number of partitions (of // any size) of a set of N distinguishable objects. // // A partition of a set is a division of the objects of the set into // subsets. // // Examples: // // There are 15 partitions of a set of 4 objects: // // (1234), // (123) (4), // (124) (3), // (12) (34), // (12) (3) (4), // (134) (2), // (13) (24), // (13) (2) (4), // (14) (23), // (1) (234), // (1) (23) (4), // (14) (2) (3), // (1) (24) (3), // (1) (2) (34), // (1) (2) (3) (4). // // and so B(4) = 15. // // First values: // // N B(N) // 0 1 // 1 1 // 2 2 // 3 5 // 4 15 // 5 52 // 6 203 // 7 877 // 8 4140 // 9 21147 // 10 115975 // // Recursion: // // B(I) = sum ( 1 <= J <=I ) Binomial ( I-1, J-1 ) * B(I-J) // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 06 February 2003 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // 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 of the Bell number. // // Output, int &C, the value of the Bell number. // { # define N_MAX 11 static int c_vec[N_MAX] = { 1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 }; static int n_vec[N_MAX] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; n = 0; c = 0; } else { n = n_vec[n_data-1]; c = c_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void ber0_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BER0_VALUES returns some values of the Kelvin BER function of order NU = 0. // // Discussion: // // The function is defined by: // // BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4)) // // where J(NU,X) is the J Bessel function. // // In Mathematica, BER(NU,X) can be defined by: // // Re [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 27 June 2006 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 11 static double fx_vec[N_MAX] = { 1.0000000000000000, 0.9990234639908383, 0.9843817812130869, 0.9210721835462558, 0.7517341827138082, 0.3999684171295313, -0.2213802495986939, -1.193598179589928, -2.563416557258580, -4.299086551599756, -6.230082478666358 }; static double x_vec[N_MAX] = { 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void ber1_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BER1_VALUES returns some values of the Kelvin BER function of order NU = 1. // // Discussion: // // The function is defined by: // // BER(NU,X) + i * BEI(NU,X) = exp(NU*Pi*I) * J(NU,X*exp(-PI*I/4)) // // where J(NU,X) is the J Bessel function. // // In Mathematica, BER(NU,X) can be defined by: // // Re [ Exp [ NU * Pi * I ] * BesselJ [ NU, X * Exp[ -Pi * I / 4 ] ] ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 28 June 2006 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 11 static double fx_vec[N_MAX] = { 0.0000000000000000, -0.1822431237551121, -0.3958682610197114, -0.6648654179597691, -0.9970776519264285, -1.373096897645111, -1.732644221128481, -1.959644131289749, -1.869248459031899, -1.202821631480086, 0.3597766667766728 }; static double x_vec[N_MAX] = { 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bernoulli_number_values ( int &n_data, int &n, double &c ) //****************************************************************************80 // // Purpose: // // BERNOULLI_NUMBER_VALUES returns some values of the Bernoulli numbers. // // Discussion: // // The Bernoulli numbers are rational. // // If we define the sum of the M-th powers of the first N integers as: // // SIGMA(M,N) = sum ( 0 <= I <= N ) I**M // // and let C(I,J) be the combinatorial coefficient: // // C(I,J) = I! / ( ( I - J )! * J! ) // // then the Bernoulli numbers B(J) satisfy: // // SIGMA(M,N) = 1/(M+1) * sum ( 0 <= J <= M ) C(M+1,J) B(J) * (N+1)**(M+1-J) // // In Mathematica, the function can be evaluated by: // // BernoulliB[n] // // First values: // // B0 1 = 1.00000000000 // B1 -1/2 = -0.50000000000 // B2 1/6 = 1.66666666666 // B3 0 = 0 // B4 -1/30 = -0.03333333333 // B5 0 = 0 // B6 1/42 = 0.02380952380 // B7 0 = 0 // B8 -1/30 = -0.03333333333 // B9 0 = 0 // B10 5/66 = 0.07575757575 // B11 0 = 0 // B12 -691/2730 = -0.25311355311 // B13 0 = 0 // B14 7/6 = 1.16666666666 // B15 0 = 0 // B16 -3617/510 = -7.09215686274 // B17 0 = 0 // B18 43867/798 = 54.97117794486 // B19 0 = 0 // B20 -174611/330 = -529.12424242424 // B21 0 = 0 // B22 854,513/138 = 6192.123 // B23 0 = 0 // B24 -236364091/2730 = -86580.257 // B25 0 = 0 // B26 8553103/6 = 1425517.16666 // B27 0 = 0 // B28 -23749461029/870 = -27298231.0678 // B29 0 = 0 // B30 8615841276005/14322 = 601580873.901 // // Recursion: // // With C(N+1,K) denoting the standard binomial coefficient, // // B(0) = 1.0 // B(N) = - ( sum ( 0 <= K < N ) C(N+1,K) * B(K) ) / C(N+1,N) // // Special Values: // // Except for B(1), all Bernoulli numbers of odd index are 0. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 09 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 of the Bernoulli number. // // Output, double &C, the value of the Bernoulli number. // { # define N_MAX 10 static double c_vec[N_MAX] = { 0.1000000000000000E+01, -0.5000000000000000E+00, 0.1666666666666667E+00, 0.0000000000000000E+00, -0.3333333333333333E-01, -0.2380952380952380E-01, -0.3333333333333333E-01, 0.7575757575757575E-01, -0.5291242424242424E+03, 0.6015808739006424E+09 }; static int n_vec[N_MAX] = { 0, 1, 2, 3, 4, 6, 8, 10, 20, 30 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; n = 0; c = 0.0E+00; } else { n = n_vec[n_data-1]; c = c_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bernoulli_poly_values ( int &n_data, int &n, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BERNOULLI_POLY_VALUES returns some values of the Bernoulli polynomials. // // Discussion: // // In Mathematica, the function can be evaluated by: // // BernoulliB[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 09 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 of the Bernoulli polynomial. // // Output, double &X, the argument of the Bernoulli polynomial. // // Output, double &FX, the value of the Bernoulli polynomial. // { # define N_MAX 27 static double fx_vec[N_MAX] = { 0.1000000000000000E+01, -0.3000000000000000E+00, 0.6666666666666667E-02, 0.4800000000000000E-01, -0.7733333333333333E-02, -0.2368000000000000E-01, 0.6913523809523810E-02, 0.2490880000000000E-01, -0.1014997333333333E-01, -0.4527820800000000E-01, 0.2332631815757576E-01, -0.3125000000000000E+00, -0.1142400000000000E+00, -0.0176800000000000E+00, 0.0156800000000000E+00, 0.0147400000000000E+00, 0.0000000000000000E+00, -0.1524000000000000E-01, -0.2368000000000000E-01, -0.2282000000000000E-01, -0.1376000000000000E-01, 0.0000000000000000E+01, 0.1376000000000000E-01, 0.2282000000000000E-01, 0.2368000000000000E-01, 0.1524000000000000E-01, 0.0000000000000000E+01 }; static int n_vec[N_MAX] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 }; static double x_vec[N_MAX] = { 0.2E+00, 0.2E+00, 0.2E+00, 0.2E+00, 0.2E+00, 0.2E+00, 0.2E+00, 0.2E+00, 0.2E+00, 0.2E+00, 0.2E+00, -0.5E+00, -0.4E+00, -0.3E+00, -0.2E+00, -0.1E+00, 0.0E+00, 0.1E+00, 0.2E+00, 0.3E+00, 0.4E+00, 0.5E+00, 0.6E+00, 0.7E+00, 0.8E+00, 0.9E+00, 1.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; n = 0; x = 0.0; fx = 0.0; } else { n = n_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bernstein_poly_01_values ( int &n_data, int &n, int &k, double &x, double &b ) //****************************************************************************80 // // Purpose: // // BERNSTEIN_POLY_01_VALUES returns some values of the Bernstein polynomials. // // Discussion: // // The Bernstein polynomials are assumed to be based on [0,1]. // // The formula for the Bernstein polynomials is // // B(N,I)(X) = [N!/(I//(N-I)!)] * (1-X)^(N-I) * X^I // // In Mathematica, the function can be evaluated by: // // Binomial[n,i] * (1-x)^(n-i) * x^i // // First values: // // B(0,0)(X) = 1 // // B(1,0)(X) = 1-X // B(1,1)(X) = X // // B(2,0)(X) = (1-X)^2 // B(2,1)(X) = 2 * (1-X) * X // B(2,2)(X) = X^2 // // B(3,0)(X) = (1-X)^3 // B(3,1)(X) = 3 * (1-X)^2 * X // B(3,2)(X) = 3 * (1-X) * X^2 // B(3,3)(X) = X^3 // // B(4,0)(X) = (1-X)^4 // B(4,1)(X) = 4 * (1-X)^3 * X // B(4,2)(X) = 6 * (1-X)^2 * X^2 // B(4,3)(X) = 4 * (1-X) * X^3 // B(4,4)(X) = X^4 // // Special values: // // B(N,I)(X) has a unique maximum value at X = I/N. // // B(N,I)(X) has an I-fold zero at 0 and and N-I fold zero at 1. // // B(N,I)(1/2) = C(N,K) / 2^N // // For a fixed X and N, the polynomials add up to 1: // // Sum ( 0 <= I <= N ) B(N,I)(X) = 1 // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 19 August 2004 // // Author: // // John Burkardt // // Reference: // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/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 degree of the polynomial. // // Output, int &K, the index of the polynomial. // // Output, double &X, the argument of the polynomial. // // Output, double &B, the value of the polynomial B(N,K)(X). // { # define N_MAX 15 static double b_vec[N_MAX] = { 0.1000000000000000E+01, 0.7500000000000000E+00, 0.2500000000000000E+00, 0.5625000000000000E+00, 0.3750000000000000E+00, 0.6250000000000000E-01, 0.4218750000000000E+00, 0.4218750000000000E+00, 0.1406250000000000E+00, 0.1562500000000000E-01, 0.3164062500000000E+00, 0.4218750000000000E+00, 0.2109375000000000E+00, 0.4687500000000000E-01, 0.3906250000000000E-02 }; static int k_vec[N_MAX] = { 0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4 }; static int n_vec[N_MAX] = { 0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4 }; static double x_vec[N_MAX] = { 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; n = 0; k = 0; x = 0.0; b = 0.0; } else { n = n_vec[n_data-1]; k = k_vec[n_data-1]; x = x_vec[n_data-1]; b = b_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_i0_int_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_I0_INT_VALUES returns some values of the Bessel I0 integral. // // Discussion: // // The function is defined by: // // I0_INT(x) = Integral ( 0 <= t <= x ) I0(t) dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 29 August 2004 // // Author: // // John Burkardt // // Reference: // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.19531256208818052282E-02, -0.39062549670565734544E-02, 0.62520348032546565850E-01, 0.12516285581366971819E+00, -0.51051480879740303760E+00, 0.10865210970235898158E+01, 0.27750019054282535299E+01, -0.13775208868039716639E+02, 0.46424372058106108576E+03, 0.64111867658021584522E+07, -0.10414860803175857953E+08, 0.44758598913855743089E+08, -0.11852985311558287888E+09, 0.31430078220715992752E+09, -0.83440212900794309620E+09, 0.22175367579074298261E+10, 0.58991731842803636487E+10, -0.41857073244691522147E+11, 0.79553885818472357663E+12, 0.15089715082719201025E+17 }; static double x_vec[N_MAX] = { 0.0019531250E+00, -0.0039062500E+00, 0.0625000000E+00, 0.1250000000E+00, -0.5000000000E+00, 1.0000000000E+00, 2.0000000000E+00, -4.0000000000E+00, 8.0000000000E+00, 18.0000000000E+00, -18.5000000000E+00, 20.0000000000E+00, -21.0000000000E+00, 22.0000000000E+00, -23.0000000000E+00, 24.0000000000E+00, 25.0000000000E+00, -27.0000000000E+00, 30.0000000000E+00, 40.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_i0_spherical_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_I0_SPHERICAL_VALUES returns some values of the Spherical Bessel function i0. // // Discussion: // // In Mathematica, the function can be evaluated by: // // Sqrt[Pi/(2*x)] * BesselI[1/2,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 06 January 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 21 static double fx_vec[N_MAX] = { 1.001667500198440E+00, 1.006680012705470E+00, 1.026880814507039E+00, 1.061089303580402E+00, 1.110132477734529E+00, 1.175201193643801E+00, 1.257884462843477E+00, 1.360215358179667E+00, 1.484729970750144E+00, 1.634541271164267E+00, 1.813430203923509E+00, 2.025956895698133E+00, 2.277595505698373E+00, 2.574897010920645E+00, 2.925685126512827E+00, 3.339291642469967E+00, 3.826838748926716E+00, 4.401577467270101E+00, 5.079293155726485E+00, 5.878791279137455E+00, 6.822479299281938E+00 }; static double x_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.2E+00, 1.4E+00, 1.6E+00, 1.8E+00, 2.0E+00, 2.2E+00, 2.4E+00, 2.6E+00, 2.8E+00, 3.0E+00, 3.2E+00, 3.4E+00, 3.6E+00, 3.8E+00, 4.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_i0_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_I0_VALUES returns some values of the I0 Bessel function. // // Discussion: // // The modified Bessel functions In(Z) and Kn(Z) are solutions of // the differential equation // // Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. // // The modified Bessel function I0(Z) corresponds to N = 0. // // In Mathematica, the function can be evaluated by: // // BesselI[0,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 20 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.1000000000000000E+01, 0.1010025027795146E+01, 0.1040401782229341E+01, 0.1092045364317340E+01, 0.1166514922869803E+01, 0.1266065877752008E+01, 0.1393725584134064E+01, 0.1553395099731217E+01, 0.1749980639738909E+01, 0.1989559356618051E+01, 0.2279585302336067E+01, 0.3289839144050123E+01, 0.4880792585865024E+01, 0.7378203432225480E+01, 0.1130192195213633E+02, 0.1748117185560928E+02, 0.2723987182360445E+02, 0.6723440697647798E+02, 0.4275641157218048E+03, 0.2815716628466254E+04 }; static double x_vec[N_MAX] = { 0.00E+00, 0.20E+00, 0.40E+00, 0.60E+00, 0.80E+00, 0.10E+01, 0.12E+01, 0.14E+01, 0.16E+01, 0.18E+01, 0.20E+01, 0.25E+01, 0.30E+01, 0.35E+01, 0.40E+01, 0.45E+01, 0.50E+01, 0.60E+01, 0.80E+01, 0.10E+02 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_i1_spherical_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_I1_SPHERICAL_VALUES returns some values of the Spherical Bessel function i1. // // Discussion: // // In Mathematica, the function can be evaluated by: // // Sqrt[Pi/(2*x)] * BesselI[3/2,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 06 January 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 21 static double fx_vec[N_MAX] = { 0.03336667857363341E+00, 0.06693371456802954E+00, 0.1354788933285401E+00, 0.2072931911031093E+00, 0.2841280857128948E+00, 0.3678794411714423E+00, 0.4606425870674146E+00, 0.5647736480096238E+00, 0.6829590627779635E+00, 0.8182955028627777E+00, 0.9743827435800610E+00, 1.155432469636406E+00, 1.366396525527973E+00, 1.613118767572064E+00, 1.902515460838681E+00, 2.242790117769266E+00, 2.643689828630357E+00, 3.116811526884873E+00, 3.675968313148932E+00, 4.337627987747642E+00, 5.121438384183637E+00 }; static double x_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.2E+00, 1.4E+00, 1.6E+00, 1.8E+00, 2.0E+00, 2.2E+00, 2.4E+00, 2.6E+00, 2.8E+00, 3.0E+00, 3.2E+00, 3.4E+00, 3.6E+00, 3.8E+00, 4.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_i1_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_I1_VALUES returns some values of the I1 Bessel function. // // Discussion: // // The modified Bessel functions In(Z) and Kn(Z) are solutions of // the differential equation // // Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. // // In Mathematica, the function can be evaluated by: // // BesselI[1,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 20 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.0000000000000000E+00, 0.1005008340281251E+00, 0.2040267557335706E+00, 0.3137040256049221E+00, 0.4328648026206398E+00, 0.5651591039924850E+00, 0.7146779415526431E+00, 0.8860919814143274E+00, 0.1084810635129880E+01, 0.1317167230391899E+01, 0.1590636854637329E+01, 0.2516716245288698E+01, 0.3953370217402609E+01, 0.6205834922258365E+01, 0.9759465153704450E+01, 0.1538922275373592E+02, 0.2433564214245053E+02, 0.6134193677764024E+02, 0.3998731367825601E+03, 0.2670988303701255E+04 }; static double x_vec[N_MAX] = { 0.00E+00, 0.20E+00, 0.40E+00, 0.60E+00, 0.80E+00, 0.10E+01, 0.12E+01, 0.14E+01, 0.16E+01, 0.18E+01, 0.20E+01, 0.25E+01, 0.30E+01, 0.35E+01, 0.40E+01, 0.45E+01, 0.50E+01, 0.60E+01, 0.80E+01, 0.10E+02 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_in_values ( int &n_data, int &nu, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_IN_VALUES returns some values of the In Bessel function. // // Discussion: // // The modified Bessel functions In(Z) and Kn(Z) are solutions of // the differential equation // // Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. // // In Mathematica, the function can be evaluated by: // // BesselI[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 21 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 &NU, the order of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 28 static double fx_vec[N_MAX] = { 0.5016687513894678E-02, 0.1357476697670383E+00, 0.6889484476987382E+00, 0.1276466147819164E+01, 0.2245212440929951E+01, 0.1750561496662424E+02, 0.2281518967726004E+04, 0.3931278522104076E+08, 0.2216842492433190E-01, 0.2127399592398527E+00, 0.1033115016915114E+02, 0.1758380716610853E+04, 0.2677764138883941E+21, 0.2714631559569719E-03, 0.9825679323131702E-02, 0.2157974547322546E+01, 0.7771882864032600E+03, 0.2278548307911282E+21, 0.2752948039836874E-09, 0.3016963879350684E-06, 0.4580044419176051E-02, 0.2189170616372337E+02, 0.1071597159477637E+21, 0.3966835985819020E-24, 0.4310560576109548E-18, 0.5024239357971806E-10, 0.1250799735644948E-03, 0.5442008402752998E+19 }; static int nu_vec[N_MAX] = { 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20 }; static double x_vec[N_MAX] = { 0.2E+00, 1.0E+00, 2.0E+00, 2.5E+00, 3.0E+00, 5.0E+00, 10.0E+00, 20.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; nu = 0; x = 0.0; fx = 0.0; } else { nu = nu_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_ix_values ( int &n_data, double &nu, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_IX_VALUES returns some values of the Ix Bessel function. // // Discussion: // // This set of data considers the less common case in which the // index of the Bessel function In is actually not an integer. // We may suggest this case by occasionally replacing the symbol // "In" by "Ix". // // The modified Bessel functions In(Z) and Kn(Z) are solutions of // the differential equation // // Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. // // In Mathematica, the function can be evaluated by: // // BesselI[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 02 March 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &NU, the order of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 28 static double fx_vec[N_MAX] = { 0.3592084175833614E+00, 0.9376748882454876E+00, 2.046236863089055E+00, 3.053093538196718E+00, 4.614822903407601E+00, 26.47754749755907E+00, 2778.784603874571E+00, 4.327974627242893E+07, 0.2935253263474798E+00, 1.099473188633110E+00, 21.18444226479414E+00, 2500.906154942118E+00, 2.866653715931464E+20, 0.05709890920304825E+00, 0.3970270801393905E+00, 13.76688213868258E+00, 2028.512757391936E+00, 2.753157630035402E+20, 0.4139416015642352E+00, 1.340196758982897E+00, 22.85715510364670E+00, 2593.006763432002E+00, 2.886630075077766E+20, 0.03590910483251082E+00, 0.2931108636266483E+00, 11.99397010023068E+00, 1894.575731562383E+00, 2.716911375760483E+20 }; static double nu_vec[N_MAX] = { 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 1.25E+00, 1.25E+00, 1.25E+00, 1.25E+00, 1.25E+00, 2.75E+00, 2.75E+00, 2.75E+00, 2.75E+00, 2.75E+00 }; static double x_vec[N_MAX] = { 0.2E+00, 1.0E+00, 2.0E+00, 2.5E+00, 3.0E+00, 5.0E+00, 10.0E+00, 20.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; nu = 0.0; x = 0.0; fx = 0.0; } else { nu = nu_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_j0_int_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_J0_INT_VALUES returns some values of the Bessel J0 integral. // // Discussion: // // The function is defined by: // // J0_INT(x) = Integral ( 0 <= t <= x ) J0(t) dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 29 August 2004 // // Author: // // John Burkardt // // Reference: // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.97656242238978822427E-03, 0.39062450329491108875E-02, -0.62479657927917933620E-01, 0.12483733492120479139E+00, -0.48968050664604505505E+00, 0.91973041008976023931E+00, -0.14257702931970265690E+01, 0.10247341594606064818E+01, -0.12107468348304501655E+01, 0.11008652032736190799E+01, -0.10060334829904124192E+01, 0.81330572662485953519E+00, -0.10583788214211277585E+01, 0.87101492116545875169E+00, -0.88424908882547488420E+00, 0.11257761503599914603E+01, -0.90141212258183461184E+00, 0.91441344369647797803E+00, -0.94482281938334394886E+00, 0.92266255696016607257E+00 }; static double x_vec[N_MAX] = { 0.0009765625E+00, 0.0039062500E+00, -0.0625000000E+00, 0.1250000000E+00, -0.5000000000E+00, 1.0000000000E+00, -2.0000000000E+00, 4.0000000000E+00, -8.0000000000E+00, 16.0000000000E+00, -16.5000000000E+00, 18.0000000000E+00, -20.0000000000E+00, 25.0000000000E+00, -30.0000000000E+00, 40.0000000000E+00, -50.0000000000E+00, 75.0000000000E+00, -80.0000000000E+00, 100.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_j0_spherical_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_J0_SPHERICAL_VALUES returns some values of the Spherical Bessel function j0. // // Discussion: // // In Mathematica, the function can be evaluated by: // // Sqrt[Pi/(2*x)] * BesselJ[1/2,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 21 static double fx_vec[N_MAX] = { 0.9983341664682815E+00, 0.9933466539753061E+00, 0.9735458557716262E+00, 0.9410707889917256E+00, 0.8966951136244035E+00, 0.8414709848078965E+00, 0.7766992383060220E+00, 0.7038926642774716E+00, 0.6247335019009407E+00, 0.5410264615989973E+00, 0.4546487134128408E+00, 0.3674983653725410E+00, 0.2814429918963129E+00, 0.1982697583928709E+00, 0.1196386250556803E+00, 0.4704000268662241E-01, -0.1824191982111872E-01, -0.7515914765495039E-01, -0.1229223453596812E+00, -0.1610152344586103E+00, -0.1892006238269821E+00 }; static double x_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.2E+00, 1.4E+00, 1.6E+00, 1.8E+00, 2.0E+00, 2.2E+00, 2.4E+00, 2.6E+00, 2.8E+00, 3.0E+00, 3.2E+00, 3.4E+00, 3.6E+00, 3.8E+00, 4.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_j0_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_J0_VALUES returns some values of the J0 Bessel function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // BesselJ[0,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 10 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 21 static double fx_vec[N_MAX] = { -0.1775967713143383E+00, -0.3971498098638474E+00, -0.2600519549019334E+00, 0.2238907791412357E+00, 0.7651976865579666E+00, 0.1000000000000000E+01, 0.7651976865579666E+00, 0.2238907791412357E+00, -0.2600519549019334E+00, -0.3971498098638474E+00, -0.1775967713143383E+00, 0.1506452572509969E+00, 0.3000792705195556E+00, 0.1716508071375539E+00, -0.9033361118287613E-01, -0.2459357644513483E+00, -0.1711903004071961E+00, 0.4768931079683354E-01, 0.2069261023770678E+00, 0.1710734761104587E+00, -0.1422447282678077E-01 }; static double x_vec[N_MAX] = { -5.0E+00, -4.0E+00, -3.0E+00, -2.0E+00, -1.0E+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, 11.0E+00, 12.0E+00, 13.0E+00, 14.0E+00, 15.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_j1_spherical_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_J1_SPHERICAL_VALUES returns some values of the Spherical Bessel function j1. // // Discussion: // // In Mathematica, the function can be evaluated by: // // Sqrt[Pi/(2*x)] * BesselJ[3/2,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 23 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 21 static double fx_vec[N_MAX] = { 0.3330001190255757E-01, 0.6640038067032223E-01, 0.1312121544218529E+00, 0.1928919568034122E+00, 0.2499855053465475E+00, 0.3011686789397568E+00, 0.3452845698577903E+00, 0.3813753724123076E+00, 0.4087081401263934E+00, 0.4267936423844913E+00, 0.4353977749799916E+00, 0.4345452193763121E+00, 0.4245152947656493E+00, 0.4058301968314685E+00, 0.3792360591872637E+00, 0.3456774997623560E+00, 0.3062665174917607E+00, 0.2622467779189737E+00, 0.2149544641595738E+00, 0.1657769677515280E+00, 0.1161107492591575E+00 }; static double x_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.2E+00, 1.4E+00, 1.6E+00, 1.8E+00, 2.0E+00, 2.2E+00, 2.4E+00, 2.6E+00, 2.8E+00, 3.0E+00, 3.2E+00, 3.4E+00, 3.6E+00, 3.8E+00, 4.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_j1_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_J1_VALUES returns some values of the J1 Bessel function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // BesselJ[1,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 12 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 21 static double fx_vec[N_MAX] = { 0.3275791375914652E+00, 0.6604332802354914E-01, -0.3390589585259365E+00, -0.5767248077568734E+00, -0.4400505857449335E+00, 0.0000000000000000E+00, 0.4400505857449335E+00, 0.5767248077568734E+00, 0.3390589585259365E+00, -0.6604332802354914E-01, -0.3275791375914652E+00, -0.2766838581275656E+00, -0.4682823482345833E-02, 0.2346363468539146E+00, 0.2453117865733253E+00, 0.4347274616886144E-01, -0.1767852989567215E+00, -0.2234471044906276E+00, -0.7031805212177837E-01, 0.1333751546987933E+00, 0.2051040386135228E+00 }; static double x_vec[N_MAX] = { -5.0E+00, -4.0E+00, -3.0E+00, -2.0E+00, -1.0E+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, 11.0E+00, 12.0E+00, 13.0E+00, 14.0E+00, 15.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_jn_values ( int &n_data, int &nu, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_JN_VALUES returns some values of the Jn Bessel function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // BesselJ[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 29 April 2001 // // 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 &NU, the order of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.1149034849319005E+00, 0.3528340286156377E+00, 0.4656511627775222E-01, 0.2546303136851206E+00, -0.5971280079425882E-01, 0.2497577302112344E-03, 0.7039629755871685E-02, 0.2611405461201701E+00, -0.2340615281867936E+00, -0.8140024769656964E-01, 0.2630615123687453E-09, 0.2515386282716737E-06, 0.1467802647310474E-02, 0.2074861066333589E+00, -0.1138478491494694E+00, 0.3873503008524658E-24, 0.3918972805090754E-18, 0.2770330052128942E-10, 0.1151336924781340E-04, -0.1167043527595797E+00 }; static int nu_vec[N_MAX] = { 2, 2, 2, 2, 2, 5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20 }; static double x_vec[N_MAX] = { 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; nu = 0; x = 0.0; fx = 0.0; } else { nu = nu_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_jx_values ( int &n_data, double &nu, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_JX_VALUES returns some values of the Jx Bessel function. // // Discussion: // // This set of data considers the less common case in which the // index of the Bessel function Jn is actually not an integer. // We may suggest this case by occasionally replacing the symbol // "Jn" by "Jx". // // In Mathematica, the function can be evaluated by: // // BesselJ[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 01 April 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &NU, the order of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 28 static double fx_vec[N_MAX] = { 0.3544507442114011E+00, 0.6713967071418031E+00, 0.5130161365618278E+00, 0.3020049060623657E+00, 0.06500818287737578E+00, -0.3421679847981618E+00, -0.1372637357550505E+00, 0.1628807638550299E+00, 0.2402978391234270E+00, 0.4912937786871623E+00, -0.1696513061447408E+00, 0.1979824927558931E+00, -0.1094768729883180E+00, 0.04949681022847794E+00, 0.2239245314689158E+00, 0.2403772011113174E+00, 0.1966584835818184E+00, 0.02303721950962553E+00, 0.3314145508558904E+00, 0.5461734240402840E+00, -0.2616584152094124E+00, 0.1296035513791289E+00, -0.1117432171933552E+00, 0.03142623570527935E+00, 0.1717922192746527E+00, 0.3126634069544786E+00, 0.1340289119304364E+00, 0.06235967135106445E+00 }; static double nu_vec[N_MAX] = { 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 1.25E+00, 1.25E+00, 1.25E+00, 1.25E+00, 1.25E+00, 2.75E+00, 2.75E+00, 2.75E+00, 2.75E+00, 2.75E+00 }; static double x_vec[N_MAX] = { 0.2E+00, 1.0E+00, 2.0E+00, 2.5E+00, 3.0E+00, 5.0E+00, 10.0E+00, 20.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; nu = 0.0; x = 0.0; fx = 0.0; } else { nu = nu_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_k0_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_K0_VALUES returns some values of the K0 Bessel function. // // Discussion: // // The modified Bessel functions In(Z) and Kn(Z) are solutions of // the differential equation // // Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. // // The modified Bessel function K0(Z) corresponds to N = 0. // // In Mathematica, the function can be evaluated by: // // BesselK[0,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 15 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.2427069024702017E+01, 0.1752703855528146E+01, 0.1114529134524434E+01, 0.7775220919047293E+00, 0.5653471052658957E+00, 0.4210244382407083E+00, 0.3185082202865936E+00, 0.2436550611815419E+00, 0.1879547519693323E+00, 0.1459314004898280E+00, 0.1138938727495334E+00, 0.6234755320036619E-01, 0.3473950438627925E-01, 0.1959889717036849E-01, 0.1115967608585302E-01, 0.6399857243233975E-02, 0.3691098334042594E-02, 0.1243994328013123E-02, 0.1464707052228154E-03, 0.1778006231616765E-04 }; static double x_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.2E+00, 1.4E+00, 1.6E+00, 1.8E+00, 2.0E+00, 2.5E+00, 3.0E+00, 3.5E+00, 4.0E+00, 4.5E+00, 5.0E+00, 6.0E+00, 8.0E+00, 10.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_k0_int_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_K0_INT_VALUES returns some values of the Bessel K0 integral. // // Discussion: // // The function is defined by: // // K0_INT(x) = Integral ( 0 <= t <= x ) K0(t) dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 29 August 2004 // // Author: // // John Burkardt // // Reference: // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.78587929563466784589E-02, 0.26019991617330578111E-01, 0.24311842237541167904E+00, 0.39999633750480508861E+00, 0.92710252093114907345E+00, 0.12425098486237782662E+01, 0.14736757343168286825E+01, 0.15606495706051741364E+01, 0.15673873907283660493E+01, 0.15696345532693743714E+01, 0.15701153443250786355E+01, 0.15706574852894436220E+01, 0.15707793116159788598E+01, 0.15707942066465767196E+01, 0.15707962315469192247E+01, 0.15707963262340149876E+01, 0.15707963267948756308E+01, 0.15707963267948966192E+01, 0.15707963267948966192E+01, 0.15707963267948966192E+01 }; static double x_vec[N_MAX] = { 0.0009765625E+00, 0.0039062500E+00, 0.0625000000E+00, 0.1250000000E+00, 0.5000000000E+00, 1.0000000000E+00, 2.0000000000E+00, 4.0000000000E+00, 5.0000000000E+00, 6.0000000000E+00, 6.5000000000E+00, 8.0000000000E+00, 10.0000000000E+00, 12.0000000000E+00, 15.0000000000E+00, 20.0000000000E+00, 30.0000000000E+00, 50.0000000000E+00, 80.0000000000E+00, 100.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_k1_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_K1_VALUES returns some values of the K1 Bessel function. // // Discussion: // // The modified Bessel functions In(Z) and Kn(Z) are solutions of // the differential equation // // Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. // // The modified Bessel function K1(Z) corresponds to N = 1. // // In Mathematica, the function can be evaluated by: // // BesselK[1,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 15 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { 0.9853844780870606E+01, 0.4775972543220472E+01, 0.2184354424732687E+01, 0.1302834939763502E+01, 0.8617816344721803E+00, 0.6019072301972346E+00, 0.4345923910607150E+00, 0.3208359022298758E+00, 0.2406339113576119E+00, 0.1826230998017470E+00, 0.1398658818165224E+00, 0.7389081634774706E-01, 0.4015643112819418E-01, 0.2223939292592383E-01, 0.1248349888726843E-01, 0.7078094908968090E-02, 0.4044613445452164E-02, 0.1343919717735509E-02, 0.1553692118050011E-03, 0.1864877345382558E-04 }; static double x_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.2E+00, 1.4E+00, 1.6E+00, 1.8E+00, 2.0E+00, 2.5E+00, 3.0E+00, 3.5E+00, 4.0E+00, 4.5E+00, 5.0E+00, 6.0E+00, 8.0E+00, 10.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_kn_values ( int &n_data, int &nu, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_KN_VALUES returns some values of the Kn Bessel function. // // Discussion: // // The modified Bessel functions In(Z) and Kn(Z) are solutions of // the differential equation // // Z^2 * W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. // // In Mathematica, the function can be evaluated by: // // BesselK[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 24 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 &NU, the order of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 28 static double fx_vec[N_MAX] = { 0.4951242928773287E+02, 0.1624838898635177E+01, 0.2537597545660559E+00, 0.1214602062785638E+00, 0.6151045847174204E-01, 0.5308943712223460E-02, 0.2150981700693277E-04, 0.6329543612292228E-09, 0.7101262824737945E+01, 0.6473853909486342E+00, 0.8291768415230932E-02, 0.2725270025659869E-04, 0.3727936773826211E-22, 0.3609605896012407E+03, 0.9431049100596467E+01, 0.3270627371203186E-01, 0.5754184998531228E-04, 0.4367182254100986E-22, 0.1807132899010295E+09, 0.1624824039795591E+06, 0.9758562829177810E+01, 0.1614255300390670E-02, 0.9150988209987996E-22, 0.6294369360424535E+23, 0.5770856852700241E+17, 0.4827000520621485E+09, 0.1787442782077055E+03, 0.1706148379722035E-20 }; static int nu_vec[N_MAX] = { 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20 }; static double x_vec[N_MAX] = { 0.2E+00, 1.0E+00, 2.0E+00, 2.5E+00, 3.0E+00, 5.0E+00, 10.0E+00, 20.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; nu = 0; x = 0.0; fx = 0.0; } else { nu = nu_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_kx_values ( int &n_data, double &nu, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_KX_VALUES returns some values of the Kx Bessel function. // // Discussion: // // This set of data considers the less common case in which the // index of the Bessel function Kn is actually not an integer. // We may suggest this case by occasionally replacing the symbol // "Kn" by "Kx". // // The modified Bessel functions In(Z) and Kn(Z) are solutions of // the differential equation // // Z^2 W'' + Z * W' - ( Z^2 + N^2 ) * W = 0. // // In Mathematica, the function can be evaluated by: // // BesselK[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 01 April 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &NU, the order of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 28 static double fx_vec[N_MAX] = { 2.294489339798475E+00, 0.4610685044478946E+00, 0.1199377719680614E+00, 0.06506594315400999E+00, 0.03602598513176459E+00, 0.003776613374642883E+00, 0.00001799347809370518E+00, 5.776373974707445E-10, 0.9221370088957891E+00, 0.1799066579520922E+00, 0.004531936049571459E+00, 0.00001979282590307570E+00, 3.486992497366216E-23, 3.227479531135262E+00, 0.3897977588961997E+00, 0.006495775004385758E+00, 0.00002393132586462789E+00, 3.627839645299048E-23, 0.7311451879202114E+00, 0.1567475478393932E+00, 0.004257389528177461E+00, 0.00001915541065869563E+00, 3.463337593569306E-23, 4.731184839919541E+00, 0.4976876225514758E+00, 0.007300864610941163E+00, 0.00002546421294106458E+00, 3.675275677913656E-23 }; static double nu_vec[N_MAX] = { 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 1.25E+00, 1.25E+00, 1.25E+00, 1.25E+00, 1.25E+00, 2.75E+00, 2.75E+00, 2.75E+00, 2.75E+00, 2.75E+00 }; static double x_vec[N_MAX] = { 0.2E+00, 1.0E+00, 2.0E+00, 2.5E+00, 3.0E+00, 5.0E+00, 10.0E+00, 20.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; nu = 0.0; x = 0.0; fx = 0.0; } else { nu = nu_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_y0_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_Y0_VALUES returns some values of the Y0 Bessel function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // BesselY[0,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 16 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 16 static double fx_vec[N_MAX] = { -0.1534238651350367E+01, 0.8825696421567696E-01, 0.5103756726497451E+00, 0.3768500100127904E+00, -0.1694073932506499E-01, -0.3085176252490338E+00, -0.2881946839815792E+00, -0.2594974396720926E-01, 0.2235214893875662E+00, 0.2499366982850247E+00, 0.5567116728359939E-01, -0.1688473238920795E+00, -0.2252373126343614E+00, -0.7820786452787591E-01, 0.1271925685821837E+00, 0.2054642960389183E+00 }; static double x_vec[N_MAX] = { 0.1E+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, 11.0E+00, 12.0E+00, 13.0E+00, 14.0E+00, 15.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_y0_int_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_Y0_INT_VALUES returns some values of the Bessel Y0 integral. // // Discussion: // // The function is defined by: // // Y0_INT(x) = Integral ( 0 <= t <= x ) Y0(t) dt // // The data was reported by McLeod. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 30 August 2004 // // Author: // // John Burkardt // // Reference: // // Allan McLeod, // Algorithm 757: // MISCFUN: A software package to compute uncommon special functions, // ACM Transactions on Mathematical Software, // Volume 22, Number 3, September 1996, pages 288-301. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 20 static double fx_vec[N_MAX] = { -0.91442642860172110926E-02, -0.29682047390397591290E-01, -0.25391431276585388961E+00, -0.56179545591464028187E+00, -0.63706937660742309754E+00, -0.28219285008510084123E+00, 0.38366964785312561103E+00, -0.12595061285798929390E+00, 0.24129031832266684828E+00, 0.17138069757627037938E+00, 0.18958142627134083732E+00, 0.17203846136449706946E+00, -0.16821597677215029611E+00, -0.93607927351428988679E-01, 0.88229711948036648408E-01, -0.89324662736274161841E-02, -0.54814071000063488284E-01, -0.94958246003466381588E-01, -0.19598064853404969850E-01, -0.83084772357154773468E-02 }; static double x_vec[N_MAX] = { 0.0019531250E+00, 0.0078125000E+00, 0.1250000000E+00, 0.5000000000E+00, 1.0000000000E+00, 2.0000000000E+00, 4.0000000000E+00, 6.0000000000E+00, 10.0000000000E+00, 16.0000000000E+00, 16.2500000000E+00, 17.0000000000E+00, 20.0000000000E+00, 25.0000000000E+00, 30.0000000000E+00, 40.0000000000E+00, 50.0000000000E+00, 70.0000000000E+00, 100.0000000000E+00, 125.0000000000E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_y0_spherical_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_Y0_SPHERICAL_VALUES returns some values of the Spherical Bessel function y0. // // Discussion: // // In Mathematica, the function can be evaluated by: // // Sqrt[Pi/(2*x)] * BesselY[1/2,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 24 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 21 static double fx_vec[N_MAX] = { -0.9950041652780258E+01, -0.4900332889206208E+01, -0.2302652485007213E+01, -0.1375559358182797E+01, -0.8708833866839568E+00, -0.5403023058681397E+00, -0.3019647953972280E+00, -0.1214051020716007E+00, 0.1824970143830545E-01, 0.1262233859406039E+00, 0.2080734182735712E+00, 0.2675005078433390E+00, 0.3072473814755190E+00, 0.3295725974495951E+00, 0.3365079788102351E+00, 0.3299974988668152E+00, 0.3119671174358603E+00, 0.2843524095821944E+00, 0.2490995600928186E+00, 0.2081493978722149E+00, 0.1634109052159030E+00 }; static double x_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.2E+00, 1.4E+00, 1.6E+00, 1.8E+00, 2.0E+00, 2.2E+00, 2.4E+00, 2.6E+00, 2.8E+00, 3.0E+00, 3.2E+00, 3.4E+00, 3.6E+00, 3.8E+00, 4.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_y1_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_Y1_VALUES returns some values of the Y1 Bessel function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // BesselY[1,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 16 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 16 static double fx_vec[N_MAX] = { -0.6458951094702027E+01, -0.7812128213002887E+00, -0.1070324315409375E+00, 0.3246744247918000E+00, 0.3979257105571000E+00, 0.1478631433912268E+00, -0.1750103443003983E+00, -0.3026672370241849E+00, -0.1580604617312475E+00, 0.1043145751967159E+00, 0.2490154242069539E+00, 0.1637055374149429E+00, -0.5709921826089652E-01, -0.2100814084206935E+00, -0.1666448418561723E+00, 0.2107362803687351E-01 }; static double x_vec[N_MAX] = { 0.1E+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, 11.0E+00, 12.0E+00, 13.0E+00, 14.0E+00, 15.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_y1_spherical_values ( int &n_data, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_Y1_SPHERICAL_VALUES returns some values of the Spherical Bessel function y1. // // Discussion: // // In Mathematica, the function can be evaluated by: // // Sqrt[Pi/(2*x)] * BesselY[3/2,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 22 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, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 21 static double fx_vec[N_MAX] = { -0.1004987506942709E+03, -0.2549501110000635E+02, -0.6730177068289658E+01, -0.3233669719296388E+01, -0.1985299346979349E+01, -0.1381773290676036E+01, -0.1028336567803712E+01, -0.7906105943286149E+00, -0.6133274385019998E+00, -0.4709023582986618E+00, -0.3506120042760553E+00, -0.2459072254437506E+00, -0.1534232496148467E+00, -0.7151106706610352E-01, 0.5427959479750482E-03, 0.6295916360231598E-01, 0.1157316440198251E+00, 0.1587922092967723E+00, 0.1921166676076864E+00, 0.2157913917934037E+00, 0.2300533501309578E+00 }; static double x_vec[N_MAX] = { 0.1E+00, 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.2E+00, 1.4E+00, 1.6E+00, 1.8E+00, 2.0E+00, 2.2E+00, 2.4E+00, 2.6E+00, 2.8E+00, 3.0E+00, 3.2E+00, 3.4E+00, 3.6E+00, 3.8E+00, 4.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; fx = 0.0; } else { x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_yn_values ( int &n_data, int &nu, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_YN_VALUES returns some values of the Yn Bessel function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // BesselY[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 26 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 &NU, the order of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { // # define N_MAX 20 static double fx_vec[N_MAX] = { -0.1650682606816254E+01, -0.6174081041906827E+00, 0.3676628826055245E+00, -0.5868082442208615E-02, 0.9579316872759649E-01, -0.2604058666258122E+03, -0.9935989128481975E+01, -0.4536948224911019E+00, 0.1354030476893623E+00, -0.7854841391308165E-01, -0.1216180142786892E+09, -0.1291845422080393E+06, -0.2512911009561010E+02, -0.3598141521834027E+00, 0.5723897182053514E-02, -0.4113970314835505E+23, -0.4081651388998367E+17, -0.5933965296914321E+09, -0.1597483848269626E+04, 0.1644263394811578E-01 }; static int nu_vec[N_MAX] = { 2, 2, 2, 2, 2, 5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 20, 20, 20, 20, 20 }; static double x_vec[N_MAX] = { 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; nu = 0; x = 0.0; fx = 0.0; } else { nu = nu_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bessel_yx_values ( int &n_data, double &nu, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BESSEL_YX_VALUES returns some values of the Yx Bessel function. // // Discussion: // // This set of data considers the less common case in which the // index of the Bessel function Yn is actually not an integer. // We may suggest this case by occasionally replacing the symbol // "Yn" by "Yx". // // In Mathematica, the function can be evaluated by: // // BesselY[n,x] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 01 April 2007 // // Author: // // John Burkardt // // Reference: // // Milton Abramowitz, Irene Stegun, // Handbook of Mathematical Functions, // National Bureau of Standards, 1964, // ISBN: 0-486-61272-4, // LC: QA47.A34. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &NU, the order of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 28 static double fx_vec[N_MAX] = { -1.748560416961876E+00, -0.4310988680183761E+00, 0.2347857104062485E+00, 0.4042783022390569E+00, 0.4560488207946332E+00, -0.1012177091851084E+00, 0.2117088663313982E+00, -0.07280690478506185E+00, -1.102495575160179E+00, -0.3956232813587035E+00, 0.3219244429611401E+00, 0.1584346223881903E+00, 0.02742813676191382E+00, -2.876387857462161E+00, -0.8282206324443037E+00, 0.2943723749617925E+00, -0.1641784796149411E+00, 0.1105304445562544E+00, -0.9319659251969881E+00, -0.2609445010948933E+00, 0.2492796362185881E+00, 0.2174410301416733E+00, -0.01578576650557229E+00, -4.023453301501028E+00, -0.9588998694752389E+00, 0.2264260361047367E+00, -0.2193617736566760E+00, 0.09413988344515077E+00 }; static double nu_vec[N_MAX] = { 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 0.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 1.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 2.50E+00, 1.25E+00, 1.25E+00, 1.25E+00, 1.25E+00, 1.25E+00, 2.75E+00, 2.75E+00, 2.75E+00, 2.75E+00, 2.75E+00 }; static double x_vec[N_MAX] = { 0.2E+00, 1.0E+00, 2.0E+00, 2.5E+00, 3.0E+00, 5.0E+00, 10.0E+00, 20.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00, 1.0E+00, 2.0E+00, 5.0E+00, 10.0E+00, 50.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; nu = 0.0; x = 0.0; fx = 0.0; } else { nu = nu_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void beta_cdf_values ( int &n_data, double &a, double &b, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BETA_CDF_VALUES returns some values of the Beta CDF. // // Discussion: // // The incomplete Beta function may be written // // BETA_INC(A,B,X) = Integral (0 <= t <= X) T^(A-1) * (1-T)^(B-1) dT // / Integral (0 <= t <= 1) T^(A-1) * (1-T)^(B-1) dT // // Thus, // // BETA_INC(A,B,0.0) = 0.0; // BETA_INC(A,B,1.0) = 1.0 // // The incomplete Beta function is also sometimes called the // "modified" Beta function, or the "normalized" Beta function // or the Beta CDF (cumulative density function. // // In Mathematica, the function can be evaluated by: // // BETA[X,A,B] / BETA[A,B] // // The function can also be evaluated by using the Statistics package: // // Needs["Statistics`ContinuousDistributions`"] // dist = BetaDistribution [ a, b ] // CDF [ dist, x ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 28 April 2013 // // 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. // // Karl Pearson, // Tables of the Incomplete Beta Function, // Cambridge University Press, 1968. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &A, &B, the parameters of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 45 static double a_vec[N_MAX] = { 0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 5.5E+00, 10.0E+00, 10.0E+00, 10.0E+00, 10.0E+00, 20.0E+00, 20.0E+00, 20.0E+00, 20.0E+00, 20.0E+00, 30.0E+00, 30.0E+00, 40.0E+00, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.2E+01, 0.3E+01, 0.4E+01, 0.5E+01, 1.30625, 1.30625, 1.30625 }; static double b_vec[N_MAX] = { 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 5.0E+00, 0.5E+00, 5.0E+00, 5.0E+00, 10.0E+00, 5.0E+00, 10.0E+00, 10.0E+00, 20.0E+00, 20.0E+00, 10.0E+00, 10.0E+00, 20.0E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.2E+01, 0.3E+01, 0.4E+01, 0.5E+01, 0.2E+01, 0.2E+01, 0.2E+01, 0.2E+01, 11.7562, 11.7562, 11.7562 }; static double fx_vec[N_MAX] = { 0.6376856085851985E-01, 0.2048327646991335E+00, 0.1000000000000000E+01, 0.0000000000000000E+00, 0.5012562893380045E-02, 0.5131670194948620E-01, 0.2928932188134525E+00, 0.5000000000000000E+00, 0.2800000000000000E-01, 0.1040000000000000E+00, 0.2160000000000000E+00, 0.3520000000000000E+00, 0.5000000000000000E+00, 0.6480000000000000E+00, 0.7840000000000000E+00, 0.8960000000000000E+00, 0.9720000000000000E+00, 0.4361908850559777E+00, 0.1516409096347099E+00, 0.8978271484375000E-01, 0.1000000000000000E+01, 0.5000000000000000E+00, 0.4598773297575791E+00, 0.2146816102371739E+00, 0.9507364826957875E+00, 0.5000000000000000E+00, 0.8979413687105918E+00, 0.2241297491808366E+00, 0.7586405487192086E+00, 0.7001783247477069E+00, 0.5131670194948620E-01, 0.1055728090000841E+00, 0.1633399734659245E+00, 0.2254033307585166E+00, 0.3600000000000000E+00, 0.4880000000000000E+00, 0.5904000000000000E+00, 0.6723200000000000E+00, 0.2160000000000000E+00, 0.8370000000000000E-01, 0.3078000000000000E-01, 0.1093500000000000E-01, 0.918884684620518, 0.21052977489419, 0.1824130512500673 }; static double x_vec[N_MAX] = { 0.01E+00, 0.10E+00, 1.00E+00, 0.00E+00, 0.01E+00, 0.10E+00, 0.50E+00, 0.50E+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, 0.50E+00, 0.90E+00, 0.50E+00, 1.00E+00, 0.50E+00, 0.80E+00, 0.60E+00, 0.80E+00, 0.50E+00, 0.60E+00, 0.70E+00, 0.80E+00, 0.70E+00, 0.10E+00, 0.20E+00, 0.30E+00, 0.40E+00, 0.20E+00, 0.20E+00, 0.20E+00, 0.20E+00, 0.30E+00, 0.30E+00, 0.30E+00, 0.30E+00, 0.225609, 0.0335568, 0.0295222 }; 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; x = 0.0; fx = 0.0; } else { a = a_vec[n_data-1]; b = b_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void beta_inc_values ( int &n_data, double &a, double &b, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BETA_INC_VALUES returns some values of the incomplete Beta function. // // Discussion: // // The incomplete Beta function may be written // // BETA_INC(A,B,X) = Integral (0 <= t <= X) T^(A-1) * (1-T)^(B-1) dT // / Integral (0 <= t <= 1) T^(A-1) * (1-T)^(B-1) dT // // Thus, // // BETA_INC(A,B,0.0) = 0.0; // BETA_INC(A,B,1.0) = 1.0 // // The incomplete Beta function is also sometimes called the // "modified" Beta function, or the "normalized" Beta function // or the Beta CDF (cumulative density function. // // In Mathematica, the function can be evaluated by: // // BETA[X,A,B] / BETA[A,B] // // The function can also be evaluated by using the Statistics package: // // Needs["Statistics`ContinuousDistributions`"] // dist = BetaDistribution [ a, b ] // CDF [ dist, x ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 28 April 2013 // // 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. // // Karl Pearson, // Tables of the Incomplete Beta Function, // Cambridge University Press, 1968. // // Stephen Wolfram, // The Mathematica Book, // Fourth Edition, // Cambridge University Press, 1999, // ISBN: 0-521-64314-7, // LC: QA76.95.W65. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &A, &B, the parameters of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 45 static double a_vec[N_MAX] = { 0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 5.5E+00, 10.0E+00, 10.0E+00, 10.0E+00, 10.0E+00, 20.0E+00, 20.0E+00, 20.0E+00, 20.0E+00, 20.0E+00, 30.0E+00, 30.0E+00, 40.0E+00, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.2E+01, 0.3E+01, 0.4E+01, 0.5E+01, 1.30625, 1.30625, 1.30625 }; static double b_vec[N_MAX] = { 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 1.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 2.0E+00, 5.0E+00, 0.5E+00, 5.0E+00, 5.0E+00, 10.0E+00, 5.0E+00, 10.0E+00, 10.0E+00, 20.0E+00, 20.0E+00, 10.0E+00, 10.0E+00, 20.0E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.5E+00, 0.2E+01, 0.3E+01, 0.4E+01, 0.5E+01, 0.2E+01, 0.2E+01, 0.2E+01, 0.2E+01, 11.7562, 11.7562, 11.7562 }; static double fx_vec[N_MAX] = { 0.6376856085851985E-01, 0.2048327646991335E+00, 0.1000000000000000E+01, 0.0000000000000000E+00, 0.5012562893380045E-02, 0.5131670194948620E-01, 0.2928932188134525E+00, 0.5000000000000000E+00, 0.2800000000000000E-01, 0.1040000000000000E+00, 0.2160000000000000E+00, 0.3520000000000000E+00, 0.5000000000000000E+00, 0.6480000000000000E+00, 0.7840000000000000E+00, 0.8960000000000000E+00, 0.9720000000000000E+00, 0.4361908850559777E+00, 0.1516409096347099E+00, 0.8978271484375000E-01, 0.1000000000000000E+01, 0.5000000000000000E+00, 0.4598773297575791E+00, 0.2146816102371739E+00, 0.9507364826957875E+00, 0.5000000000000000E+00, 0.8979413687105918E+00, 0.2241297491808366E+00, 0.7586405487192086E+00, 0.7001783247477069E+00, 0.5131670194948620E-01, 0.1055728090000841E+00, 0.1633399734659245E+00, 0.2254033307585166E+00, 0.3600000000000000E+00, 0.4880000000000000E+00, 0.5904000000000000E+00, 0.6723200000000000E+00, 0.2160000000000000E+00, 0.8370000000000000E-01, 0.3078000000000000E-01, 0.1093500000000000E-01, 0.918884684620518, 0.21052977489419, 0.1824130512500673 }; static double x_vec[N_MAX] = { 0.01E+00, 0.10E+00, 1.00E+00, 0.00E+00, 0.01E+00, 0.10E+00, 0.50E+00, 0.50E+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, 0.50E+00, 0.90E+00, 0.50E+00, 1.00E+00, 0.50E+00, 0.80E+00, 0.60E+00, 0.80E+00, 0.50E+00, 0.60E+00, 0.70E+00, 0.80E+00, 0.70E+00, 0.10E+00, 0.20E+00, 0.30E+00, 0.40E+00, 0.20E+00, 0.20E+00, 0.20E+00, 0.20E+00, 0.30E+00, 0.30E+00, 0.30E+00, 0.30E+00, 0.225609, 0.0335568, 0.0295222 }; 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; x = 0.0; fx = 0.0; } else { a = a_vec[n_data-1]; b = b_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void beta_log_values ( int &n_data, double &x, double &y, double &fxy ) //****************************************************************************80 // // Purpose: // // BETA_LOG_VALUES returns some values of the logarithm of the Beta function. // // Discussion: // // In Mathematica, the function can be evaluated by: // // Log[Beta[x]] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 14 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, double &X, &Y, the arguments of the function. // // Output, double &FXY, the value of the function. // { # define N_MAX 17 static double fxy_vec[N_MAX] = { 0.1609437912434100E+01, 0.9162907318741551E+00, 0.5108256237659907E+00, 0.2231435513142098E+00, 0.1609437912434100E+01, 0.9162907318741551E+00, 0.0000000000000000E+00, -0.1791759469228055E+01, -0.3401197381662155E+01, -0.4941642422609304E+01, -0.6445719819385578E+01, -0.3737669618283368E+01, -0.5123963979403259E+01, -0.6222576268071369E+01, -0.7138866999945524E+01, -0.7927324360309794E+01, -0.9393661429103221E+01 }; static double x_vec[N_MAX] = { 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.0E+00, 1.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, 6.0E+00, 6.0E+00, 6.0E+00, 6.0E+00, 7.0E+00 }; static double y_vec[N_MAX] = { 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 0.2E+00, 0.4E+00, 1.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, 7.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; y = 0.0; fxy = 0.0; } else { x = x_vec[n_data-1]; y = y_vec[n_data-1]; fxy = fxy_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void beta_noncentral_cdf_values ( int &n_data, double &a, double &b, double &lambda, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // BETA_NONCENTRAL_CDF_VALUES returns some values of the noncentral Beta CDF. // // Discussion: // // The values presented here are taken from the reference, where they // were given to a limited number of decimal places. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 24 January 2008 // // Author: // // John Burkardt // // Reference: // // R Chattamvelli, R Shanmugam, // Algorithm AS 310: // Computing the Non-central Beta Distribution Function, // Applied Statistics, // Volume 46, Number 1, 1997, pages 146-156. // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 // before the first call. On each call, the routine increments N_DATA by 1, // and returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &A, &B, the shape parameters. // // Output, double &LAMBDA, the noncentrality parameter. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 25 static double a_vec[N_MAX] = { 5.0, 5.0, 5.0, 10.0, 10.0, 10.0, 20.0, 20.0, 20.0, 10.0, 10.0, 15.0, 20.0, 20.0, 20.0, 30.0, 30.0, 10.0, 10.0, 10.0, 15.0, 10.0, 12.0, 30.0, 35.0 }; static double b_vec[N_MAX] = { 5.0, 5.0, 5.0, 10.0, 10.0, 10.0, 20.0, 20.0, 20.0, 20.0, 10.0, 5.0, 10.0, 30.0, 50.0, 20.0, 40.0, 5.0, 10.0, 30.0, 20.0, 5.0, 17.0, 30.0, 30.0 }; static double fx_vec[N_MAX] = { 0.4563021, 0.1041337, 0.6022353, 0.9187770, 0.6008106, 0.0902850, 0.9998655, 0.9925997, 0.9641112, 0.9376626573, 0.7306817858, 0.1604256918, 0.1867485313, 0.6559386874, 0.9796881486, 0.1162386423, 0.9930430054, 0.0506899273, 0.1030959706, 0.9978417832, 0.2555552369, 0.0668307064, 0.0113601067, 0.7813366615, 0.8867126477 }; static double lambda_vec[N_MAX] = { 54.0, 140.0, 170.0, 54.0, 140.0, 250.0, 54.0, 140.0, 250.0, 150.0, 120.0, 80.0, 110.0, 65.0, 130.0, 80.0, 130.0, 20.0, 54.0, 80.0, 120.0, 55.0, 64.0, 140.0, 20.0 }; static double x_vec[N_MAX] = { 0.8640, 0.9000, 0.9560, 0.8686, 0.9000, 0.9000, 0.8787, 0.9000, 0.9220, 0.868, 0.900, 0.880, 0.850, 0.660, 0.720, 0.720, 0.800, 0.644, 0.700, 0.780, 0.760, 0.795, 0.560, 0.800, 0.670 }; 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; lambda = 0.0; x = 0.0; fx = 0.0; } else { a = a_vec[n_data-1]; b = b_vec[n_data-1]; lambda = lambda_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void beta_values ( int &n_data, double &x, double &y, double &fxy ) //****************************************************************************80 // // Purpose: // // BETA_VALUES returns some values of the Beta function. // // Discussion: // // Beta(X,Y) = ( Gamma(X) * Gamma(Y) ) / Gamma(X+Y) // // Both X and Y must be greater than 0. // // In Mathematica, the function can be evaluated by: // // Beta[X,Y] // // Properties: // // Beta(X,Y) = Beta(Y,X). // Beta(X,Y) = Integral ( 0 <= T <= 1 ) T^(X-1) (1-T)^(Y-1) dT. // Beta(X,Y) = Gamma(X) * Gamma(Y) / Gamma(X+Y) // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 13 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, double &X, &Y, the arguments of the function. // // Output, double &FXY, the value of the function. // { # define N_MAX 17 static double b_vec[N_MAX] = { 0.5000000000000000E+01, 0.2500000000000000E+01, 0.1666666666666667E+01, 0.1250000000000000E+01, 0.5000000000000000E+01, 0.2500000000000000E+01, 0.1000000000000000E+01, 0.1666666666666667E+00, 0.3333333333333333E-01, 0.7142857142857143E-02, 0.1587301587301587E-02, 0.2380952380952381E-01, 0.5952380952380952E-02, 0.1984126984126984E-02, 0.7936507936507937E-03, 0.3607503607503608E-03, 0.8325008325008325E-04 }; static double x_vec[N_MAX] = { 0.2E+00, 0.4E+00, 0.6E+00, 0.8E+00, 1.0E+00, 1.0E+00, 1.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, 6.0E+00, 6.0E+00, 6.0E+00, 6.0E+00, 7.0E+00 }; static double y_vec[N_MAX] = { 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 0.2E+00, 0.4E+00, 1.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, 7.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; x = 0.0; y = 0.0; fxy = 0.0; } else { x = x_vec[n_data-1]; y = y_vec[n_data-1]; fxy = b_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void binomial_values ( int &n_data, int &a, int &b, int &fx ) //****************************************************************************80 // // Purpose: // // BINOMIAL_VALUES returns some values of the binomial coefficients. // // Discussion: // // The formula for the binomial coefficient is // // C(N,K) = N! / ( K! * (N-K)! ) // // In Mathematica, the function can be evaluated by: // // Binomial[n,k] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 18 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 &A, &B, the arguments of the function. // // Output, int &FX, the value of the function. // { # define N_MAX 20 static int a_vec[N_MAX] = { 1, 6, 6, 6, 15, 15, 15, 15, 15, 15, 15, 25, 25, 25, 25, 25, 25, 25, 25, 25 }; static int b_vec[N_MAX] = { 0, 1, 3, 5, 1, 3, 5, 7, 9, 11, 13, 1, 3, 5, 7, 9, 11, 13, 15, 17 }; static int fx_vec[N_MAX] = { 1, 6, 20, 6, 15, 455, 3003, 6435, 5005, 1365, 105, 25, 2300, 53130, 480700, 2042975, 4457400, 5200300, 3268760, 1081575 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; a = 0; b = 0; fx = 0; } else { a = a_vec[n_data-1]; b = b_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void binomial_cdf_values ( int &n_data, int &a, double &b, int &x, double &fx ) //****************************************************************************80 // // Purpose: // // BINOMIAL_CDF_VALUES returns some values of the binomial CDF. // // Discussion: // // CDF(X)(A,B) is the probability of at most X successes in A trials, // given that the probability of success on a single trial is B. // // In Mathematica, the function can be evaluated by: // // Needs["Statistics`DiscreteDistributions`] // dist = BinomialDistribution [ n, p ] // CDF [ dist, x ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 15 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. // // 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, int &A, a parameter of the function. // // Output, double &B, a parameter of the function. // // Output, int &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 17 static int a_vec[N_MAX] = { 2, 2, 2, 2, 2, 4, 4, 4, 4, 10, 10, 10, 10, 10, 10, 10, 10 }; static double b_vec[N_MAX] = { 0.05E+00, 0.05E+00, 0.05E+00, 0.50E+00, 0.50E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.25E+00, 0.05E+00, 0.10E+00, 0.15E+00, 0.20E+00, 0.25E+00, 0.30E+00, 0.40E+00, 0.50E+00 }; static double fx_vec[N_MAX] = { 0.9025000000000000E+00, 0.9975000000000000E+00, 0.1000000000000000E+01, 0.2500000000000000E+00, 0.7500000000000000E+00, 0.3164062500000000E+00, 0.7382812500000000E+00, 0.9492187500000000E+00, 0.9960937500000000E+00, 0.9999363101685547E+00, 0.9983650626000000E+00, 0.9901259090013672E+00, 0.9672065024000000E+00, 0.9218730926513672E+00, 0.8497316674000000E+00, 0.6331032576000000E+00, 0.3769531250000000E+00 }; static int x_vec[N_MAX] = { 0, 1, 2, 0, 1, 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; a = 0; b = 0.0; x = 0; fx = 0.0; } else { a = a_vec[n_data-1]; b = b_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void bivariate_normal_cdf_values ( int &n_data, double &x, double &y, double &r, double &fxy ) //****************************************************************************80 // // Purpose: // // BIVARIATE_NORMAL_CDF_VALUES returns some values of the bivariate normal CDF. // // Discussion: // // FXY is the probability that two variables A and B, which are // related by a bivariate normal distribution with correlation R, // respectively satisfy A <= X and B <= Y. // // Mathematica can evaluate the bivariate normal CDF via the commands: // // < .0 = 0.0; // 1 = 1 => .1 = 0.5 // 2 = 10 => .01 = 0.25 // 3 = 11 => .11 = 0.75 // 4 = 100 => .001 = 0.125 // 5 = 101 => .101 = 0.625 // 6 = 110 => .011 = 0.375 // 7 = 111 => .111 = 0.875 // 8 = 1000 => .0001 = 0.0625 // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 28 June 2004 // // Author: // // John Burkardt // // Reference: // // John Halton, // On the efficiency of certain quasi-random sequences of points // in evaluating multi-dimensional integrals, // Numerische Mathematik, // Volume 2, pages 84-90, 1960. // // John Hammersley, // Monte Carlo methods for solving multivariable problems, // Proceedings of the New York Academy of Science, // Volume 86, pages 844-874, 1960. // // J G van der Corput, // Verteilungsfunktionen, // Proc Akad Amsterdam, // Volume 38, 1935, // Volume 39, 1936. // // 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 &BASE, the base of the sequence. // // Output, int &SEED, the index of the element of the sequence. // // Output, double &VALUE, the value of the SEED-th element of the // van der Corput sequence in base BASE. // { # define N_MAX 75 static int base_vec[N_MAX] = { 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 3, 4, 5, 7, 11, 13, 2, 3, 4, 5, 7, 11, 13, 2, 3, 4, 5, 7, 11, 13, 2, 3, 4, 5, 7, 11, 13, 29, 29, 29, 29, 29, 71, 71, 71, 71, 71, 173, 173, 173, 173, 173, 409, 409, 409, 409, 409 }; static int seed_vec[N_MAX] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 10, 10, 10, 10, 10, 10, 100, 100, 100, 100, 100, 100, 100, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 1000, 1001, 1002, 1003, 1004, 1000, 1001, 1002, 1003, 1004, 1000, 1001, 1002, 1003, 1004, 1000, 1001, 1002, 1003, 1004 }; static double value_vec[N_MAX] = { 0.0000000000000000E+00, 0.5000000000000000E+00, 0.2500000000000000E+00, 0.7500000000000000E+00, 0.1250000000000000E+00, 0.6250000000000000E+00, 0.3750000000000000E+00, 0.8750000000000000E+00, 0.0625000000000000E+00, 0.0000000000000000E+00, 0.3333333333333333E+00, 0.6666666666666666E+00, 0.1111111111111111E+00, 0.4444444444444444E+00, 0.7777777777777777E+00, 0.2222222222222222E+00, 0.5555555555555556E+00, 0.8888888888888888E+00, 0.0000000000000000E+00, 0.2500000000000000E+00, 0.5000000000000000E+00, 0.7500000000000000E+00, 0.0625000000000000E+00, 0.3125000000000000E+00, 0.5625000000000000E+00, 0.8125000000000000E+00, 0.1250000000000000E+00, 0.3125000000000000E+00, 0.3703703703703703E+00, 0.6250000000000000E+00, 0.0800000000000000E+00, 0.4489795918367347E+00, 0.9090909090909092E+00, 0.7692307692307693E+00, 0.1484375000000000E+00, 0.4115226337448559E+00, 0.0976562500000000E+00, 0.0320000000000000E+00, 0.2915451895043731E+00, 0.1652892561983471E+00, 0.7337278106508875E+00, 0.0927734375000000E+00, 0.3475080018289895E+00, 0.1708984375000000E+00, 0.0051200000000000E+00, 0.9162848812994586E+00, 0.9316303531179565E+00, 0.9904415111515704E+00, 0.0347290039062500E+00, 0.3861200020322105E+00, 0.0189208984375000E+00, 0.0005120000000000E+00, 0.5749985125245433E+00, 0.1529950140017758E+00, 0.2459297643639929E+00, 0.4887449259912255E+00, 0.5232276846119153E+00, 0.5577104432326049E+00, 0.5921932018532945E+00, 0.6266759604739842E+00, 0.0872842689942472E+00, 0.1013687760365007E+00, 0.1154532830787542E+00, 0.1295377901210077E+00, 0.1436222971632613E+00, 0.7805138828560928E+00, 0.7862942296769020E+00, 0.7920745764977113E+00, 0.7978549233185205E+00, 0.8036352701393298E+00, 0.4449997309915651E+00, 0.4474447187666262E+00, 0.4498897065416874E+00, 0.4523346943167484E+00, 0.4547796820918096E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; base = 0; seed = 0; value = 0.0; } else { base = base_vec[n_data-1]; seed = seed_vec[n_data-1]; value = value_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void viscosity_values ( int &n_data, double &tc, double &p, double &eta ) //****************************************************************************80 // // Purpose: // // VISCOSITY_VALUES returns some values of the viscosity function. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 04 February 2002 // // Author: // // John Burkardt // // Reference: // // Lester Haar, John Gallagher and George Kell, // NBS/NRC Steam Tables: // Thermodynamic and Transport Properties and Computer Programs // for Vapor and Liquid States of Water in SI Units, // Hemisphere Publishing Corporation, Washington, 1984, // TJ270.H3, page 263. // // 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 &TC, the temperature, in degrees Celsius. // // Output, double &P, the pressure, in bar. // // Output, double &ETA, the viscosity, in MegaPascal seconds. // { # define N_MAX 34 static double eta_vec[N_MAX] = { 1792.0E+00, 1791.0E+00, 1790.0E+00, 1786.0E+00, 1780.0E+00, 1775.0E+00, 1769.0E+00, 1764.0E+00, 1759.0E+00, 1754.0E+00, 1749.0E+00, 1744.0E+00, 1739.0E+00, 1735.0E+00, 1731.0E+00, 1722.0E+00, 1714.0E+00, 1707.0E+00, 1700.0E+00, 1694.0E+00, 1687.0E+00, 1682.0E+00, 1676.0E+00, 1667.0E+00, 1659.0E+00, 1653.0E+00, 890.8E+00, 547.1E+00, 378.4E+00, 12.28E+00, 16.18E+00, 24.45E+00, 32.61E+00, 40.38E+00 }; static double p_vec[N_MAX] = { 1.0E+00, 5.0E+00, 10.0E+00, 25.0E+00, 50.0E+00, 75.0E+00, 100.0E+00, 125.0E+00, 150.0E+00, 175.0E+00, 200.0E+00, 225.0E+00, 250.0E+00, 275.0E+00, 300.0E+00, 350.0E+00, 400.0E+00, 450.0E+00, 500.0E+00, 550.0E+00, 600.0E+00, 650.0E+00, 700.0E+00, 800.0E+00, 900.0E+00, 1000.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00, 1.0E+00 }; static double tc_vec[N_MAX] = { 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 25.0E+00, 50.0E+00, 75.0E+00, 100.0E+00, 200.0E+00, 400.0E+00, 600.0E+00, 800.0E+00 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; tc = 0.0; p = 0.0; eta = 0.0; } else { tc = tc_vec[n_data-1]; p = p_vec[n_data-1]; eta = eta_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void von_mises_cdf_values ( int &n_data, double &a, double &b, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // VON_MISES_CDF_VALUES returns some values of the von Mises CDF. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 08 December 2004 // // Author: // // John Burkardt // // Reference: // // Kanti Mardia, Peter Jupp, // Directional Statistics, // Wiley, 2000, QA276.M335 // // Parameters: // // Input/output, int &N_DATA. The user sets N_DATA to 0 before the // first call. On each call, the routine increments N_DATA by 1, and // returns the corresponding data; when there is no more data, the // output value of N_DATA will be 0 again. // // Output, double &A, &B, the parameters of the function. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 23 static double a_vec[N_MAX] = { 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, -0.2E+01, -0.1E+01, 0.0E+01, 0.1E+01, 0.2E+01, 0.3E+01, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00, 0.0E+00 }; static double b_vec[N_MAX] = { 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.1E+01, 0.2E+01, 0.2E+01, 0.2E+01, 0.2E+01, 0.2E+01, 0.2E+01, 0.3E+01, 0.3E+01, 0.3E+01, 0.3E+01, 0.3E+01, 0.3E+01, 0.0E+00, 0.1E+01, 0.2E+01, 0.3E+01, 0.4E+01, 0.5E+01 }; static double fx_vec[N_MAX] = { 0.2535089956281180E-01, 0.1097539041177346E+00, 0.5000000000000000E+00, 0.8043381312498558E+00, 0.9417460124555197E+00, 0.5000000000000000E+00, 0.6018204118446155E+00, 0.6959356933122230E+00, 0.7765935901304593E+00, 0.8410725934916615E+00, 0.8895777369550366E+00, 0.9960322705517925E+00, 0.9404336090170247E+00, 0.5000000000000000E+00, 0.5956639098297530E-01, 0.3967729448207649E-02, 0.2321953958111930E-03, 0.6250000000000000E+00, 0.7438406999109122E+00, 0.8369224904294019E+00, 0.8941711407897124E+00, 0.9291058600568743E+00, 0.9514289900655436E+00 }; static double x_vec[N_MAX] = { -0.2617993977991494E+01, -0.1570796326794897E+01, 0.0000000000000000E+00, 0.1047197551196598E+01, 0.2094395102393195E+01, 0.1000000000000000E+01, 0.1200000000000000E+01, 0.1400000000000000E+01, 0.1600000000000000E+01, 0.1800000000000000E+01, 0.2000000000000000E+01, 0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 0.7853981633974483E+00, 0.7853981633974483E+00, 0.7853981633974483E+00, 0.7853981633974483E+00, 0.7853981633974483E+00, 0.7853981633974483E+00 }; 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; x = 0.0; fx = 0.0; } else { a = a_vec[n_data-1]; b = b_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void weekday_values ( int &n_data, int &y, int &m, int &d, int &w ) //****************************************************************************80 // // Purpose: // // WEEKDAY_VALUES returns the day of the week for various dates. // // Discussion: // // The CE or Common Era calendar is used, under the // hybrid Julian/Gregorian Calendar, with a transition from Julian // to Gregorian. The day after 04 October 1582 was 15 October 1582. // // The year before 1 AD or CE is 1 BC or BCE. In this data set, // years BC/BCE are indicated by a negative year value. // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 26 May 2012 // // Author: // // John Burkardt // // Reference: // // Edward Reingold, Nachum Dershowitz, // Calendrical Calculations: The Millennium Edition, // Cambridge University Press, 2001, // ISBN: 0 521 77752 6 // LC: CE12.R45. // // 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 &Y, &M, &D, the Common Era date. // // Output, int &W, the day of the week. Sunday = 1. // { # define N_MAX 34 static int d_vec[N_MAX] = { 30, 8, 26, 3, 7, 18, 7, 19, 14, 18, 16, 3, 26, 20, 4, 25, 31, 9, 24, 10, 30, 24, 19, 2, 27, 19, 25, 29, 19, 7, 17, 25, 10, 18 }; static int m_vec[N_MAX] = { 7, 12, 9, 10, 1, 5, 11, 4, 10, 5, 3, 3, 3, 4, 6, 1, 3, 9, 2, 6, 6, 7, 6, 8, 3, 4, 8, 9, 4, 10, 3, 2, 11, 7 }; static int w_vec[N_MAX] = { 1, 4, 4, 1, 4, 2, 7, 1, 7, 1, 6, 7, 6, 1, 1, 4, 7, 7, 7, 4, 1, 6, 1, 2, 4, 1, 1, 2, 2, 5, 3, 1, 4, 1 }; static int y_vec[N_MAX] = { - 587, - 169, 70, 135, 470, 576, 694, 1013, 1066, 1096, 1190, 1240, 1288, 1298, 1391, 1436, 1492, 1553, 1560, 1648, 1680, 1716, 1768, 1819, 1839, 1903, 1929, 1941, 1943, 1943, 1992, 1996, 2038, 2094 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; y = 0; m = 0; d = 0; w = 0; } else { y = y_vec[n_data-1]; m = m_vec[n_data-1]; d = d_vec[n_data-1]; w = w_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void weibull_cdf_values ( int &n_data, double &alpha, double &beta, double &x, double &fx ) //****************************************************************************80 // // Purpose: // // WEIBULL_CDF_VALUES returns some values of the Weibull CDF. // // Discussion: // // In Mathematica, the function can be evaluated by: // // Needs["Statistics`ContinuousDistributions`"] // dist = WeibullDistribution [ alpha, beta ] // CDF [ dist, x ] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 31 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, double &ALPHA, the first parameter of the distribution. // // Output, double &BETA, the second parameter of the distribution. // // Output, double &X, the argument of the function. // // Output, double &FX, the value of the function. // { # define N_MAX 12 static double alpha_vec[N_MAX] = { 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.1000000000000000E+01, 0.2000000000000000E+01, 0.3000000000000000E+01, 0.4000000000000000E+01, 0.5000000000000000E+01 }; static double beta_vec[N_MAX] = { 0.5000000000000000E+00, 0.5000000000000000E+00, 0.5000000000000000E+00, 0.5000000000000000E+00, 0.2000000000000000E+01, 0.3000000000000000E+01, 0.4000000000000000E+01, 0.5000000000000000E+01, 0.2000000000000000E+01, 0.2000000000000000E+01, 0.2000000000000000E+01, 0.2000000000000000E+01 }; static double fx_vec[N_MAX] = { 0.8646647167633873E+00, 0.9816843611112658E+00, 0.9975212478233336E+00, 0.9996645373720975E+00, 0.6321205588285577E+00, 0.4865828809674080E+00, 0.3934693402873666E+00, 0.3296799539643607E+00, 0.8946007754381357E+00, 0.9657818816883340E+00, 0.9936702845725143E+00, 0.9994964109502630E+00 }; static double x_vec[N_MAX] = { 0.1000000000000000E+01, 0.2000000000000000E+01, 0.3000000000000000E+01, 0.4000000000000000E+01, 0.2000000000000000E+01, 0.2000000000000000E+01, 0.2000000000000000E+01, 0.2000000000000000E+01, 0.3000000000000000E+01, 0.3000000000000000E+01, 0.3000000000000000E+01, 0.3000000000000000E+01 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; alpha = 0.0; beta = 0.0; x = 0.0; fx = 0.0; } else { alpha = alpha_vec[n_data-1]; beta = beta_vec[n_data-1]; x = x_vec[n_data-1]; fx = fx_vec[n_data-1]; } return; # undef N_MAX } //****************************************************************************80 void zeta_values ( int &n_data, int &n, double &zeta ) //****************************************************************************80 // // Purpose: // // ZETA_VALUES returns some values of the Riemann Zeta function. // // Discussion: // // ZETA(N) = sum ( 1 <= I < +oo ) 1 / I^N // // In Mathematica, the function can be evaluated by: // // Zeta[n] // // Licensing: // // This code is distributed under the GNU LGPL license. // // Modified: // // 12 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 argument of the Zeta function. // // Output, double &ZETA, the value of the Zeta function. // { # define N_MAX 15 static int n_vec[N_MAX] = { 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 16, 20, 30, 40 }; static double zeta_vec[N_MAX] = { 0.164493406684822643647E+01, 0.120205690315959428540E+01, 0.108232323371113819152E+01, 0.103692775514336992633E+01, 0.101734306198444913971E+01, 0.100834927738192282684E+01, 0.100407735619794433939E+01, 0.100200839292608221442E+01, 0.100099457512781808534E+01, 0.100049418860411946456E+01, 0.100024608655330804830E+01, 0.100001528225940865187E+01, 0.100000095396203387280E+01, 0.100000000093132743242E+01, 0.100000000000090949478E+01 }; if ( n_data < 0 ) { n_data = 0; } n_data = n_data + 1; if ( N_MAX < n_data ) { n_data = 0; n = 0; zeta = 0.0; } else { n = n_vec[n_data-1]; zeta = zeta_vec[n_data-1]; } return; # undef N_MAX }