| /* 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)); |
| } |
