diff options
Diffstat (limited to 'gsl-1.9/specfunc/sinint.c')
| -rw-r--r-- | gsl-1.9/specfunc/sinint.c | 402 |
1 files changed, 402 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/sinint.c b/gsl-1.9/specfunc/sinint.c new file mode 100644 index 0000000..0d56f12 --- /dev/null +++ b/gsl-1.9/specfunc/sinint.c @@ -0,0 +1,402 @@ +/* specfunc/sinint.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 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_trig.h> +#include <gsl/gsl_sf_expint.h> + +#include "error.h" + +#include "chebyshev.h" +#include "cheb_eval.c" + +/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/ + +/* based on SLATEC r9sifg.f, W. Fullerton */ + +/* + series for f1 on the interval 2.00000e-02 to 6.25000e-02 + with weighted error 2.82e-17 + log weighted error 16.55 + significant figures required 15.36 + decimal places required 17.20 +*/ +static double f1_data[20] = { + -0.1191081969051363610, + -0.0247823144996236248, + 0.0011910281453357821, + -0.0000927027714388562, + 0.0000093373141568271, + -0.0000011058287820557, + 0.0000001464772071460, + -0.0000000210694496288, + 0.0000000032293492367, + -0.0000000005206529618, + 0.0000000000874878885, + -0.0000000000152176187, + 0.0000000000027257192, + -0.0000000000005007053, + 0.0000000000000940241, + -0.0000000000000180014, + 0.0000000000000035063, + -0.0000000000000006935, + 0.0000000000000001391, + -0.0000000000000000282 +}; +static cheb_series f1_cs = { + f1_data, + 19, + -1, 1, + 10 +}; + +/* + + series for f2 on the interval 0.00000e+00 to 2.00000e-02 + with weighted error 4.32e-17 + log weighted error 16.36 + significant figures required 14.75 + decimal places required 17.10 +*/ +static double f2_data[29] = { + -0.0348409253897013234, + -0.0166842205677959686, + 0.0006752901241237738, + -0.0000535066622544701, + 0.0000062693421779007, + -0.0000009526638801991, + 0.0000001745629224251, + -0.0000000368795403065, + 0.0000000087202677705, + -0.0000000022601970392, + 0.0000000006324624977, + -0.0000000001888911889, + 0.0000000000596774674, + -0.0000000000198044313, + 0.0000000000068641396, + -0.0000000000024731020, + 0.0000000000009226360, + -0.0000000000003552364, + 0.0000000000001407606, + -0.0000000000000572623, + 0.0000000000000238654, + -0.0000000000000101714, + 0.0000000000000044259, + -0.0000000000000019634, + 0.0000000000000008868, + -0.0000000000000004074, + 0.0000000000000001901, + -0.0000000000000000900, + 0.0000000000000000432 +}; +static cheb_series f2_cs = { + f2_data, + 28, + -1, 1, + 14 +}; + +/* + + series for g1 on the interval 2.00000e-02 to 6.25000e-02 + with weighted error 5.48e-17 + log weighted error 16.26 + significant figures required 15.47 + decimal places required 16.92 +*/ +static double g1_data[21] = { + -0.3040578798253495954, + -0.0566890984597120588, + 0.0039046158173275644, + -0.0003746075959202261, + 0.0000435431556559844, + -0.0000057417294453025, + 0.0000008282552104503, + -0.0000001278245892595, + 0.0000000207978352949, + -0.0000000035313205922, + 0.0000000006210824236, + -0.0000000001125215474, + 0.0000000000209088918, + -0.0000000000039715832, + 0.0000000000007690431, + -0.0000000000001514697, + 0.0000000000000302892, + -0.0000000000000061400, + 0.0000000000000012601, + -0.0000000000000002615, + 0.0000000000000000548 +}; +static cheb_series g1_cs = { + g1_data, + 20, + -1, 1, + 13 +}; + +/* + + series for g2 on the interval 0.00000e+00 to 2.00000e-02 + with weighted error 5.01e-17 + log weighted error 16.30 + significant figures required 15.12 + decimal places required 17.07 +*/ +static double g2_data[34] = { + -0.0967329367532432218, + -0.0452077907957459871, + 0.0028190005352706523, + -0.0002899167740759160, + 0.0000407444664601121, + -0.0000071056382192354, + 0.0000014534723163019, + -0.0000003364116512503, + 0.0000000859774367886, + -0.0000000238437656302, + 0.0000000070831906340, + -0.0000000022318068154, + 0.0000000007401087359, + -0.0000000002567171162, + 0.0000000000926707021, + -0.0000000000346693311, + 0.0000000000133950573, + -0.0000000000053290754, + 0.0000000000021775312, + -0.0000000000009118621, + 0.0000000000003905864, + -0.0000000000001708459, + 0.0000000000000762015, + -0.0000000000000346151, + 0.0000000000000159996, + -0.0000000000000075213, + 0.0000000000000035970, + -0.0000000000000017530, + 0.0000000000000008738, + -0.0000000000000004487, + 0.0000000000000002397, + -0.0000000000000001347, + 0.0000000000000000801, + -0.0000000000000000501 +}; +static cheb_series g2_cs = { + g2_data, + 33, + -1, 1, + 20 +}; + + +/* x >= 4.0 */ +static void fg_asymp(const double x, gsl_sf_result * f, gsl_sf_result * g) +{ + /* + xbig = sqrt (1.0/r1mach(3)) + xmaxf = exp (amin1(-alog(r1mach(1)), alog(r1mach(2))) - 0.01) + xmaxg = 1.0/sqrt(r1mach(1)) + xbnd = sqrt(50.0) + */ + const double xbig = 1.0/GSL_SQRT_DBL_EPSILON; + const double xmaxf = 1.0/GSL_DBL_MIN; + const double xmaxg = 1.0/GSL_SQRT_DBL_MIN; + const double xbnd = 7.07106781187; + + const double x2 = x*x; + + if(x <= xbnd) { + gsl_sf_result result_c1; + gsl_sf_result result_c2; + cheb_eval_e(&f1_cs, (1.0/x2-0.04125)/0.02125, &result_c1); + cheb_eval_e(&g1_cs, (1.0/x2-0.04125)/0.02125, &result_c2); + f->val = (1.0 + result_c1.val)/x; + g->val = (1.0 + result_c2.val)/x2; + f->err = result_c1.err/x + 2.0 * GSL_DBL_EPSILON * fabs(f->val); + g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val); + } + else if(x <= xbig) { + gsl_sf_result result_c1; + gsl_sf_result result_c2; + cheb_eval_e(&f2_cs, 100.0/x2-1.0, &result_c1); + cheb_eval_e(&g2_cs, 100.0/x2-1.0, &result_c2); + f->val = (1.0 + result_c1.val)/x; + g->val = (1.0 + result_c2.val)/x2; + f->err = result_c1.err/x + 2.0 * GSL_DBL_EPSILON * fabs(f->val); + g->err = result_c2.err/x2 + 2.0 * GSL_DBL_EPSILON * fabs(g->val); + } + else { + f->val = (x < xmaxf ? 1.0/x : 0.0); + g->val = (x < xmaxg ? 1.0/x2 : 0.0); + f->err = 2.0 * GSL_DBL_EPSILON * fabs(f->val); + g->err = 2.0 * GSL_DBL_EPSILON * fabs(g->val); + } + + return; +} + + +/* based on SLATEC si.f, W. Fullerton + + series for si on the interval 0.00000e+00 to 1.60000e+01 + with weighted error 1.22e-17 + log weighted error 16.91 + significant figures required 16.37 + decimal places required 17.45 +*/ + +static double si_data[12] = { + -0.1315646598184841929, + -0.2776578526973601892, + 0.0354414054866659180, + -0.0025631631447933978, + 0.0001162365390497009, + -0.0000035904327241606, + 0.0000000802342123706, + -0.0000000013562997693, + 0.0000000000179440722, + -0.0000000000001908387, + 0.0000000000000016670, + -0.0000000000000000122 +}; + +static cheb_series si_cs = { + si_data, + 11, + -1, 1, + 9 +}; + +/* + series for ci on the interval 0.00000e+00 to 1.60000e+01 + with weighted error 1.94e-18 + log weighted error 17.71 + significant figures required 17.74 + decimal places required 18.27 +*/ +static double ci_data[13] = { + -0.34004281856055363156, + -1.03302166401177456807, + 0.19388222659917082877, + -0.01918260436019865894, + 0.00110789252584784967, + -0.00004157234558247209, + 0.00000109278524300229, + -0.00000002123285954183, + 0.00000000031733482164, + -0.00000000000376141548, + 0.00000000000003622653, + -0.00000000000000028912, + 0.00000000000000000194 +}; +static cheb_series ci_cs = { + ci_data, + 12, + -1, 1, + 9 +}; + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int gsl_sf_Si_e(const double x, gsl_sf_result * result) +{ + double ax = fabs(x); + + /* CHECK_POINTER(result) */ + + if(ax < GSL_SQRT_DBL_EPSILON) { + result->val = x; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(ax <= 4.0) { + gsl_sf_result result_c; + cheb_eval_e(&si_cs, (x*x-8.0)*0.125, &result_c); + result->val = x * (0.75 + result_c.val); + result->err = ax * result_c.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + /* Note there is no loss of precision + * here bcause of the leading constant. + */ + gsl_sf_result f; + gsl_sf_result g; + fg_asymp(ax, &f, &g); + result->val = 0.5 * M_PI - f.val*cos(ax) - g.val*sin(ax); + result->err = f.err + g.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + if(x < 0.0) result->val = -result->val; + return GSL_SUCCESS; + } +} + + +int gsl_sf_Ci_e(const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x <= 0.0) { + DOMAIN_ERROR(result); + } + else if(x <= 4.0) { + const double lx = log(x); + const double y = (x*x-8.0)*0.125; + gsl_sf_result result_c; + cheb_eval_e(&ci_cs, y, &result_c); + result->val = lx - 0.5 + result_c.val; + result->err = 2.0 * GSL_DBL_EPSILON * (fabs(lx) + 0.5) + result_c.err; + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + gsl_sf_result sin_result; + gsl_sf_result cos_result; + int stat_sin = gsl_sf_sin_e(x, &sin_result); + int stat_cos = gsl_sf_cos_e(x, &cos_result); + gsl_sf_result f; + gsl_sf_result g; + fg_asymp(x, &f, &g); + result->val = f.val*sin_result.val - g.val*cos_result.val; + result->err = fabs(f.err*sin_result.val); + result->err += fabs(g.err*cos_result.val); + result->err += fabs(f.val*sin_result.err); + result->err += fabs(g.val*cos_result.err); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_ERROR_SELECT_2(stat_sin, stat_cos); + } +} + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_Si(const double x) +{ + EVAL_RESULT(gsl_sf_Si_e(x, &result)); +} + +double gsl_sf_Ci(const double x) +{ + EVAL_RESULT(gsl_sf_Ci_e(x, &result)); +} |