diff options
Diffstat (limited to 'gsl-1.9/specfunc/expint.c')
| -rw-r--r-- | gsl-1.9/specfunc/expint.c | 518 |
1 files changed, 518 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/expint.c b/gsl-1.9/specfunc/expint.c new file mode 100644 index 0000000..45957ea --- /dev/null +++ b/gsl-1.9/specfunc/expint.c @@ -0,0 +1,518 @@ +/* specfunc/expint.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 Gerard Jungman + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +/* Author: G. Jungman */ + +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sf_expint.h> + +#include "error.h" + +#include "chebyshev.h" +#include "cheb_eval.c" + +/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ + +/* + Chebyshev expansions: based on SLATEC e1.f, W. Fullerton + + Series for AE11 on the interval -1.00000D-01 to 0. + with weighted error 1.76E-17 + log weighted error 16.75 + significant figures required 15.70 + decimal places required 17.55 + + + Series for AE12 on the interval -2.50000D-01 to -1.00000D-01 + with weighted error 5.83E-17 + log weighted error 16.23 + significant figures required 15.76 + decimal places required 16.93 + + + Series for E11 on the interval -4.00000D+00 to -1.00000D+00 + with weighted error 1.08E-18 + log weighted error 17.97 + significant figures required 19.02 + decimal places required 18.61 + + + Series for E12 on the interval -1.00000D+00 to 1.00000D+00 + with weighted error 3.15E-18 + log weighted error 17.50 + approx significant figures required 15.8 + decimal places required 18.10 + + + Series for AE13 on the interval 2.50000D-01 to 1.00000D+00 + with weighted error 2.34E-17 + log weighted error 16.63 + significant figures required 16.14 + decimal places required 17.33 + + + Series for AE14 on the interval 0. to 2.50000D-01 + with weighted error 5.41E-17 + log weighted error 16.27 + significant figures required 15.38 + decimal places required 16.97 +*/ + +static double AE11_data[39] = { + 0.121503239716065790, + -0.065088778513550150, + 0.004897651357459670, + -0.000649237843027216, + 0.000093840434587471, + 0.000000420236380882, + -0.000008113374735904, + 0.000002804247688663, + 0.000000056487164441, + -0.000000344809174450, + 0.000000058209273578, + 0.000000038711426349, + -0.000000012453235014, + -0.000000005118504888, + 0.000000002148771527, + 0.000000000868459898, + -0.000000000343650105, + -0.000000000179796603, + 0.000000000047442060, + 0.000000000040423282, + -0.000000000003543928, + -0.000000000008853444, + -0.000000000000960151, + 0.000000000001692921, + 0.000000000000607990, + -0.000000000000224338, + -0.000000000000200327, + -0.000000000000006246, + 0.000000000000045571, + 0.000000000000016383, + -0.000000000000005561, + -0.000000000000006074, + -0.000000000000000862, + 0.000000000000001223, + 0.000000000000000716, + -0.000000000000000024, + -0.000000000000000201, + -0.000000000000000082, + 0.000000000000000017 +}; +static cheb_series AE11_cs = { + AE11_data, + 38, + -1, 1, + 20 +}; + +static double AE12_data[25] = { + 0.582417495134726740, + -0.158348850905782750, + -0.006764275590323141, + 0.005125843950185725, + 0.000435232492169391, + -0.000143613366305483, + -0.000041801320556301, + -0.000002713395758640, + 0.000001151381913647, + 0.000000420650022012, + 0.000000066581901391, + 0.000000000662143777, + -0.000000002844104870, + -0.000000000940724197, + -0.000000000177476602, + -0.000000000015830222, + 0.000000000002905732, + 0.000000000001769356, + 0.000000000000492735, + 0.000000000000093709, + 0.000000000000010707, + -0.000000000000000537, + -0.000000000000000716, + -0.000000000000000244, + -0.000000000000000058 +}; +static cheb_series AE12_cs = { + AE12_data, + 24, + -1, 1, + 15 +}; + +static double E11_data[19] = { + -16.11346165557149402600, + 7.79407277874268027690, + -1.95540581886314195070, + 0.37337293866277945612, + -0.05692503191092901938, + 0.00721107776966009185, + -0.00078104901449841593, + 0.00007388093356262168, + -0.00000620286187580820, + 0.00000046816002303176, + -0.00000003209288853329, + 0.00000000201519974874, + -0.00000000011673686816, + 0.00000000000627627066, + -0.00000000000031481541, + 0.00000000000001479904, + -0.00000000000000065457, + 0.00000000000000002733, + -0.00000000000000000108 +}; +static cheb_series E11_cs = { + E11_data, + 18, + -1, 1, + 13 +}; + +static double E12_data[16] = { + -0.03739021479220279500, + 0.04272398606220957700, + -0.13031820798497005440, + 0.01441912402469889073, + -0.00134617078051068022, + 0.00010731029253063780, + -0.00000742999951611943, + 0.00000045377325690753, + -0.00000002476417211390, + 0.00000000122076581374, + -0.00000000005485141480, + 0.00000000000226362142, + -0.00000000000008635897, + 0.00000000000000306291, + -0.00000000000000010148, + 0.00000000000000000315 +}; +static cheb_series E12_cs = { + E12_data, + 15, + -1, 1, + 10 +}; + +static double AE13_data[25] = { + -0.605773246640603460, + -0.112535243483660900, + 0.013432266247902779, + -0.001926845187381145, + 0.000309118337720603, + -0.000053564132129618, + 0.000009827812880247, + -0.000001885368984916, + 0.000000374943193568, + -0.000000076823455870, + 0.000000016143270567, + -0.000000003466802211, + 0.000000000758754209, + -0.000000000168864333, + 0.000000000038145706, + -0.000000000008733026, + 0.000000000002023672, + -0.000000000000474132, + 0.000000000000112211, + -0.000000000000026804, + 0.000000000000006457, + -0.000000000000001568, + 0.000000000000000383, + -0.000000000000000094, + 0.000000000000000023 +}; +static cheb_series AE13_cs = { + AE13_data, + 24, + -1, 1, + 15 +}; + +static double AE14_data[26] = { + -0.18929180007530170, + -0.08648117855259871, + 0.00722410154374659, + -0.00080975594575573, + 0.00010999134432661, + -0.00001717332998937, + 0.00000298562751447, + -0.00000056596491457, + 0.00000011526808397, + -0.00000002495030440, + 0.00000000569232420, + -0.00000000135995766, + 0.00000000033846628, + -0.00000000008737853, + 0.00000000002331588, + -0.00000000000641148, + 0.00000000000181224, + -0.00000000000052538, + 0.00000000000015592, + -0.00000000000004729, + 0.00000000000001463, + -0.00000000000000461, + 0.00000000000000148, + -0.00000000000000048, + 0.00000000000000016, + -0.00000000000000005 +}; +static cheb_series AE14_cs = { + AE14_data, + 25, + -1, 1, + 13 +}; + + + +/* implementation for E1, allowing for scaling by exp(x) */ +static +int expint_E1_impl(const double x, gsl_sf_result * result, const int scale) +{ + const double xmaxt = -GSL_LOG_DBL_MIN; /* XMAXT = -LOG (R1MACH(1)) */ + const double xmax = xmaxt - log(xmaxt); /* XMAX = XMAXT - LOG(XMAXT) */ + + /* CHECK_POINTER(result) */ + + if(x < -xmax && !scale) { + OVERFLOW_ERROR(result); + } + else if(x <= -10.0) { + const double s = 1.0/x * ( scale ? 1.0 : exp(-x) ); + gsl_sf_result result_c; + cheb_eval_e(&AE11_cs, 20.0/x+1.0, &result_c); + result->val = s * (1.0 + result_c.val); + result->err = s * result_c.err; + result->err += 2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) * fabs(result->val); + return GSL_SUCCESS; + } + else if(x <= -4.0) { + const double s = 1.0/x * ( scale ? 1.0 : exp(-x) ); + gsl_sf_result result_c; + cheb_eval_e(&AE12_cs, (40.0/x+7.0)/3.0, &result_c); + result->val = s * (1.0 + result_c.val); + result->err = s * result_c.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x <= -1.0) { + const double ln_term = -log(fabs(x)); + const double scale_factor = ( scale ? exp(x) : 1.0 ); + gsl_sf_result result_c; + cheb_eval_e(&E11_cs, (2.0*x+5.0)/3.0, &result_c); + result->val = scale_factor * (ln_term + result_c.val); + result->err = scale_factor * (result_c.err + GSL_DBL_EPSILON * fabs(ln_term)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x == 0.0) { + DOMAIN_ERROR(result); + } + else if(x <= 1.0) { + const double ln_term = -log(fabs(x)); + const double scale_factor = ( scale ? exp(x) : 1.0 ); + gsl_sf_result result_c; + cheb_eval_e(&E12_cs, x, &result_c); + result->val = scale_factor * (ln_term - 0.6875 + x + result_c.val); + result->err = scale_factor * (result_c.err + GSL_DBL_EPSILON * fabs(ln_term)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x <= 4.0) { + const double s = 1.0/x * ( scale ? 1.0 : exp(-x) ); + gsl_sf_result result_c; + cheb_eval_e(&AE13_cs, (8.0/x-5.0)/3.0, &result_c); + result->val = s * (1.0 + result_c.val); + result->err = s * result_c.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(x <= xmax || scale) { + const double s = 1.0/x * ( scale ? 1.0 : exp(-x) ); + gsl_sf_result result_c; + cheb_eval_e(&AE14_cs, 8.0/x-1.0, &result_c); + result->val = s * (1.0 + result_c.val); + result->err = s * (GSL_DBL_EPSILON + result_c.err); + result->err += 2.0 * (x + 1.0) * GSL_DBL_EPSILON * fabs(result->val); + if(result->val == 0.0) + UNDERFLOW_ERROR(result); + else + return GSL_SUCCESS; + } + else { + UNDERFLOW_ERROR(result); + } +} + + +static +int expint_E2_impl(const double x, gsl_sf_result * result, const int scale) +{ + const double xmaxt = -GSL_LOG_DBL_MIN; + const double xmax = xmaxt - log(xmaxt); + + /* CHECK_POINTER(result) */ + + if(x < -xmax && !scale) { + OVERFLOW_ERROR(result); + } + else if (x == 0.0) { + result->val = (scale ? 1.0 : 1.0); + result->err = 0.0; + return GSL_SUCCESS; + } else if(x < 100.0) { + const double ex = ( scale ? 1.0 : exp(-x) ); + gsl_sf_result result_E1; + int stat_E1 = expint_E1_impl(x, &result_E1, scale); + result->val = ex - x*result_E1.val; + result->err = GSL_DBL_EPSILON*ex + fabs(x) * result_E1.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_E1; + } + else if(x < xmax || scale) { + const double s = ( scale ? 1.0 : exp(-x) ); + const double c1 = -2.0; + const double c2 = 6.0; + const double c3 = -24.0; + const double c4 = 120.0; + const double c5 = -720.0; + const double c6 = 5040.0; + const double c7 = -40320.0; + const double c8 = 362880.0; + const double c9 = -3628800.0; + const double c10 = 39916800.0; + const double c11 = -479001600.0; + const double c12 = 6227020800.0; + const double c13 = -87178291200.0; + const double y = 1.0/x; + const double sum6 = c6+y*(c7+y*(c8+y*(c9+y*(c10+y*(c11+y*(c12+y*c13)))))); + const double sum = y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*sum6))))); + result->val = s * (1.0 + sum)/x; + result->err = 2.0 * (x + 1.0) * GSL_DBL_EPSILON * result->val; + if(result->val == 0.0) + UNDERFLOW_ERROR(result); + else + return GSL_SUCCESS; + } + else { + UNDERFLOW_ERROR(result); + } +} + + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + + +int gsl_sf_expint_E1_e(const double x, gsl_sf_result * result) +{ + return expint_E1_impl(x, result, 0); +} + + +int gsl_sf_expint_E1_scaled_e(const double x, gsl_sf_result * result) +{ + return expint_E1_impl(x, result, 1); +} + + +int gsl_sf_expint_E2_e(const double x, gsl_sf_result * result) +{ + return expint_E2_impl(x, result, 0); +} + + +int gsl_sf_expint_E2_scaled_e(const double x, gsl_sf_result * result) +{ + return expint_E2_impl(x, result, 1); +} + + +int gsl_sf_expint_Ei_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + { + int status = gsl_sf_expint_E1_e(-x, result); + result->val = -result->val; + return status; + } +} + + +int gsl_sf_expint_Ei_scaled_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + { + int status = gsl_sf_expint_E1_scaled_e(-x, result); + result->val = -result->val; + return status; + } +} + + +#if 0 +static double recurse_En(int n, double x, double E1) +{ + int i; + double En = E1; + double ex = exp(-x); + for(i=2; i<=n; i++) { + En = 1./(i-1) * (ex - x * En); + } + return En; +} +#endif + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_expint_E1(const double x) +{ + EVAL_RESULT(gsl_sf_expint_E1_e(x, &result)); +} + +double gsl_sf_expint_E1_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_expint_E1_scaled_e(x, &result)); +} + +double gsl_sf_expint_E2(const double x) +{ + EVAL_RESULT(gsl_sf_expint_E2_e(x, &result)); +} + +double gsl_sf_expint_E2_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_expint_E2_scaled_e(x, &result)); +} + +double gsl_sf_expint_Ei(const double x) +{ + EVAL_RESULT(gsl_sf_expint_Ei_e(x, &result)); +} + +double gsl_sf_expint_Ei_scaled(const double x) +{ + EVAL_RESULT(gsl_sf_expint_Ei_scaled_e(x, &result)); +} |