/* specfunc/hyperg_1F1.c * * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman * Copyright (C) 2010 Brian Gough * * 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 3 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 #include #include #include #include #include #include #include #include #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 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; nval = 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; nval = 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; nval = 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 = (int) 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; nval = 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; nval = 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; apval = 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; nval = 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; nval = 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) < 2*GSL_LOG_DBL_MAX) { /* 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)); }