diff options
Diffstat (limited to 'gsl-1.9/specfunc/hyperg.c')
| -rw-r--r-- | gsl-1.9/specfunc/hyperg.c | 291 |
1 files changed, 291 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/hyperg.c b/gsl-1.9/specfunc/hyperg.c new file mode 100644 index 0000000..624e05e --- /dev/null +++ b/gsl-1.9/specfunc/hyperg.c @@ -0,0 +1,291 @@ +/* specfunc/hyperg.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 */ + +/* Miscellaneous implementations of use + * for evaluation of hypergeometric functions. + */ +#include <config.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_sf_exp.h> +#include <gsl/gsl_sf_gamma.h> + +#include "error.h" +#include "hyperg.h" + +#define SUM_LARGE (1.0e-5*GSL_DBL_MAX) + + +int +gsl_sf_hyperg_1F1_series_e(const double a, const double b, const double x, + gsl_sf_result * result + ) +{ + double an = a; + double bn = b; + double n = 1.0; + double del = 1.0; + double abs_del = 1.0; + double max_abs_del = 1.0; + double sum_val = 1.0; + double sum_err = 0.0; + + while(abs_del/fabs(sum_val) > 0.25*GSL_DBL_EPSILON) { + double u, abs_u; + + if(bn == 0.0) { + DOMAIN_ERROR(result); + } + + if(an == 0.0) { + result->val = sum_val; + result->err = sum_err; + result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val); + return GSL_SUCCESS; + } + + if (n > 10000.0) { + result->val = sum_val; + result->err = sum_err; + GSL_ERROR ("hypergeometric series failed to converge", GSL_EFAILED); + } + + u = x * (an/(bn*n)); + abs_u = fabs(u); + if(abs_u > 1.0 && max_abs_del > GSL_DBL_MAX/abs_u) { + result->val = sum_val; + result->err = fabs(sum_val); + GSL_ERROR ("overflow", GSL_EOVRFLW); + } + del *= u; + sum_val += del; + if(fabs(sum_val) > SUM_LARGE) { + result->val = sum_val; + result->err = fabs(sum_val); + GSL_ERROR ("overflow", GSL_EOVRFLW); + } + + abs_del = fabs(del); + max_abs_del = GSL_MAX_DBL(abs_del, max_abs_del); + sum_err += 2.0*GSL_DBL_EPSILON*abs_del; + + an += 1.0; + bn += 1.0; + n += 1.0; + } + + result->val = sum_val; + result->err = sum_err; + result->err += abs_del; + result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val); + + return GSL_SUCCESS; +} + + +int +gsl_sf_hyperg_1F1_large_b_e(const double a, const double b, const double x, gsl_sf_result * result) +{ + if(fabs(x/b) < 1.0) { + const double u = x/b; + const double v = 1.0/(1.0-u); + const double pre = pow(v,a); + const double uv = u*v; + const double uv2 = uv*uv; + const double t1 = a*(a+1.0)/(2.0*b)*uv2; + const double t2a = a*(a+1.0)/(24.0*b*b)*uv2; + const double t2b = 12.0 + 16.0*(a+2.0)*uv + 3.0*(a+2.0)*(a+3.0)*uv2; + const double t2 = t2a*t2b; + result->val = pre * (1.0 - t1 + t2); + result->err = pre * GSL_DBL_EPSILON * (1.0 + fabs(t1) + fabs(t2)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else { + DOMAIN_ERROR(result); + } +} + + +int +gsl_sf_hyperg_U_large_b_e(const double a, const double b, const double x, + gsl_sf_result * result, + double * ln_multiplier + ) +{ + double N = floor(b); /* b = N + eps */ + double eps = b - N; + + if(fabs(eps) < GSL_SQRT_DBL_EPSILON) { + double lnpre_val; + double lnpre_err; + gsl_sf_result M; + if(b > 1.0) { + double tmp = (1.0-b)*log(x); + gsl_sf_result lg_bm1; + gsl_sf_result lg_a; + gsl_sf_lngamma_e(b-1.0, &lg_bm1); + gsl_sf_lngamma_e(a, &lg_a); + lnpre_val = tmp + x + lg_bm1.val - lg_a.val; + lnpre_err = lg_bm1.err + lg_a.err + GSL_DBL_EPSILON * (fabs(x) + fabs(tmp)); + gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, -x, &M); + } + else { + gsl_sf_result lg_1mb; + gsl_sf_result lg_1pamb; + gsl_sf_lngamma_e(1.0-b, &lg_1mb); + gsl_sf_lngamma_e(1.0+a-b, &lg_1pamb); + lnpre_val = lg_1mb.val - lg_1pamb.val; + lnpre_err = lg_1mb.err + lg_1pamb.err; + gsl_sf_hyperg_1F1_large_b_e(a, b, x, &M); + } + + if(lnpre_val > GSL_LOG_DBL_MAX-10.0) { + result->val = M.val; + result->err = M.err; + *ln_multiplier = lnpre_val; + GSL_ERROR ("overflow", GSL_EOVRFLW); + } + else { + gsl_sf_result epre; + int stat_e = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &epre); + result->val = epre.val * M.val; + result->err = epre.val * M.err + epre.err * fabs(M.val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + *ln_multiplier = 0.0; + return stat_e; + } + } + else { + double omb_lnx = (1.0-b)*log(x); + gsl_sf_result lg_1mb; double sgn_1mb; + gsl_sf_result lg_1pamb; double sgn_1pamb; + gsl_sf_result lg_bm1; double sgn_bm1; + gsl_sf_result lg_a; double sgn_a; + gsl_sf_result M1, M2; + double lnpre1_val, lnpre2_val; + double lnpre1_err, lnpre2_err; + double sgpre1, sgpre2; + gsl_sf_hyperg_1F1_large_b_e( a, b, x, &M1); + gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, x, &M2); + + gsl_sf_lngamma_sgn_e(1.0-b, &lg_1mb, &sgn_1mb); + gsl_sf_lngamma_sgn_e(1.0+a-b, &lg_1pamb, &sgn_1pamb); + + gsl_sf_lngamma_sgn_e(b-1.0, &lg_bm1, &sgn_bm1); + gsl_sf_lngamma_sgn_e(a, &lg_a, &sgn_a); + + lnpre1_val = lg_1mb.val - lg_1pamb.val; + lnpre1_err = lg_1mb.err + lg_1pamb.err; + lnpre2_val = lg_bm1.val - lg_a.val - omb_lnx - x; + lnpre2_err = lg_bm1.err + lg_a.err + GSL_DBL_EPSILON * (fabs(omb_lnx)+fabs(x)); + sgpre1 = sgn_1mb * sgn_1pamb; + sgpre2 = sgn_bm1 * sgn_a; + + if(lnpre1_val > GSL_LOG_DBL_MAX-10.0 || lnpre2_val > GSL_LOG_DBL_MAX-10.0) { + double max_lnpre_val = GSL_MAX(lnpre1_val,lnpre2_val); + double max_lnpre_err = GSL_MAX(lnpre1_err,lnpre2_err); + double lp1 = lnpre1_val - max_lnpre_val; + double lp2 = lnpre2_val - max_lnpre_val; + double t1 = sgpre1*exp(lp1); + double t2 = sgpre2*exp(lp2); + result->val = t1*M1.val + t2*M2.val; + result->err = fabs(t1)*M1.err + fabs(t2)*M2.err; + result->err += GSL_DBL_EPSILON * exp(max_lnpre_err) * (fabs(t1*M1.val) + fabs(t2*M2.val)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + *ln_multiplier = max_lnpre_val; + GSL_ERROR ("overflow", GSL_EOVRFLW); + } + else { + double t1 = sgpre1*exp(lnpre1_val); + double t2 = sgpre2*exp(lnpre2_val); + result->val = t1*M1.val + t2*M2.val; + result->err = fabs(t1) * M1.err + fabs(t2)*M2.err; + result->err += GSL_DBL_EPSILON * (exp(lnpre1_err)*fabs(t1*M1.val) + exp(lnpre2_err)*fabs(t2*M2.val)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + *ln_multiplier = 0.0; + return GSL_SUCCESS; + } + } +} + + + +/* [Carlson, p.109] says the error in truncating this asymptotic series + * is less than the absolute value of the first neglected term. + * + * A termination argument is provided, so that the series will + * be summed at most up to n=n_trunc. If n_trunc is set negative, + * then the series is summed until it appears to start diverging. + */ +int +gsl_sf_hyperg_2F0_series_e(const double a, const double b, const double x, + int n_trunc, + gsl_sf_result * result + ) +{ + const int maxiter = 2000; + double an = a; + double bn = b; + double n = 1.0; + double sum = 1.0; + double del = 1.0; + double abs_del = 1.0; + double max_abs_del = 1.0; + double last_abs_del = 1.0; + + while(abs_del/fabs(sum) > GSL_DBL_EPSILON && n < maxiter) { + + double u = an * (bn/n * x); + double abs_u = fabs(u); + + if(abs_u > 1.0 && (max_abs_del > GSL_DBL_MAX/abs_u)) { + result->val = sum; + result->err = fabs(sum); + GSL_ERROR ("overflow", GSL_EOVRFLW); + } + + del *= u; + sum += del; + + abs_del = fabs(del); + + if(abs_del > last_abs_del) break; /* series is probably starting to grow */ + + last_abs_del = abs_del; + max_abs_del = GSL_MAX(abs_del, max_abs_del); + + an += 1.0; + bn += 1.0; + n += 1.0; + + if(an == 0.0 || bn == 0.0) break; /* series terminated */ + + if(n_trunc >= 0 && n >= n_trunc) break; /* reached requested timeout */ + } + + result->val = sum; + result->err = GSL_DBL_EPSILON * n + abs_del; + if(n >= maxiter) + GSL_ERROR ("error", GSL_EMAXITER); + else + return GSL_SUCCESS; +} |