diff options
Diffstat (limited to 'gsl-1.9/specfunc/clausen.c')
-rw-r--r-- | gsl-1.9/specfunc/clausen.c | 111 |
1 files changed, 111 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/clausen.c b/gsl-1.9/specfunc/clausen.c new file mode 100644 index 0000000..87632f1 --- /dev/null +++ b/gsl-1.9/specfunc/clausen.c @@ -0,0 +1,111 @@ +/* specfunc/clausen.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_clausen.h> + +#include "chebyshev.h" +#include "cheb_eval.c" + +static double aclaus_data[15] = { + 2.142694363766688447e+00, + 0.723324281221257925e-01, + 0.101642475021151164e-02, + 0.3245250328531645e-04, + 0.133315187571472e-05, + 0.6213240591653e-07, + 0.313004135337e-08, + 0.16635723056e-09, + 0.919659293e-11, + 0.52400462e-12, + 0.3058040e-13, + 0.18197e-14, + 0.1100e-15, + 0.68e-17, + 0.4e-18 +}; +static cheb_series aclaus_cs = { + aclaus_data, + 14, + -1, 1, + 8 /* FIXME: this is a guess, correct value needed here BJG */ +}; + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int gsl_sf_clausen_e(double x, gsl_sf_result *result) +{ + const double x_cut = M_PI * GSL_SQRT_DBL_EPSILON; + + double sgn = 1.0; + int status_red; + + if(x < 0.0) { + x = -x; + sgn = -1.0; + } + + /* Argument reduction to [0, 2pi) */ + status_red = gsl_sf_angle_restrict_pos_e(&x); + + /* Further reduction to [0,pi) */ + if(x > M_PI) { + /* simulated extra precision: 2PI = p0 + p1 */ + const double p0 = 6.28125; + const double p1 = 0.19353071795864769253e-02; + x = (p0 - x) + p1; + sgn = -sgn; + } + + if(x == 0.0) { + result->val = 0.0; + result->err = 0.0; + } + else if(x < x_cut) { + result->val = x * (1.0 - log(x)); + result->err = x * GSL_DBL_EPSILON; + } + else { + const double t = 2.0*(x*x / (M_PI*M_PI) - 0.5); + gsl_sf_result result_c; + cheb_eval_e(&aclaus_cs, t, &result_c); + result->val = x * (result_c.val - log(x)); + result->err = x * (result_c.err + GSL_DBL_EPSILON); + } + + result->val *= sgn; + + return status_red; +} + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_clausen(const double x) +{ + EVAL_RESULT(gsl_sf_clausen_e(x, &result)); +} |