diff options
Diffstat (limited to 'gsl-1.9/specfunc/hyperg_1F1.c')
| -rw-r--r-- | gsl-1.9/specfunc/hyperg_1F1.c | 2064 |
1 files changed, 2064 insertions, 0 deletions
diff --git a/gsl-1.9/specfunc/hyperg_1F1.c b/gsl-1.9/specfunc/hyperg_1F1.c new file mode 100644 index 0000000..15ffd12 --- /dev/null +++ b/gsl-1.9/specfunc/hyperg_1F1.c @@ -0,0 +1,2064 @@ +/* specfunc/hyperg_1F1.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_elementary.h> +#include <gsl/gsl_sf_exp.h> +#include <gsl/gsl_sf_bessel.h> +#include <gsl/gsl_sf_gamma.h> +#include <gsl/gsl_sf_laguerre.h> +#include <gsl/gsl_sf_hyperg.h> + +#include "error.h" +#include "hyperg.h" + +#define _1F1_INT_THRESHOLD (100.0*GSL_DBL_EPSILON) + + +/* Asymptotic result for 1F1(a, b, x) x -> -Infinity. + * Assumes b-a != neg integer and b != neg integer. + */ +static +int +hyperg_1F1_asymp_negx(const double a, const double b, const double x, + gsl_sf_result * result) +{ + gsl_sf_result lg_b; + gsl_sf_result lg_bma; + double sgn_b; + double sgn_bma; + + int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b); + int stat_bma = gsl_sf_lngamma_sgn_e(b-a, &lg_bma, &sgn_bma); + + if(stat_b == GSL_SUCCESS && stat_bma == GSL_SUCCESS) { + gsl_sf_result F; + int stat_F = gsl_sf_hyperg_2F0_series_e(a, 1.0+a-b, -1.0/x, -1, &F); + if(F.val != 0) { + double ln_term_val = a*log(-x); + double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(ln_term_val)); + double ln_pre_val = lg_b.val - lg_bma.val - ln_term_val; + double ln_pre_err = lg_b.err + lg_bma.err + ln_term_err; + int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err, + sgn_bma*sgn_b*F.val, F.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_F); + } + else { + result->val = 0.0; + result->err = 0.0; + return stat_F; + } + } + else { + DOMAIN_ERROR(result); + } +} + + +/* Asymptotic result for 1F1(a, b, x) x -> +Infinity + * Assumes b != neg integer and a != neg integer + */ +static +int +hyperg_1F1_asymp_posx(const double a, const double b, const double x, + gsl_sf_result * result) +{ + gsl_sf_result lg_b; + gsl_sf_result lg_a; + double sgn_b; + double sgn_a; + + int stat_b = gsl_sf_lngamma_sgn_e(b, &lg_b, &sgn_b); + int stat_a = gsl_sf_lngamma_sgn_e(a, &lg_a, &sgn_a); + + if(stat_a == GSL_SUCCESS && stat_b == GSL_SUCCESS) { + gsl_sf_result F; + int stat_F = gsl_sf_hyperg_2F0_series_e(b-a, 1.0-a, 1.0/x, -1, &F); + if(stat_F == GSL_SUCCESS && F.val != 0) { + double lnx = log(x); + double ln_term_val = (a-b)*lnx; + double ln_term_err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + fabs(b)) * fabs(lnx) + + 2.0 * GSL_DBL_EPSILON * fabs(a-b); + double ln_pre_val = lg_b.val - lg_a.val + ln_term_val + x; + double ln_pre_err = lg_b.err + lg_a.err + ln_term_err + 2.0 * GSL_DBL_EPSILON * fabs(x); + int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err, + sgn_a*sgn_b*F.val, F.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_F); + } + else { + result->val = 0.0; + result->err = 0.0; + return stat_F; + } + } + else { + DOMAIN_ERROR(result); + } +} + +/* Asymptotic result from Slater 4.3.7 + * + * To get the general series, write M(a,b,x) as + * + * M(a,b,x)=sum ((a)_n/(b)_n) (x^n / n!) + * + * and expand (b)_n in inverse powers of b as follows + * + * -log(1/(b)_n) = sum_(k=0)^(n-1) log(b+k) + * = n log(b) + sum_(k=0)^(n-1) log(1+k/b) + * + * Do a taylor expansion of the log in 1/b and sum the resulting terms + * using the standard algebraic formulas for finite sums of powers of + * k. This should then give + * + * M(a,b,x) = sum_(n=0)^(inf) (a_n/n!) (x/b)^n * (1 - n(n-1)/(2b) + * + (n-1)n(n+1)(3n-2)/(24b^2) + ... + * + * which can be summed explicitly. The trick for summing it is to take + * derivatives of sum_(i=0)^(inf) a_n*y^n/n! = (1-y)^(-a); + * + * [BJG 16/01/2007] + */ + +static +int +hyperg_1F1_largebx(const double a, const double b, const double x, gsl_sf_result * result) +{ + double y = x/b; + double f = exp(-a*log1p(-y)); + double t1 = -((a*(a+1.0))/(2*b))*pow((y/(1.0-y)),2.0); + double t2 = (1/(24*b*b))*((a*(a+1)*y*y)/pow(1-y,4))*(12+8*(2*a+1)*y+(3*a*a-a-2)*y*y); + double t3 = (-1/(48*b*b*b*pow(1-y,6)))*a*((a + 1)*((y*((a + 1)*(a*(y*(y*((y*(a - 2) + 16)*(a - 1)) + 72)) + 96)) + 24)*pow(y, 2))); + result->val = f * (1 + t1 + t2 + t3); + result->err = 2*fabs(f*t3) + 2*GSL_DBL_EPSILON*fabs(result->val); + return GSL_SUCCESS; +} + +/* Asymptotic result for x < 2b-4a, 2b-4a large. + * [Abramowitz+Stegun, 13.5.21] + * + * assumes 0 <= x/(2b-4a) <= 1 + */ +static +int +hyperg_1F1_large2bm4a(const double a, const double b, const double x, gsl_sf_result * result) +{ + double eta = 2.0*b - 4.0*a; + double cos2th = x/eta; + double sin2th = 1.0 - cos2th; + double th = acos(sqrt(cos2th)); + double pre_h = 0.25*M_PI*M_PI*eta*eta*cos2th*sin2th; + gsl_sf_result lg_b; + int stat_lg = gsl_sf_lngamma_e(b, &lg_b); + double t1 = 0.5*(1.0-b)*log(0.25*x*eta); + double t2 = 0.25*log(pre_h); + double lnpre_val = lg_b.val + 0.5*x + t1 - t2; + double lnpre_err = lg_b.err + 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + fabs(t1) + fabs(t2)); +#if SMALL_ANGLE + const double eps = asin(sqrt(cos2th)); /* theta = pi/2 - eps */ + double s1 = (fmod(a, 1.0) == 0.0) ? 0.0 : sin(a*M_PI); + double eta_reduc = (fmod(eta + 1, 4.0) == 0.0) ? 0.0 : fmod(eta + 1, 8.0); + double phi1 = 0.25*eta_reduc*M_PI; + double phi2 = 0.25*eta*(2*eps + sin(2.0*eps)); + double s2 = sin(phi1 - phi2); +#else + double s1 = sin(a*M_PI); + double s2 = sin(0.25*eta*(2.0*th - sin(2.0*th)) + 0.25*M_PI); +#endif + double ser_val = s1 + s2; + double ser_err = 2.0 * GSL_DBL_EPSILON * (fabs(s1) + fabs(s2)); + int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err, + ser_val, ser_err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_lg); +} + + +/* Luke's rational approximation. + * See [Luke, Algorithms for the Computation of Mathematical Functions, p.182] + * + * Like the case of the 2F1 rational approximations, these are + * probably guaranteed to converge for x < 0, barring gross + * numerical instability in the pre-asymptotic regime. + */ +static +int +hyperg_1F1_luke(const double a, const double c, const double xin, + gsl_sf_result * result) +{ + const double RECUR_BIG = 1.0e+50; + const int nmax = 5000; + int n = 3; + const double x = -xin; + const double x3 = x*x*x; + const double t0 = a/c; + const double t1 = (a+1.0)/(2.0*c); + const double t2 = (a+2.0)/(2.0*(c+1.0)); + double F = 1.0; + double prec; + + double Bnm3 = 1.0; /* B0 */ + double Bnm2 = 1.0 + t1 * x; /* B1 */ + double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x); /* B2 */ + + double Anm3 = 1.0; /* A0 */ + double Anm2 = Bnm2 - t0 * x; /* A1 */ + double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x; /* A2 */ + + while(1) { + double npam1 = n + a - 1; + double npcm1 = n + c - 1; + double npam2 = n + a - 2; + double npcm2 = n + c - 2; + double tnm1 = 2*n - 1; + double tnm3 = 2*n - 3; + double tnm5 = 2*n - 5; + double F1 = (n-a-2) / (2*tnm3*npcm1); + double F2 = (n+a)*npam1 / (4*tnm1*tnm3*npcm2*npcm1); + double F3 = -npam2*npam1*(n-a-2) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1); + double E = -npam1*(n-c-1) / (2*tnm3*npcm2*npcm1); + + double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3; + double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3; + double r = An/Bn; + + prec = fabs((F - r)/F); + F = r; + + if(prec < GSL_DBL_EPSILON || n > nmax) break; + + if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) { + An /= RECUR_BIG; + Bn /= RECUR_BIG; + Anm1 /= RECUR_BIG; + Bnm1 /= RECUR_BIG; + Anm2 /= RECUR_BIG; + Bnm2 /= RECUR_BIG; + Anm3 /= RECUR_BIG; + Bnm3 /= RECUR_BIG; + } + else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) { + An *= RECUR_BIG; + Bn *= RECUR_BIG; + Anm1 *= RECUR_BIG; + Bnm1 *= RECUR_BIG; + Anm2 *= RECUR_BIG; + Bnm2 *= RECUR_BIG; + Anm3 *= RECUR_BIG; + Bnm3 *= RECUR_BIG; + } + + n++; + Bnm3 = Bnm2; + Bnm2 = Bnm1; + Bnm1 = Bn; + Anm3 = Anm2; + Anm2 = Anm1; + Anm1 = An; + } + + result->val = F; + result->err = 2.0 * fabs(F * prec); + result->err += 2.0 * GSL_DBL_EPSILON * (n-1.0) * fabs(F); + + return GSL_SUCCESS; +} + +/* Series for 1F1(1,b,x) + * b > 0 + */ +static +int +hyperg_1F1_1_series(const double b, const double x, gsl_sf_result * result) +{ + double sum_val = 1.0; + double sum_err = 0.0; + double term = 1.0; + double n = 1.0; + while(fabs(term/sum_val) > 0.25*GSL_DBL_EPSILON) { + term *= x/(b+n-1); + sum_val += term; + sum_err += 8.0*GSL_DBL_EPSILON*fabs(term) + GSL_DBL_EPSILON*fabs(sum_val); + n += 1.0; + } + result->val = sum_val; + result->err = sum_err; + result->err += 2.0 * fabs(term); + return GSL_SUCCESS; +} + + +/* 1F1(1,b,x) + * b >= 1, b integer + */ +static +int +hyperg_1F1_1_int(const int b, const double x, gsl_sf_result * result) +{ + if(b < 1) { + DOMAIN_ERROR(result); + } + else if(b == 1) { + return gsl_sf_exp_e(x, result); + } + else if(b == 2) { + return gsl_sf_exprel_e(x, result); + } + else if(b == 3) { + return gsl_sf_exprel_2_e(x, result); + } + else { + return gsl_sf_exprel_n_e(b-1, x, result); + } +} + + +/* 1F1(1,b,x) + * b >=1, b real + * + * checked OK: [GJ] Thu Oct 1 16:46:35 MDT 1998 + */ +static +int +hyperg_1F1_1(const double b, const double x, gsl_sf_result * result) +{ + double ax = fabs(x); + double ib = floor(b + 0.1); + + if(b < 1.0) { + DOMAIN_ERROR(result); + } + else if(b == 1.0) { + return gsl_sf_exp_e(x, result); + } + else if(b >= 1.4*ax) { + return hyperg_1F1_1_series(b, x, result); + } + else if(fabs(b - ib) < _1F1_INT_THRESHOLD && ib < INT_MAX) { + return hyperg_1F1_1_int((int)ib, x, result); + } + else if(x > 0.0) { + if(x > 100.0 && b < 0.75*x) { + return hyperg_1F1_asymp_posx(1.0, b, x, result); + } + else if(b < 1.0e+05) { + /* Recurse backward on b, from a + * chosen offset point. For x > 0, + * which holds here, this should + * be a stable direction. + */ + const double off = ceil(1.4*x-b) + 1.0; + double bp = b + off; + gsl_sf_result M; + int stat_s = hyperg_1F1_1_series(bp, x, &M); + const double err_rat = M.err / fabs(M.val); + while(bp > b+0.1) { + /* M(1,b-1) = x/(b-1) M(1,b) + 1 */ + bp -= 1.0; + M.val = 1.0 + x/bp * M.val; + } + result->val = M.val; + result->err = err_rat * fabs(M.val); + result->err += 2.0 * GSL_DBL_EPSILON * (fabs(off)+1.0) * fabs(M.val); + return stat_s; + } else if (fabs(x) < fabs(b) && fabs(x) < sqrt(fabs(b)) * fabs(b-x)) { + return hyperg_1F1_largebx(1.0, b, x, result); + } else if (fabs(x) > fabs(b)) { + return hyperg_1F1_1_series(b, x, result); + } else { + return hyperg_1F1_large2bm4a(1.0, b, x, result); + } + } + else { + /* x <= 0 and b not large compared to |x| + */ + if(ax < 10.0 && b < 10.0) { + return hyperg_1F1_1_series(b, x, result); + } + else if(ax >= 100.0 && GSL_MAX_DBL(fabs(2.0-b),1.0) < 0.99*ax) { + return hyperg_1F1_asymp_negx(1.0, b, x, result); + } + else { + return hyperg_1F1_luke(1.0, b, x, result); + } + } +} + + +/* 1F1(a,b,x)/Gamma(b) for b->0 + * [limit of Abramowitz+Stegun 13.3.7] + */ +static +int +hyperg_1F1_renorm_b0(const double a, const double x, gsl_sf_result * result) +{ + double eta = a*x; + if(eta > 0.0) { + double root_eta = sqrt(eta); + gsl_sf_result I1_scaled; + int stat_I = gsl_sf_bessel_I1_scaled_e(2.0*root_eta, &I1_scaled); + if(I1_scaled.val <= 0.0) { + result->val = 0.0; + result->err = 0.0; + return GSL_ERROR_SELECT_2(stat_I, GSL_EDOM); + } + else { + /* Note that 13.3.7 contains higher terms which are zeroth order + in b. These make a non-negligible contribution to the sum. + With the first correction term, the I1 above is replaced by + I1 + (2/3)*a*(x/(4a))**(3/2)*I2(2*root_eta). We will add + this as part of the result and error estimate. */ + + const double corr1 =(2.0/3.0)*a*pow(x/(4.0*a),1.5)*gsl_sf_bessel_In_scaled(2, 2.0*root_eta) + ; + const double lnr_val = 0.5*x + 0.5*log(eta) + fabs(2.0*root_eta) + log(I1_scaled.val+corr1); + const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs((I1_scaled.err+corr1)/I1_scaled.val); + return gsl_sf_exp_err_e(lnr_val, lnr_err, result); + } + } + else if(eta == 0.0) { + result->val = 0.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else { + /* eta < 0 */ + double root_eta = sqrt(-eta); + gsl_sf_result J1; + int stat_J = gsl_sf_bessel_J1_e(2.0*root_eta, &J1); + if(J1.val <= 0.0) { + result->val = 0.0; + result->err = 0.0; + return GSL_ERROR_SELECT_2(stat_J, GSL_EDOM); + } + else { + const double t1 = 0.5*x; + const double t2 = 0.5*log(-eta); + const double t3 = fabs(x); + const double t4 = log(J1.val); + const double lnr_val = t1 + t2 + t3 + t4; + const double lnr_err = GSL_DBL_EPSILON * (1.5*fabs(x) + 1.0) + fabs(J1.err/J1.val); + gsl_sf_result ex; + int stat_e = gsl_sf_exp_err_e(lnr_val, lnr_err, &ex); + result->val = -ex.val; + result->err = ex.err; + return stat_e; + } + } + +} + + +/* 1F1'(a,b,x)/1F1(a,b,x) + * Uses Gautschi's version of the CF. + * [Gautschi, Math. Comp. 31, 994 (1977)] + * + * Supposedly this suffers from the "anomalous convergence" + * problem when b < x. I have seen anomalous convergence + * in several of the continued fractions associated with + * 1F1(a,b,x). This particular CF formulation seems stable + * for b > x. However, it does display a painful artifact + * of the anomalous convergence; the convergence plateaus + * unless b >>> x. For example, even for b=1000, x=1, this + * method locks onto a ratio which is only good to about + * 4 digits. Apparently the rest of the digits are hiding + * way out on the plateau, but finite-precision lossage + * means you will never get them. + */ +#if 0 +static +int +hyperg_1F1_CF1_p(const double a, const double b, const double x, double * result) +{ + const double RECUR_BIG = GSL_SQRT_DBL_MAX; + const int maxiter = 5000; + int n = 1; + double Anm2 = 1.0; + double Bnm2 = 0.0; + double Anm1 = 0.0; + double Bnm1 = 1.0; + double a1 = 1.0; + double b1 = 1.0; + double An = b1*Anm1 + a1*Anm2; + double Bn = b1*Bnm1 + a1*Bnm2; + double an, bn; + double fn = An/Bn; + + while(n < maxiter) { + double old_fn; + double del; + n++; + Anm2 = Anm1; + Bnm2 = Bnm1; + Anm1 = An; + Bnm1 = Bn; + an = (a+n)*x/((b-x+n-1)*(b-x+n)); + bn = 1.0; + An = bn*Anm1 + an*Anm2; + Bn = bn*Bnm1 + an*Bnm2; + + if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) { + An /= RECUR_BIG; + Bn /= RECUR_BIG; + Anm1 /= RECUR_BIG; + Bnm1 /= RECUR_BIG; + Anm2 /= RECUR_BIG; + Bnm2 /= RECUR_BIG; + } + + old_fn = fn; + fn = An/Bn; + del = old_fn/fn; + + if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break; + } + + *result = a/(b-x) * fn; + + if(n == maxiter) + GSL_ERROR ("error", GSL_EMAXITER); + else + return GSL_SUCCESS; +} +#endif /* 0 */ + + +/* 1F1'(a,b,x)/1F1(a,b,x) + * Uses Gautschi's series transformation of the + * continued fraction. This is apparently the best + * method for getting this ratio in the stable region. + * The convergence is monotone and supergeometric + * when b > x. + * Assumes a >= -1. + */ +static +int +hyperg_1F1_CF1_p_ser(const double a, const double b, const double x, double * result) +{ + if(a == 0.0) { + *result = 0.0; + return GSL_SUCCESS; + } + else { + const int maxiter = 5000; + double sum = 1.0; + double pk = 1.0; + double rhok = 0.0; + int k; + for(k=1; k<maxiter; k++) { + double ak = (a + k)*x/((b-x+k-1.0)*(b-x+k)); + rhok = -ak*(1.0 + rhok)/(1.0 + ak*(1.0+rhok)); + pk *= rhok; + sum += pk; + if(fabs(pk/sum) < 2.0*GSL_DBL_EPSILON) break; + } + *result = a/(b-x) * sum; + if(k == maxiter) + GSL_ERROR ("error", GSL_EMAXITER); + else + return GSL_SUCCESS; + } +} + + +/* 1F1(a+1,b,x)/1F1(a,b,x) + * + * I think this suffers from typical "anomalous convergence". + * I could not find a region where it was truly useful. + */ +#if 0 +static +int +hyperg_1F1_CF1(const double a, const double b, const double x, double * result) +{ + const double RECUR_BIG = GSL_SQRT_DBL_MAX; + const int maxiter = 5000; + int n = 1; + double Anm2 = 1.0; + double Bnm2 = 0.0; + double Anm1 = 0.0; + double Bnm1 = 1.0; + double a1 = b - a - 1.0; + double b1 = b - x - 2.0*(a+1.0); + double An = b1*Anm1 + a1*Anm2; + double Bn = b1*Bnm1 + a1*Bnm2; + double an, bn; + double fn = An/Bn; + + while(n < maxiter) { + double old_fn; + double del; + n++; + Anm2 = Anm1; + Bnm2 = Bnm1; + Anm1 = An; + Bnm1 = Bn; + an = (a + n - 1.0) * (b - a - n); + bn = b - x - 2.0*(a+n); + An = bn*Anm1 + an*Anm2; + Bn = bn*Bnm1 + an*Bnm2; + + if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) { + An /= RECUR_BIG; + Bn /= RECUR_BIG; + Anm1 /= RECUR_BIG; + Bnm1 /= RECUR_BIG; + Anm2 /= RECUR_BIG; + Bnm2 /= RECUR_BIG; + } + + old_fn = fn; + fn = An/Bn; + del = old_fn/fn; + + if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break; + } + + *result = fn; + if(n == maxiter) + GSL_ERROR ("error", GSL_EMAXITER); + else + return GSL_SUCCESS; +} +#endif /* 0 */ + + +/* 1F1(a,b+1,x)/1F1(a,b,x) + * + * This seemed to suffer from "anomalous convergence". + * However, I have no theory for this recurrence. + */ +#if 0 +static +int +hyperg_1F1_CF1_b(const double a, const double b, const double x, double * result) +{ + const double RECUR_BIG = GSL_SQRT_DBL_MAX; + const int maxiter = 5000; + int n = 1; + double Anm2 = 1.0; + double Bnm2 = 0.0; + double Anm1 = 0.0; + double Bnm1 = 1.0; + double a1 = b + 1.0; + double b1 = (b + 1.0) * (b - x); + double An = b1*Anm1 + a1*Anm2; + double Bn = b1*Bnm1 + a1*Bnm2; + double an, bn; + double fn = An/Bn; + + while(n < maxiter) { + double old_fn; + double del; + n++; + Anm2 = Anm1; + Bnm2 = Bnm1; + Anm1 = An; + Bnm1 = Bn; + an = (b + n) * (b + n - 1.0 - a) * x; + bn = (b + n) * (b + n - 1.0 - x); + An = bn*Anm1 + an*Anm2; + Bn = bn*Bnm1 + an*Bnm2; + + if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) { + An /= RECUR_BIG; + Bn /= RECUR_BIG; + Anm1 /= RECUR_BIG; + Bnm1 /= RECUR_BIG; + Anm2 /= RECUR_BIG; + Bnm2 /= RECUR_BIG; + } + + old_fn = fn; + fn = An/Bn; + del = old_fn/fn; + + if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break; + } + + *result = fn; + if(n == maxiter) + GSL_ERROR ("error", GSL_EMAXITER); + else + return GSL_SUCCESS; +} +#endif /* 0 */ + + +/* 1F1(a,b,x) + * |a| <= 1, b > 0 + */ +static +int +hyperg_1F1_small_a_bgt0(const double a, const double b, const double x, gsl_sf_result * result) +{ + const double bma = b-a; + const double oma = 1.0-a; + const double ap1mb = 1.0+a-b; + const double abs_bma = fabs(bma); + const double abs_oma = fabs(oma); + const double abs_ap1mb = fabs(ap1mb); + + const double ax = fabs(x); + + if(a == 0.0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(a == 1.0 && b >= 1.0) { + return hyperg_1F1_1(b, x, result); + } + else if(a == -1.0) { + result->val = 1.0 + a/b * x; + result->err = GSL_DBL_EPSILON * (1.0 + fabs(a/b * x)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return GSL_SUCCESS; + } + else if(b >= 1.4*ax) { + return gsl_sf_hyperg_1F1_series_e(a, b, x, result); + } + else if(x > 0.0) { + if(x > 100.0 && abs_bma*abs_oma < 0.5*x) { + return hyperg_1F1_asymp_posx(a, b, x, result); + } + else if(b < 5.0e+06) { + /* Recurse backward on b from + * a suitably high point. + */ + const double b_del = ceil(1.4*x-b) + 1.0; + double bp = b + b_del; + gsl_sf_result r_Mbp1; + gsl_sf_result r_Mb; + double Mbp1; + double Mb; + double Mbm1; + int stat_0 = gsl_sf_hyperg_1F1_series_e(a, bp+1.0, x, &r_Mbp1); + int stat_1 = gsl_sf_hyperg_1F1_series_e(a, bp, x, &r_Mb); + const double err_rat = fabs(r_Mbp1.err/r_Mbp1.val) + fabs(r_Mb.err/r_Mb.val); + Mbp1 = r_Mbp1.val; + Mb = r_Mb.val; + while(bp > b+0.1) { + /* Do backward recursion. */ + Mbm1 = ((x+bp-1.0)*Mb - x*(bp-a)/bp*Mbp1)/(bp-1.0); + bp -= 1.0; + Mbp1 = Mb; + Mb = Mbm1; + } + result->val = Mb; + result->err = err_rat * (fabs(b_del)+1.0) * fabs(Mb); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mb); + return GSL_ERROR_SELECT_2(stat_0, stat_1); + } + else if (fabs(x) < fabs(b) && fabs(a*x) < sqrt(fabs(b)) * fabs(b-x)) { + return hyperg_1F1_largebx(a, b, x, result); + } else { + return hyperg_1F1_large2bm4a(a, b, x, result); + } + } + else { + /* x < 0 and b not large compared to |x| + */ + if(ax < 10.0 && b < 10.0) { + return gsl_sf_hyperg_1F1_series_e(a, b, x, result); + } + else if(ax >= 100.0 && GSL_MAX(abs_ap1mb,1.0) < 0.99*ax) { + return hyperg_1F1_asymp_negx(a, b, x, result); + } + else { + return hyperg_1F1_luke(a, b, x, result); + } + } +} + + +/* 1F1(b+eps,b,x) + * |eps|<=1, b > 0 + */ +static +int +hyperg_1F1_beps_bgt0(const double eps, const double b, const double x, gsl_sf_result * result) +{ + if(b > fabs(x) && fabs(eps) < GSL_SQRT_DBL_EPSILON) { + /* If b-a is very small and x/b is not too large we can + * use this explicit approximation. + * + * 1F1(b+eps,b,x) = exp(ax/b) (1 - eps x^2 (v2 + v3 x + ...) + ...) + * + * v2 = a/(2b^2(b+1)) + * v3 = a(b-2a)/(3b^3(b+1)(b+2)) + * ... + * + * See [Luke, Mathematical Functions and Their Approximations, p.292] + * + * This cannot be used for b near a negative integer or zero. + * Also, if x/b is large the deviation from exp(x) behaviour grows. + */ + double a = b + eps; + gsl_sf_result exab; + int stat_e = gsl_sf_exp_e(a*x/b, &exab); + double v2 = a/(2.0*b*b*(b+1.0)); + double v3 = a*(b-2.0*a)/(3.0*b*b*b*(b+1.0)*(b+2.0)); + double v = v2 + v3 * x; + double f = (1.0 - eps*x*x*v); + result->val = exab.val * f; + result->err = exab.err * fabs(f); + result->err += fabs(exab.val) * GSL_DBL_EPSILON * (1.0 + fabs(eps*x*x*v)); + result->err += 4.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_e; + } + else { + /* Otherwise use a Kummer transformation to reduce + * it to the small a case. + */ + gsl_sf_result Kummer_1F1; + int stat_K = hyperg_1F1_small_a_bgt0(-eps, b, -x, &Kummer_1F1); + if(Kummer_1F1.val != 0.0) { + int stat_e = gsl_sf_exp_mult_err_e(x, 2.0*GSL_DBL_EPSILON*fabs(x), + Kummer_1F1.val, Kummer_1F1.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } + else { + result->val = 0.0; + result->err = 0.0; + return stat_K; + } + } +} + + +/* 1F1(a,2a,x) = Gamma(a + 1/2) E(x) (|x|/4)^(-a+1/2) scaled_I(a-1/2,|x|/2) + * + * E(x) = exp(x) x > 0 + * = 1 x < 0 + * + * a >= 1/2 + */ +static +int +hyperg_1F1_beq2a_pos(const double a, const double x, gsl_sf_result * result) +{ + if(x == 0.0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else { + gsl_sf_result I; + int stat_I = gsl_sf_bessel_Inu_scaled_e(a-0.5, 0.5*fabs(x), &I); + gsl_sf_result lg; + int stat_g = gsl_sf_lngamma_e(a + 0.5, &lg); + double ln_term = (0.5-a)*log(0.25*fabs(x)); + double lnpre_val = lg.val + GSL_MAX_DBL(x,0.0) + ln_term; + double lnpre_err = lg.err + GSL_DBL_EPSILON * (fabs(ln_term) + fabs(x)); + int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err, + I.val, I.err, + result); + return GSL_ERROR_SELECT_3(stat_e, stat_g, stat_I); + } +} + + +/* Determine middle parts of diagonal recursion along b=2a + * from two endpoints, i.e. + * + * given: M(a,b) and M(a+1,b+2) + * get: M(a+1,b+1) and M(a,b+1) + */ +#if 0 +inline +static +int +hyperg_1F1_diag_step(const double a, const double b, const double x, + const double Mab, const double Map1bp2, + double * Map1bp1, double * Mabp1) +{ + if(a == b) { + *Map1bp1 = Mab; + *Mabp1 = Mab - x/(b+1.0) * Map1bp2; + } + else { + *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2; + *Mabp1 = (a * *Map1bp1 - b * Mab)/(a-b); + } + return GSL_SUCCESS; +} +#endif /* 0 */ + + +/* Determine endpoint of diagonal recursion. + * + * given: M(a,b) and M(a+1,b+2) + * get: M(a+1,b) and M(a+1,b+1) + */ +#if 0 +inline +static +int +hyperg_1F1_diag_end_step(const double a, const double b, const double x, + const double Mab, const double Map1bp2, + double * Map1b, double * Map1bp1) +{ + *Map1bp1 = Mab - x * (a-b)/(b*(b+1.0)) * Map1bp2; + *Map1b = Mab + x/b * *Map1bp1; + return GSL_SUCCESS; +} +#endif /* 0 */ + + +/* Handle the case of a and b both positive integers. + * Assumes a > 0 and b > 0. + */ +static +int +hyperg_1F1_ab_posint(const int a, const int b, const double x, gsl_sf_result * result) +{ + double ax = fabs(x); + + if(a == b) { + return gsl_sf_exp_e(x, result); /* 1F1(a,a,x) */ + } + else if(a == 1) { + return gsl_sf_exprel_n_e(b-1, x, result); /* 1F1(1,b,x) */ + } + else if(b == a + 1) { + gsl_sf_result K; + int stat_K = gsl_sf_exprel_n_e(a, -x, &K); /* 1F1(1,1+a,-x) */ + int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x), + K.val, K.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } + else if(a == b + 1) { + gsl_sf_result ex; + int stat_e = gsl_sf_exp_e(x, &ex); + result->val = ex.val * (1.0 + x/b); + result->err = ex.err * (1.0 + x/b); + result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b)); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_e; + } + else if(a == b + 2) { + gsl_sf_result ex; + int stat_e = gsl_sf_exp_e(x, &ex); + double poly = (1.0 + x/b*(2.0 + x/(b+1.0))); + result->val = ex.val * poly; + result->err = ex.err * fabs(poly); + result->err += ex.val * GSL_DBL_EPSILON * (1.0 + fabs(x/b) * (2.0 + fabs(x/(b+1.0)))); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); + return stat_e; + } + else if(b == 2*a) { + return hyperg_1F1_beq2a_pos(a, x, result); /* 1F1(a,2a,x) */ + } + else if( ( b < 10 && a < 10 && ax < 5.0 ) + || ( b > a*ax ) + || ( b > a && ax < 5.0 ) + ) { + return gsl_sf_hyperg_1F1_series_e(a, b, x, result); + } + else if(b > a && b >= 2*a + x) { + /* Use the Gautschi CF series, then + * recurse backward to a=0 for normalization. + * This will work for either sign of x. + */ + double rap; + int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap); + double ra = 1.0 + x/a * rap; + double Ma = GSL_SQRT_DBL_MIN; + double Map1 = ra * Ma; + double Mnp1 = Map1; + double Mn = Ma; + double Mnm1; + int n; + for(n=a; n>0; n--) { + Mnm1 = (n * Mnp1 - (2*n-b+x) * Mn) / (b-n); + Mnp1 = Mn; + Mn = Mnm1; + } + result->val = Ma/Mn; + result->err = 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ma/Mn); + return stat_CF1; + } + else if(b > a && b < 2*a + x && b > x) { + /* Use the Gautschi series representation of + * the continued fraction. Then recurse forward + * to the a=b line for normalization. This will + * work for either sign of x, although we do need + * to check for b > x, for when x is positive. + */ + double rap; + int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap); + double ra = 1.0 + x/a * rap; + gsl_sf_result ex; + int stat_ex; + + double Ma = GSL_SQRT_DBL_MIN; + double Map1 = ra * Ma; + double Mnm1 = Ma; + double Mn = Map1; + double Mnp1; + int n; + for(n=a+1; n<b; n++) { + Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n; + Mnm1 = Mn; + Mn = Mnp1; + } + + stat_ex = gsl_sf_exp_e(x, &ex); /* 1F1(b,b,x) */ + result->val = ex.val * Ma/Mn; + result->err = ex.err * fabs(Ma/Mn); + result->err += 4.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val); + return GSL_ERROR_SELECT_2(stat_ex, stat_CF1); + } + else if(x >= 0.0) { + + if(b < a) { + /* The point b,b is below the b=2a+x line. + * Forward recursion on a from b,b+1 is possible. + * Note that a > b + 1 as well, since we already tried a = b + 1. + */ + if(x + log(fabs(x/b)) < GSL_LOG_DBL_MAX-2.0) { + double ex = exp(x); + int n; + double Mnm1 = ex; /* 1F1(b,b,x) */ + double Mn = ex * (1.0 + x/b); /* 1F1(b+1,b,x) */ + double Mnp1; + for(n=b+1; n<a; n++) { + Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n; + Mnm1 = Mn; + Mn = Mnp1; + } + result->val = Mn; + result->err = (x + 1.0) * GSL_DBL_EPSILON * fabs(Mn); + result->err *= fabs(a-b)+1.0; + return GSL_SUCCESS; + } + else { + OVERFLOW_ERROR(result); + } + } + else { + /* b > a + * b < 2a + x + * b <= x (otherwise we would have finished above) + * + * Gautschi anomalous convergence region. However, we can + * recurse forward all the way from a=0,1 because we are + * always underneath the b=2a+x line. + */ + gsl_sf_result r_Mn; + double Mnm1 = 1.0; /* 1F1(0,b,x) */ + double Mn; /* 1F1(1,b,x) */ + double Mnp1; + int n; + gsl_sf_exprel_n_e(b-1, x, &r_Mn); + Mn = r_Mn.val; + for(n=1; n<a; n++) { + Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n; + Mnm1 = Mn; + Mn = Mnp1; + } + result->val = Mn; + result->err = fabs(Mn) * (1.0 + fabs(a)) * fabs(r_Mn.err / r_Mn.val); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn); + return GSL_SUCCESS; + } + } + else { + /* x < 0 + * b < a (otherwise we would have tripped one of the above) + */ + + if(a <= 0.5*(b-x) || a >= -x) { + /* Gautschi continued fraction is in the anomalous region, + * so we must find another way. We recurse down in b, + * from the a=b line. + */ + double ex = exp(x); + double Manp1 = ex; + double Man = ex * (1.0 + x/(a-1.0)); + double Manm1; + int n; + for(n=a-1; n>b; n--) { + Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0)); + Manp1 = Man; + Man = Manm1; + } + result->val = Man; + result->err = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Man); + result->err *= fabs(b-a)+1.0; + return GSL_SUCCESS; + } + else { + /* Pick a0 such that b ~= 2a0 + x, then + * recurse down in b from a0,a0 to determine + * the values near the line b=2a+x. Then recurse + * forward on a from a0. + */ + int a0 = ceil(0.5*(b-x)); + double Ma0b; /* M(a0,b) */ + double Ma0bp1; /* M(a0,b+1) */ + double Ma0p1b; /* M(a0+1,b) */ + double Mnm1; + double Mn; + double Mnp1; + int n; + { + double ex = exp(x); + double Ma0np1 = ex; + double Ma0n = ex * (1.0 + x/(a0-1.0)); + double Ma0nm1; + for(n=a0-1; n>b; n--) { + Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0)); + Ma0np1 = Ma0n; + Ma0n = Ma0nm1; + } + Ma0bp1 = Ma0np1; + Ma0b = Ma0n; + Ma0p1b = (b*(a0+x)*Ma0b + x*(a0-b)*Ma0bp1)/(a0*b); + } + + /* Initialise the recurrence correctly BJG */ + + if (a0 >= a) + { + Mn = Ma0b; + } + else if (a0 + 1>= a) + { + Mn = Ma0p1b; + } + else + { + Mnm1 = Ma0b; + Mn = Ma0p1b; + + for(n=a0+1; n<a; n++) { + Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n; + Mnm1 = Mn; + Mn = Mnp1; + } + } + + result->val = Mn; + result->err = (fabs(x) + 1.0) * GSL_DBL_EPSILON * fabs(Mn); + result->err *= fabs(b-a)+1.0; + return GSL_SUCCESS; + } + } +} + + +/* Evaluate a <= 0, a integer, cases directly. (Polynomial; Horner) + * When the terms are all positive, this + * must work. We will assume this here. + */ +static +int +hyperg_1F1_a_negint_poly(const int a, const double b, const double x, gsl_sf_result * result) +{ + if(a == 0) { + result->val = 1.0; + result->err = 1.0; + return GSL_SUCCESS; + } + else { + int N = -a; + double poly = 1.0; + int k; + for(k=N-1; k>=0; k--) { + double t = (a+k)/(b+k) * (x/(k+1)); + double r = t + 1.0/poly; + if(r > 0.9*GSL_DBL_MAX/poly) { + OVERFLOW_ERROR(result); + } + else { + poly *= r; /* P_n = 1 + t_n P_{n-1} */ + } + } + result->val = poly; + result->err = 2.0 * (sqrt(N) + 1.0) * GSL_DBL_EPSILON * fabs(poly); + return GSL_SUCCESS; + } +} + + +/* Evaluate negative integer a case by relation + * to Laguerre polynomials. This is more general than + * the direct polynomial evaluation, but is safe + * for all values of x. + * + * 1F1(-n,b,x) = n!/(b)_n Laguerre[n,b-1,x] + * = n B(b,n) Laguerre[n,b-1,x] + * + * assumes b is not a negative integer + */ +static +int +hyperg_1F1_a_negint_lag(const int a, const double b, const double x, gsl_sf_result * result) +{ + const int n = -a; + + gsl_sf_result lag; + const int stat_l = gsl_sf_laguerre_n_e(n, b-1.0, x, &lag); + if(b < 0.0) { + gsl_sf_result lnfact; + gsl_sf_result lng1; + gsl_sf_result lng2; + double s1, s2; + const int stat_f = gsl_sf_lnfact_e(n, &lnfact); + const int stat_g1 = gsl_sf_lngamma_sgn_e(b + n, &lng1, &s1); + const int stat_g2 = gsl_sf_lngamma_sgn_e(b, &lng2, &s2); + const double lnpre_val = lnfact.val - (lng1.val - lng2.val); + const double lnpre_err = lnfact.err + lng1.err + lng2.err + + 2.0 * GSL_DBL_EPSILON * fabs(lnpre_val); + const int stat_e = gsl_sf_exp_mult_err_e(lnpre_val, lnpre_err, + s1*s2*lag.val, lag.err, + result); + return GSL_ERROR_SELECT_5(stat_e, stat_l, stat_g1, stat_g2, stat_f); + } + else { + gsl_sf_result lnbeta; + gsl_sf_lnbeta_e(b, n, &lnbeta); + if(fabs(lnbeta.val) < 0.1) { + /* As we have noted, when B(x,y) is near 1, + * evaluating log(B(x,y)) is not accurate. + * Instead we evaluate B(x,y) directly. + */ + const double ln_term_val = log(1.25*n); + const double ln_term_err = 2.0 * GSL_DBL_EPSILON * ln_term_val; + gsl_sf_result beta; + int stat_b = gsl_sf_beta_e(b, n, &beta); + int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err, + lag.val, lag.err, + result); + result->val *= beta.val/1.25; + result->err *= beta.val/1.25; + return GSL_ERROR_SELECT_3(stat_e, stat_l, stat_b); + } + else { + /* B(x,y) was not near 1, so it is safe to use + * the logarithmic values. + */ + const double ln_n = log(n); + const double ln_term_val = lnbeta.val + ln_n; + const double ln_term_err = lnbeta.err + 2.0 * GSL_DBL_EPSILON * fabs(ln_n); + int stat_e = gsl_sf_exp_mult_err_e(ln_term_val, ln_term_err, + lag.val, lag.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_l); + } + } +} + + +/* Handle negative integer a case for x > 0 and + * generic b. + * + * Combine [Abramowitz+Stegun, 13.6.9 + 13.6.27] + * M(-n,b,x) = (-1)^n / (b)_n U(-n,b,x) = n! / (b)_n Laguerre^(b-1)_n(x) + */ +#if 0 +static +int +hyperg_1F1_a_negint_U(const int a, const double b, const double x, gsl_sf_result * result) +{ + const int n = -a; + const double sgn = ( GSL_IS_ODD(n) ? -1.0 : 1.0 ); + double sgpoch; + gsl_sf_result lnpoch; + gsl_sf_result U; + const int stat_p = gsl_sf_lnpoch_sgn_e(b, n, &lnpoch, &sgpoch); + const int stat_U = gsl_sf_hyperg_U_e(-n, b, x, &U); + const int stat_e = gsl_sf_exp_mult_err_e(-lnpoch.val, lnpoch.err, + sgn * sgpoch * U.val, U.err, + result); + return GSL_ERROR_SELECT_3(stat_e, stat_U, stat_p); +} +#endif + + +/* Assumes a <= -1, b <= -1, and b <= a. + */ +static +int +hyperg_1F1_ab_negint(const int a, const int b, const double x, gsl_sf_result * result) +{ + if(x == 0.0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(x > 0.0) { + return hyperg_1F1_a_negint_poly(a, b, x, result); + } + else { + /* Apply a Kummer transformation to make x > 0 so + * we can evaluate the polynomial safely. Of course, + * this assumes b <= a, which must be true for + * a<0 and b<0, since otherwise the thing is undefined. + */ + gsl_sf_result K; + int stat_K = hyperg_1F1_a_negint_poly(b-a, b, -x, &K); + int stat_e = gsl_sf_exp_mult_err_e(x, 2.0 * GSL_DBL_EPSILON * fabs(x), + K.val, K.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } +} + + +/* [Abramowitz+Stegun, 13.1.3] + * + * M(a,b,x) = Gamma(1+a-b)/Gamma(2-b) x^(1-b) * + * { Gamma(b)/Gamma(a) M(1+a-b,2-b,x) - (b-1) U(1+a-b,2-b,x) } + * + * b not an integer >= 2 + * a-b not a negative integer + */ +static +int +hyperg_1F1_U(const double a, const double b, const double x, gsl_sf_result * result) +{ + const double bp = 2.0 - b; + const double ap = a - b + 1.0; + + gsl_sf_result lg_ap, lg_bp; + double sg_ap; + int stat_lg0 = gsl_sf_lngamma_sgn_e(ap, &lg_ap, &sg_ap); + int stat_lg1 = gsl_sf_lngamma_e(bp, &lg_bp); + int stat_lg2 = GSL_ERROR_SELECT_2(stat_lg0, stat_lg1); + double t1 = (bp-1.0) * log(x); + double lnpre_val = lg_ap.val - lg_bp.val + t1; + double lnpre_err = lg_ap.err + lg_bp.err + 2.0 * GSL_DBL_EPSILON * fabs(t1); + + gsl_sf_result lg_2mbp, lg_1papmbp; + double sg_2mbp, sg_1papmbp; + int stat_lg3 = gsl_sf_lngamma_sgn_e(2.0-bp, &lg_2mbp, &sg_2mbp); + int stat_lg4 = gsl_sf_lngamma_sgn_e(1.0+ap-bp, &lg_1papmbp, &sg_1papmbp); + int stat_lg5 = GSL_ERROR_SELECT_2(stat_lg3, stat_lg4); + double lnc1_val = lg_2mbp.val - lg_1papmbp.val; + double lnc1_err = lg_2mbp.err + lg_1papmbp.err + + GSL_DBL_EPSILON * (fabs(lg_2mbp.val) + fabs(lg_1papmbp.val)); + + gsl_sf_result M; + gsl_sf_result_e10 U; + int stat_F = gsl_sf_hyperg_1F1_e(ap, bp, x, &M); + int stat_U = gsl_sf_hyperg_U_e10_e(ap, bp, x, &U); + int stat_FU = GSL_ERROR_SELECT_2(stat_F, stat_U); + + gsl_sf_result_e10 term_M; + int stat_e0 = gsl_sf_exp_mult_err_e10_e(lnc1_val, lnc1_err, + sg_2mbp*sg_1papmbp*M.val, M.err, + &term_M); + + const double ombp = 1.0 - bp; + const double Uee_val = U.e10*M_LN10; + const double Uee_err = 2.0 * GSL_DBL_EPSILON * fabs(Uee_val); + const double Mee_val = term_M.e10*M_LN10; + const double Mee_err = 2.0 * GSL_DBL_EPSILON * fabs(Mee_val); + int stat_e1; + + /* Do a little dance with the exponential prefactors + * to avoid overflows in intermediate results. + */ + if(Uee_val > Mee_val) { + const double factorM_val = exp(Mee_val-Uee_val); + const double factorM_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorM_val; + const double inner_val = term_M.val*factorM_val - ombp*U.val; + const double inner_err = + term_M.err*factorM_val + fabs(ombp) * U.err + + fabs(term_M.val) * factorM_err + + GSL_DBL_EPSILON * (fabs(term_M.val*factorM_val) + fabs(ombp*U.val)); + stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Uee_val, lnpre_err+Uee_err, + sg_ap*inner_val, inner_err, + result); + } + else { + const double factorU_val = exp(Uee_val - Mee_val); + const double factorU_err = 2.0 * GSL_DBL_EPSILON * (fabs(Mee_val-Uee_val)+1.0) * factorU_val; + const double inner_val = term_M.val - ombp*factorU_val*U.val; + const double inner_err = + term_M.err + fabs(ombp*factorU_val*U.err) + + fabs(ombp*factorU_err*U.val) + + GSL_DBL_EPSILON * (fabs(term_M.val) + fabs(ombp*factorU_val*U.val)); + stat_e1 = gsl_sf_exp_mult_err_e(lnpre_val+Mee_val, lnpre_err+Mee_err, + sg_ap*inner_val, inner_err, + result); + } + + return GSL_ERROR_SELECT_5(stat_e1, stat_e0, stat_FU, stat_lg5, stat_lg2); +} + + +/* Handle case of generic positive a, b. + * Assumes b-a is not a negative integer. + */ +static +int +hyperg_1F1_ab_pos(const double a, const double b, + const double x, + gsl_sf_result * result) +{ + const double ax = fabs(x); + + if( ( b < 10.0 && a < 10.0 && ax < 5.0 ) + || ( b > a*ax ) + || ( b > a && ax < 5.0 ) + ) { + return gsl_sf_hyperg_1F1_series_e(a, b, x, result); + } + else if( x < -100.0 + && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.7*fabs(x) + ) { + /* Large negative x asymptotic. + */ + return hyperg_1F1_asymp_negx(a, b, x, result); + } + else if( x > 100.0 + && GSL_MAX_DBL(fabs(b-a),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.7*fabs(x) + ) { + /* Large positive x asymptotic. + */ + return hyperg_1F1_asymp_posx(a, b, x, result); + } + else if(fabs(b-a) <= 1.0) { + /* Directly handle b near a. + */ + return hyperg_1F1_beps_bgt0(a-b, b, x, result); /* a = b + eps */ + } + + else if(b > a && b >= 2*a + x) { + /* Use the Gautschi CF series, then + * recurse backward to a near 0 for normalization. + * This will work for either sign of x. + */ + double rap; + int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap); + double ra = 1.0 + x/a * rap; + + double Ma = GSL_SQRT_DBL_MIN; + double Map1 = ra * Ma; + double Mnp1 = Map1; + double Mn = Ma; + double Mnm1; + gsl_sf_result Mn_true; + int stat_Mt; + double n; + for(n=a; n>0.5; n -= 1.0) { + Mnm1 = (n * Mnp1 - (2.0*n-b+x) * Mn) / (b-n); + Mnp1 = Mn; + Mn = Mnm1; + } + + stat_Mt = hyperg_1F1_small_a_bgt0(n, b, x, &Mn_true); + + result->val = (Ma/Mn) * Mn_true.val; + result->err = fabs(Ma/Mn) * Mn_true.err; + result->err += 2.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(result->val); + return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1); + } + else if(b > a && b < 2*a + x && b > x) { + /* Use the Gautschi series representation of + * the continued fraction. Then recurse forward + * to near the a=b line for normalization. This will + * work for either sign of x, although we do need + * to check for b > x, which is relevant when x is positive. + */ + gsl_sf_result Mn_true; + int stat_Mt; + double rap; + int stat_CF1 = hyperg_1F1_CF1_p_ser(a, b, x, &rap); + double ra = 1.0 + x/a * rap; + double Ma = GSL_SQRT_DBL_MIN; + double Mnm1 = Ma; + double Mn = ra * Mnm1; + double Mnp1; + double n; + for(n=a+1.0; n<b-0.5; n += 1.0) { + Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n; + Mnm1 = Mn; + Mn = Mnp1; + } + stat_Mt = hyperg_1F1_beps_bgt0(n-b, b, x, &Mn_true); + result->val = Ma/Mn * Mn_true.val; + result->err = fabs(Ma/Mn) * Mn_true.err; + result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-a)+1.0) * fabs(result->val); + return GSL_ERROR_SELECT_2(stat_Mt, stat_CF1); + } + else if(x >= 0.0) { + + if(b < a) { + /* Forward recursion on a from a=b+eps-1,b+eps. + */ + double N = floor(a-b); + double eps = a - b - N; + gsl_sf_result r_M0; + gsl_sf_result r_M1; + int stat_0 = hyperg_1F1_beps_bgt0(eps-1.0, b, x, &r_M0); + int stat_1 = hyperg_1F1_beps_bgt0(eps, b, x, &r_M1); + double M0 = r_M0.val; + double M1 = r_M1.val; + + double Mam1 = M0; + double Ma = M1; + double Map1; + double ap; + double start_pair = fabs(M0) + fabs(M1); + double minim_pair = GSL_DBL_MAX; + double pair_ratio; + double rat_0 = fabs(r_M0.err/r_M0.val); + double rat_1 = fabs(r_M1.err/r_M1.val); + for(ap=b+eps; ap<a-0.1; ap += 1.0) { + Map1 = ((b-ap)*Mam1 + (2.0*ap-b+x)*Ma)/ap; + Mam1 = Ma; + Ma = Map1; + minim_pair = GSL_MIN_DBL(fabs(Mam1) + fabs(Ma), minim_pair); + } + pair_ratio = start_pair/minim_pair; + result->val = Ma; + result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Ma); + result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Ma); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(Ma); + return GSL_ERROR_SELECT_2(stat_0, stat_1); + } + else { + /* b > a + * b < 2a + x + * b <= x + * + * Recurse forward on a from a=eps,eps+1. + */ + double eps = a - floor(a); + gsl_sf_result r_Mnm1; + gsl_sf_result r_Mn; + int stat_0 = hyperg_1F1_small_a_bgt0(eps, b, x, &r_Mnm1); + int stat_1 = hyperg_1F1_small_a_bgt0(eps+1.0, b, x, &r_Mn); + double Mnm1 = r_Mnm1.val; + double Mn = r_Mn.val; + double Mnp1; + + double n; + double start_pair = fabs(Mn) + fabs(Mnm1); + double minim_pair = GSL_DBL_MAX; + double pair_ratio; + double rat_0 = fabs(r_Mnm1.err/r_Mnm1.val); + double rat_1 = fabs(r_Mn.err/r_Mn.val); + for(n=eps+1.0; n<a-0.1; n++) { + Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n; + Mnm1 = Mn; + Mn = Mnp1; + minim_pair = GSL_MIN_DBL(fabs(Mn) + fabs(Mnm1), minim_pair); + } + pair_ratio = start_pair/minim_pair; + result->val = Mn; + result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(a)+1.0) * fabs(Mn); + result->err += 2.0 * (rat_0 + rat_1) * pair_ratio*pair_ratio * fabs(Mn); + result->err += 2.0 * GSL_DBL_EPSILON * fabs(Mn); + return GSL_ERROR_SELECT_2(stat_0, stat_1); + } + } + else { + /* x < 0 + * b < a + */ + + if(a <= 0.5*(b-x) || a >= -x) { + /* Recurse down in b, from near the a=b line, b=a+eps,a+eps-1. + */ + double N = floor(a - b); + double eps = 1.0 + N - a + b; + gsl_sf_result r_Manp1; + gsl_sf_result r_Man; + int stat_0 = hyperg_1F1_beps_bgt0(-eps, a+eps, x, &r_Manp1); + int stat_1 = hyperg_1F1_beps_bgt0(1.0-eps, a+eps-1.0, x, &r_Man); + double Manp1 = r_Manp1.val; + double Man = r_Man.val; + double Manm1; + + double n; + double start_pair = fabs(Manp1) + fabs(Man); + double minim_pair = GSL_DBL_MAX; + double pair_ratio; + double rat_0 = fabs(r_Manp1.err/r_Manp1.val); + double rat_1 = fabs(r_Man.err/r_Man.val); + for(n=a+eps-1.0; n>b+0.1; n -= 1.0) { + Manm1 = (-n*(1-n-x)*Man - x*(n-a)*Manp1)/(n*(n-1.0)); + Manp1 = Man; + Man = Manm1; + minim_pair = GSL_MIN_DBL(fabs(Manp1) + fabs(Man), minim_pair); + } + + /* FIXME: this is a nasty little hack; there is some + (transient?) instability in this recurrence for some + values. I can tell when it happens, which is when + this pair_ratio is large. But I do not know how to + measure the error in terms of it. I guessed quadratic + below, but it is probably worse than that. + */ + pair_ratio = start_pair/minim_pair; + result->val = Man; + result->err = 2.0 * (rat_0 + rat_1 + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Man); + result->err *= pair_ratio*pair_ratio + 1.0; + return GSL_ERROR_SELECT_2(stat_0, stat_1); + } + else { + /* Pick a0 such that b ~= 2a0 + x, then + * recurse down in b from a0,a0 to determine + * the values near the line b=2a+x. Then recurse + * forward on a from a0. + */ + double epsa = a - floor(a); + double a0 = floor(0.5*(b-x)) + epsa; + double N = floor(a0 - b); + double epsb = 1.0 + N - a0 + b; + double Ma0b; + double Ma0bp1; + double Ma0p1b; + int stat_a0; + double Mnm1; + double Mn; + double Mnp1; + double n; + double err_rat; + { + gsl_sf_result r_Ma0np1; + gsl_sf_result r_Ma0n; + int stat_0 = hyperg_1F1_beps_bgt0(-epsb, a0+epsb, x, &r_Ma0np1); + int stat_1 = hyperg_1F1_beps_bgt0(1.0-epsb, a0+epsb-1.0, x, &r_Ma0n); + double Ma0np1 = r_Ma0np1.val; + double Ma0n = r_Ma0n.val; + double Ma0nm1; + + err_rat = fabs(r_Ma0np1.err/r_Ma0np1.val) + fabs(r_Ma0n.err/r_Ma0n.val); + + for(n=a0+epsb-1.0; n>b+0.1; n -= 1.0) { + Ma0nm1 = (-n*(1-n-x)*Ma0n - x*(n-a0)*Ma0np1)/(n*(n-1.0)); + Ma0np1 = Ma0n; + Ma0n = Ma0nm1; + } + Ma0bp1 = Ma0np1; + Ma0b = Ma0n; + Ma0p1b = (b*(a0+x)*Ma0b+x*(a0-b)*Ma0bp1)/(a0*b); /* right-down hook */ + stat_a0 = GSL_ERROR_SELECT_2(stat_0, stat_1); + } + + + /* Initialise the recurrence correctly BJG */ + + if (a0 >= a - 0.1) + { + Mn = Ma0b; + } + else if (a0 + 1>= a - 0.1) + { + Mn = Ma0p1b; + } + else + { + Mnm1 = Ma0b; + Mn = Ma0p1b; + + for(n=a0+1.0; n<a-0.1; n += 1.0) { + Mnp1 = ((b-n)*Mnm1 + (2*n-b+x)*Mn)/n; + Mnm1 = Mn; + Mn = Mnp1; + } + } + + result->val = Mn; + result->err = (err_rat + GSL_DBL_EPSILON) * (fabs(b-a)+1.0) * fabs(Mn); + return stat_a0; + } + } +} + + +/* Assumes b != integer + * Assumes a != integer when x > 0 + * Assumes b-a != neg integer when x < 0 + */ +static +int +hyperg_1F1_ab_neg(const double a, const double b, const double x, + gsl_sf_result * result) +{ + const double bma = b - a; + const double abs_x = fabs(x); + const double abs_a = fabs(a); + const double abs_b = fabs(b); + const double size_a = GSL_MAX(abs_a, 1.0); + const double size_b = GSL_MAX(abs_b, 1.0); + const int bma_integer = ( bma - floor(bma+0.5) < _1F1_INT_THRESHOLD ); + + if( (abs_a < 10.0 && abs_b < 10.0 && abs_x < 5.0) + || (b > 0.8*GSL_MAX(fabs(a),1.0)*fabs(x)) + ) { + return gsl_sf_hyperg_1F1_series_e(a, b, x, result); + } + else if( x > 0.0 + && size_b > size_a + && size_a*log(M_E*x/size_b) < GSL_LOG_DBL_EPSILON+7.0 + ) { + /* Series terms are positive definite up until + * there is a sign change. But by then the + * terms are small due to the last condition. + */ + return gsl_sf_hyperg_1F1_series_e(a, b, x, result); + } + else if( (abs_x < 5.0 && fabs(bma) < 10.0 && abs_b < 10.0) + || (b > 0.8*GSL_MAX_DBL(fabs(bma),1.0)*abs_x) + ) { + /* Use Kummer transformation to render series safe. + */ + gsl_sf_result Kummer_1F1; + int stat_K = gsl_sf_hyperg_1F1_series_e(bma, b, -x, &Kummer_1F1); + int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x), + Kummer_1F1.val, Kummer_1F1.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } + else if( x < -30.0 + && GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x) + ) { + /* Large negative x asymptotic. + * Note that we do not check if b-a is a negative integer. + */ + return hyperg_1F1_asymp_negx(a, b, x, result); + } + else if( x > 100.0 + && GSL_MAX_DBL(fabs(bma),1.0)*GSL_MAX_DBL(fabs(1.0-a),1.0) < 0.99*fabs(x) + ) { + /* Large positive x asymptotic. + * Note that we do not check if a is a negative integer. + */ + return hyperg_1F1_asymp_posx(a, b, x, result); + } + else if(x > 0.0 && !(bma_integer && bma > 0.0)) { + return hyperg_1F1_U(a, b, x, result); + } + else { + /* FIXME: if all else fails, try the series... BJG */ + if (x < 0.0) { + /* Apply Kummer Transformation */ + int status = gsl_sf_hyperg_1F1_series_e(b-a, b, -x, result); + double K_factor = exp(x); + result->val *= K_factor; + result->err *= K_factor; + return status; + } else { + int status = gsl_sf_hyperg_1F1_series_e(a, b, x, result); + return status; + } + + /* Sadness... */ + /* result->val = 0.0; */ + /* result->err = 0.0; */ + /* GSL_ERROR ("error", GSL_EUNIMPL); */ + } +} + + +/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ + +int +gsl_sf_hyperg_1F1_int_e(const int a, const int b, const double x, gsl_sf_result * result) +{ + /* CHECK_POINTER(result) */ + + if(x == 0.0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(a == b) { + return gsl_sf_exp_e(x, result); + } + else if(b == 0) { + DOMAIN_ERROR(result); + } + else if(a == 0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(b < 0 && (a < b || a > 0)) { + /* Standard domain error due to singularity. */ + DOMAIN_ERROR(result); + } + else if(x > 100.0 && GSL_MAX_DBL(1.0,fabs(b-a))*GSL_MAX_DBL(1.0,fabs(1-a)) < 0.5 * x) { + /* x -> +Inf asymptotic */ + return hyperg_1F1_asymp_posx(a, b, x, result); + } + else if(x < -100.0 && GSL_MAX_DBL(1.0,fabs(a))*GSL_MAX_DBL(1.0,fabs(1+a-b)) < 0.5 * fabs(x)) { + /* x -> -Inf asymptotic */ + return hyperg_1F1_asymp_negx(a, b, x, result); + } + else if(a < 0 && b < 0) { + return hyperg_1F1_ab_negint(a, b, x, result); + } + else if(a < 0 && b > 0) { + /* Use Kummer to reduce it to the positive integer case. + * Note that b > a, strictly, since we already trapped b = a. + */ + gsl_sf_result Kummer_1F1; + int stat_K = hyperg_1F1_ab_posint(b-a, b, -x, &Kummer_1F1); + int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x), + Kummer_1F1.val, Kummer_1F1.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } + else { + /* a > 0 and b > 0 */ + return hyperg_1F1_ab_posint(a, b, x, result); + } +} + + +int +gsl_sf_hyperg_1F1_e(const double a, const double b, const double x, + gsl_sf_result * result + ) +{ + const double bma = b - a; + const double rinta = floor(a + 0.5); + const double rintb = floor(b + 0.5); + const double rintbma = floor(bma + 0.5); + const int a_integer = ( fabs(a-rinta) < _1F1_INT_THRESHOLD && rinta > INT_MIN && rinta < INT_MAX ); + const int b_integer = ( fabs(b-rintb) < _1F1_INT_THRESHOLD && rintb > INT_MIN && rintb < INT_MAX ); + const int bma_integer = ( fabs(bma-rintbma) < _1F1_INT_THRESHOLD && rintbma > INT_MIN && rintbma < INT_MAX ); + const int b_neg_integer = ( b < -0.1 && b_integer ); + const int a_neg_integer = ( a < -0.1 && a_integer ); + const int bma_neg_integer = ( bma < -0.1 && bma_integer ); + + /* CHECK_POINTER(result) */ + + if(x == 0.0) { + /* Testing for this before testing a and b + * is somewhat arbitrary. The result is that + * we have 1F1(a,0,0) = 1. + */ + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(b == 0.0) { + DOMAIN_ERROR(result); + } + else if(a == 0.0) { + result->val = 1.0; + result->err = 0.0; + return GSL_SUCCESS; + } + else if(a == b) { + /* case: a==b; exp(x) + * It's good to test exact equality now. + * We also test approximate equality later. + */ + return gsl_sf_exp_e(x, result); + } else if(fabs(b) < _1F1_INT_THRESHOLD && fabs(a) < _1F1_INT_THRESHOLD) { + /* a and b near zero: 1 + a/b (exp(x)-1) + */ + + /* Note that neither a nor b is zero, since + * we eliminated that with the above tests. + */ + + gsl_sf_result exm1; + int stat_e = gsl_sf_expm1_e(x, &exm1); + double sa = ( a > 0.0 ? 1.0 : -1.0 ); + double sb = ( b > 0.0 ? 1.0 : -1.0 ); + double lnab = log(fabs(a/b)); /* safe */ + gsl_sf_result hx; + int stat_hx = gsl_sf_exp_mult_err_e(lnab, GSL_DBL_EPSILON * fabs(lnab), + sa * sb * exm1.val, exm1.err, + &hx); + result->val = (hx.val == GSL_DBL_MAX ? hx.val : 1.0 + hx.val); /* FIXME: excessive paranoia ? what is DBL_MAX+1 ?*/ + result->err = hx.err; + return GSL_ERROR_SELECT_2(stat_hx, stat_e); + } else if (fabs(b) < _1F1_INT_THRESHOLD && fabs(x*a) < 1) { + /* b near zero and a not near zero + */ + const double m_arg = 1.0/(0.5*b); + gsl_sf_result F_renorm; + int stat_F = hyperg_1F1_renorm_b0(a, x, &F_renorm); + int stat_m = gsl_sf_multiply_err_e(m_arg, 2.0 * GSL_DBL_EPSILON * m_arg, + 0.5*F_renorm.val, 0.5*F_renorm.err, + result); + return GSL_ERROR_SELECT_2(stat_m, stat_F); + } + else if(a_integer && b_integer) { + /* Check for reduction to the integer case. + * Relies on the arbitrary "near an integer" test. + */ + return gsl_sf_hyperg_1F1_int_e((int)rinta, (int)rintb, x, result); + } + else if(b_neg_integer && !(a_neg_integer && a > b)) { + /* Standard domain error due to + * uncancelled singularity. + */ + DOMAIN_ERROR(result); + } + else if(a_neg_integer) { + return hyperg_1F1_a_negint_lag((int)rinta, b, x, result); + } + else if(b > 0.0) { + if(-1.0 <= a && a <= 1.0) { + /* Handle small a explicitly. + */ + return hyperg_1F1_small_a_bgt0(a, b, x, result); + } + else if(bma_neg_integer) { + /* Catch this now, to avoid problems in the + * generic evaluation code. + */ + gsl_sf_result Kummer_1F1; + int stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &Kummer_1F1); + int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x), + Kummer_1F1.val, Kummer_1F1.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } + else if(a < 0.0 && fabs(x) < 100.0) { + /* Use Kummer to reduce it to the generic positive case. + * Note that b > a, strictly, since we already trapped b = a. + * Also b-(b-a)=a, and a is not a negative integer here, + * so the generic evaluation is safe. + */ + gsl_sf_result Kummer_1F1; + int stat_K = hyperg_1F1_ab_pos(b-a, b, -x, &Kummer_1F1); + int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x), + Kummer_1F1.val, Kummer_1F1.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } + else if (a > 0) { + /* a > 0.0 */ + return hyperg_1F1_ab_pos(a, b, x, result); + } else { + return gsl_sf_hyperg_1F1_series_e(a, b, x, result); + } + } + else { + /* b < 0.0 */ + + if(bma_neg_integer && x < 0.0) { + /* Handle this now to prevent problems + * in the generic evaluation. + */ + gsl_sf_result K; + int stat_K; + int stat_e; + if(a < 0.0) { + /* Kummer transformed version of safe polynomial. + * The condition a < 0 is equivalent to b < b-a, + * which is the condition required for the series + * to be positive definite here. + */ + stat_K = hyperg_1F1_a_negint_poly((int)rintbma, b, -x, &K); + } + else { + /* Generic eval for negative integer a. */ + stat_K = hyperg_1F1_a_negint_lag((int)rintbma, b, -x, &K); + } + stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x), + K.val, K.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } + else if(a > 0.0) { + /* Use Kummer to reduce it to the generic negative case. + */ + gsl_sf_result K; + int stat_K = hyperg_1F1_ab_neg(b-a, b, -x, &K); + int stat_e = gsl_sf_exp_mult_err_e(x, GSL_DBL_EPSILON * fabs(x), + K.val, K.err, + result); + return GSL_ERROR_SELECT_2(stat_e, stat_K); + } + else { + return hyperg_1F1_ab_neg(a, b, x, result); + } + } +} + + + +#if 0 + /* Luke in the canonical case. + */ + if(x < 0.0 && !a_neg_integer && !bma_neg_integer) { + double prec; + return hyperg_1F1_luke(a, b, x, result, &prec); + } + + + /* Luke with Kummer transformation. + */ + if(x > 0.0 && !a_neg_integer && !bma_neg_integer) { + double prec; + double Kummer_1F1; + double ex; + int stat_F = hyperg_1F1_luke(b-a, b, -x, &Kummer_1F1, &prec); + int stat_e = gsl_sf_exp_e(x, &ex); + if(stat_F == GSL_SUCCESS && stat_e == GSL_SUCCESS) { + double lnr = log(fabs(Kummer_1F1)) + x; + if(lnr < GSL_LOG_DBL_MAX) { + *result = ex * Kummer_1F1; + return GSL_SUCCESS; + } + else { + *result = GSL_POSINF; + GSL_ERROR ("overflow", GSL_EOVRFLW); + } + } + else if(stat_F != GSL_SUCCESS) { + *result = 0.0; + return stat_F; + } + else { + *result = 0.0; + return stat_e; + } + } +#endif + + + +/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ + +#include "eval.h" + +double gsl_sf_hyperg_1F1_int(const int m, const int n, double x) +{ + EVAL_RESULT(gsl_sf_hyperg_1F1_int_e(m, n, x, &result)); +} + +double gsl_sf_hyperg_1F1(double a, double b, double x) +{ + EVAL_RESULT(gsl_sf_hyperg_1F1_e(a, b, x, &result)); +} |