3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
20 /* Author: G. Jungman */
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_sf_expint.h>
29 #include "chebyshev.h"
30 #include "cheb_eval.c"
32 static double expint3_data[24] = {
34 -0.248846446384140982,
35 0.80526220717231041e-01,
36 -0.25772733251968330e-01,
37 0.7599878873073774e-02,
38 -0.2030695581940405e-02,
39 0.490834586699330e-03,
40 -0.107682239142021e-03,
58 static cheb_series expint3_cs = {
65 static double expint3a_data[23] = {
66 1.9270464955068273729,
67 -0.349293565204813805e-01,
68 0.14503383718983009e-02,
69 -0.8925336718327903e-04,
90 static cheb_series expint3a_cs = {
98 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
100 int gsl_sf_expint_3_e(const double x, gsl_sf_result * result)
102 const double val_infinity = 0.892979511569249211;
104 /* CHECK_POINTER(result) */
107 DOMAIN_ERROR(result);
109 else if(x < 1.6*GSL_ROOT3_DBL_EPSILON) {
115 const double t = x*x*x/4.0 - 1.0;
116 gsl_sf_result result_c;
117 cheb_eval_e(&expint3_cs, t, &result_c);
118 result->val = x * result_c.val;
119 result->err = x * result_c.err;
122 else if(x < pow(-GSL_LOG_DBL_EPSILON, 1.0/3.0)) {
123 const double t = 16.0/(x*x*x) - 1.0;
124 const double s = exp(-x*x*x)/(3.0*x*x);
125 gsl_sf_result result_c;
126 cheb_eval_e(&expint3a_cs, t, &result_c);
127 result->val = val_infinity - result_c.val * s;
128 result->err = val_infinity * GSL_DBL_EPSILON + s * result_c.err;
132 result->val = val_infinity;
133 result->err = val_infinity * GSL_DBL_EPSILON;
139 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
143 double gsl_sf_expint_3(double x)
145 EVAL_RESULT(gsl_sf_expint_3_e(x, &result));