| /* specfunc/hyperg_U.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_exp.h> |
| #include <gsl/gsl_sf_gamma.h> |
| #include <gsl/gsl_sf_bessel.h> |
| #include <gsl/gsl_sf_laguerre.h> |
| #include <gsl/gsl_sf_pow_int.h> |
| #include <gsl/gsl_sf_hyperg.h> |
| |
| #include "error.h" |
| #include "hyperg.h" |
| |
| #define INT_THRESHOLD (1000.0*GSL_DBL_EPSILON) |
| |
| #define SERIES_EVAL_OK(a,b,x) ((fabs(a) < 5 && b < 5 && x < 2.0) || (fabs(a) < 10 && b < 10 && x < 1.0)) |
| |
| #define ASYMP_EVAL_OK(a,b,x) (GSL_MAX_DBL(fabs(a),1.0)*GSL_MAX_DBL(fabs(1.0+a-b),1.0) < 0.99*fabs(x)) |
| |
| |
| /* Log[U(a,2a,x)] |
| * [Abramowitz+stegun, 13.6.21] |
| * Assumes x > 0, a > 1/2. |
| */ |
| static |
| int |
| hyperg_lnU_beq2a(const double a, const double x, gsl_sf_result * result) |
| { |
| const double lx = log(x); |
| const double nu = a - 0.5; |
| const double lnpre = 0.5*(x - M_LNPI) - nu*lx; |
| gsl_sf_result lnK; |
| gsl_sf_bessel_lnKnu_e(nu, 0.5*x, &lnK); |
| result->val = lnpre + lnK.val; |
| result->err = 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + 0.5*M_LNPI + fabs(nu*lx)); |
| result->err += lnK.err; |
| result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
| return GSL_SUCCESS; |
| } |
| |
| |
| /* Evaluate u_{N+1}/u_N by Steed's continued fraction method. |
| * |
| * u_N := Gamma[a+N]/Gamma[a] U(a + N, b, x) |
| * |
| * u_{N+1}/u_N = (a+N) U(a+N+1,b,x)/U(a+N,b,x) |
| */ |
| static |
| int |
| hyperg_U_CF1(const double a, const double b, const int N, const double x, |
| double * result, int * count) |
| { |
| const double RECUR_BIG = GSL_SQRT_DBL_MAX; |
| const int maxiter = 20000; |
| int n = 1; |
| double Anm2 = 1.0; |
| double Bnm2 = 0.0; |
| double Anm1 = 0.0; |
| double Bnm1 = 1.0; |
| double a1 = -(a + N); |
| double b1 = (b - 2.0*a - x - 2.0*(N+1)); |
| 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 + n - b)*(a + N + n - 1.0); |
| bn = (b - 2.0*a - x - 2.0*(N+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; |
| *count = n; |
| |
| if(n == maxiter) |
| GSL_ERROR ("error", GSL_EMAXITER); |
| else |
| return GSL_SUCCESS; |
| } |
| |
| |
| /* Large x asymptotic for x^a U(a,b,x) |
| * Based on SLATEC D9CHU() [W. Fullerton] |
| * |
| * Uses a rational approximation due to Luke. |
| * See [Luke, Algorithms for the Computation of Special Functions, p. 252] |
| * [Luke, Utilitas Math. (1977)] |
| * |
| * z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z) |
| * |
| * This assumes that a is not a negative integer and |
| * that 1+a-b is not a negative integer. If one of them |
| * is, then the 2F0 actually terminates, the above |
| * relation is an equality, and the sum should be |
| * evaluated directly [see below]. |
| */ |
| static |
| int |
| d9chu(const double a, const double b, const double x, gsl_sf_result * result) |
| { |
| const double EPS = 8.0 * GSL_DBL_EPSILON; /* EPS = 4.0D0*D1MACH(4) */ |
| const int maxiter = 500; |
| double aa[4], bb[4]; |
| int i; |
| |
| double bp = 1.0 + a - b; |
| double ab = a*bp; |
| double ct2 = 2.0 * (x - ab); |
| double sab = a + bp; |
| |
| double ct3 = sab + 1.0 + ab; |
| double anbn = ct3 + sab + 3.0; |
| double ct1 = 1.0 + 2.0*x/anbn; |
| |
| bb[0] = 1.0; |
| aa[0] = 1.0; |
| |
| bb[1] = 1.0 + 2.0*x/ct3; |
| aa[1] = 1.0 + ct2/ct3; |
| |
| bb[2] = 1.0 + 6.0*ct1*x/ct3; |
| aa[2] = 1.0 + 6.0*ab/anbn + 3.0*ct1*ct2/ct3; |
| |
| for(i=4; i<maxiter; i++) { |
| int j; |
| double c2; |
| double d1z; |
| double g1, g2, g3; |
| double x2i1 = 2*i - 3; |
| ct1 = x2i1/(x2i1-2.0); |
| anbn += x2i1 + sab; |
| ct2 = (x2i1 - 1.0)/anbn; |
| c2 = x2i1*ct2 - 1.0; |
| d1z = 2.0*x2i1*x/anbn; |
| |
| ct3 = sab*ct2; |
| g1 = d1z + ct1*(c2+ct3); |
| g2 = d1z - c2; |
| g3 = ct1*(1.0 - ct3 - 2.0*ct2); |
| |
| bb[3] = g1*bb[2] + g2*bb[1] + g3*bb[0]; |
| aa[3] = g1*aa[2] + g2*aa[1] + g3*aa[0]; |
| |
| if(fabs(aa[3]*bb[0]-aa[0]*bb[3]) < EPS*fabs(bb[3]*bb[0])) break; |
| |
| for(j=0; j<3; j++) { |
| aa[j] = aa[j+1]; |
| bb[j] = bb[j+1]; |
| } |
| } |
| |
| result->val = aa[3]/bb[3]; |
| result->err = 8.0 * GSL_DBL_EPSILON * fabs(result->val); |
| |
| if(i == maxiter) { |
| GSL_ERROR ("error", GSL_EMAXITER); |
| } |
| else { |
| return GSL_SUCCESS; |
| } |
| } |
| |
| |
| /* Evaluate asymptotic for z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z) |
| * We check for termination of the 2F0 as a special case. |
| * Assumes x > 0. |
| * Also assumes a,b are not too large compared to x. |
| */ |
| static |
| int |
| hyperg_zaU_asymp(const double a, const double b, const double x, gsl_sf_result *result) |
| { |
| const double ap = a; |
| const double bp = 1.0 + a - b; |
| const double rintap = floor(ap + 0.5); |
| const double rintbp = floor(bp + 0.5); |
| const int ap_neg_int = ( ap < 0.0 && fabs(ap - rintap) < INT_THRESHOLD ); |
| const int bp_neg_int = ( bp < 0.0 && fabs(bp - rintbp) < INT_THRESHOLD ); |
| |
| if(ap_neg_int || bp_neg_int) { |
| /* Evaluate 2F0 polynomial. |
| */ |
| double mxi = -1.0/x; |
| double nmax = -(int)(GSL_MIN(ap,bp) - 0.1); |
| double tn = 1.0; |
| double sum = 1.0; |
| double n = 1.0; |
| double sum_err = 0.0; |
| while(n <= nmax) { |
| double apn = (ap+n-1.0); |
| double bpn = (bp+n-1.0); |
| tn *= ((apn/n)*mxi)*bpn; |
| sum += tn; |
| sum_err += 2.0 * GSL_DBL_EPSILON * fabs(tn); |
| n += 1.0; |
| } |
| result->val = sum; |
| result->err = sum_err; |
| result->err += 2.0 * GSL_DBL_EPSILON * (fabs(nmax)+1.0) * fabs(sum); |
| return GSL_SUCCESS; |
| } |
| else { |
| return d9chu(a,b,x,result); |
| } |
| } |
| |
| |
| /* Evaluate finite sum which appears below. |
| */ |
| static |
| int |
| hyperg_U_finite_sum(int N, double a, double b, double x, double xeps, |
| gsl_sf_result * result) |
| { |
| int i; |
| double sum_val; |
| double sum_err; |
| |
| if(N <= 0) { |
| double t_val = 1.0; |
| double t_err = 0.0; |
| gsl_sf_result poch; |
| int stat_poch; |
| |
| sum_val = 1.0; |
| sum_err = 0.0; |
| for(i=1; i<= -N; i++) { |
| const double xi1 = i - 1; |
| const double mult = (a+xi1)*x/((b+xi1)*(xi1+1.0)); |
| t_val *= mult; |
| t_err += fabs(mult) * t_err + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON; |
| sum_val += t_val; |
| sum_err += t_err; |
| } |
| |
| stat_poch = gsl_sf_poch_e(1.0+a-b, -a, &poch); |
| |
| result->val = sum_val * poch.val; |
| result->err = fabs(sum_val) * poch.err + sum_err * fabs(poch.val); |
| result->err += fabs(poch.val) * (fabs(N) + 2.0) * GSL_DBL_EPSILON * fabs(sum_val); |
| result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
| result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */ |
| return stat_poch; |
| } |
| else { |
| const int M = N - 2; |
| if(M < 0) { |
| result->val = 0.0; |
| result->err = 0.0; |
| return GSL_SUCCESS; |
| } |
| else { |
| gsl_sf_result gbm1; |
| gsl_sf_result gamr; |
| int stat_gbm1; |
| int stat_gamr; |
| double t_val = 1.0; |
| double t_err = 0.0; |
| |
| sum_val = 1.0; |
| sum_err = 0.0; |
| for(i=1; i<=M; i++) { |
| const double mult = (a-b+i)*x/((1.0-b+i)*i); |
| t_val *= mult; |
| t_err += t_err * fabs(mult) + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON; |
| sum_val += t_val; |
| sum_err += t_err; |
| } |
| |
| stat_gbm1 = gsl_sf_gamma_e(b-1.0, &gbm1); |
| stat_gamr = gsl_sf_gammainv_e(a, &gamr); |
| |
| if(stat_gbm1 == GSL_SUCCESS) { |
| gsl_sf_result powx1N; |
| int stat_p = gsl_sf_pow_int_e(x, 1-N, &powx1N); |
| double pe_val = powx1N.val * xeps; |
| double pe_err = powx1N.err * fabs(xeps) + 2.0 * GSL_DBL_EPSILON * fabs(pe_val); |
| double coeff_val = gbm1.val * gamr.val * pe_val; |
| double coeff_err = gbm1.err * fabs(gamr.val * pe_val) |
| + gamr.err * fabs(gbm1.val * pe_val) |
| + fabs(gbm1.val * gamr.val) * pe_err |
| + 2.0 * GSL_DBL_EPSILON * fabs(coeff_val); |
| |
| result->val = sum_val * coeff_val; |
| result->err = fabs(sum_val) * coeff_err + sum_err * fabs(coeff_val); |
| result->err += 2.0 * GSL_DBL_EPSILON * (M+2.0) * fabs(result->val); |
| result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */ |
| return stat_p; |
| } |
| else { |
| result->val = 0.0; |
| result->err = 0.0; |
| return stat_gbm1; |
| } |
| } |
| } |
| } |
| |
| |
| /* Based on SLATEC DCHU() [W. Fullerton] |
| * Assumes x > 0. |
| * This is just a series summation method, and |
| * it is not good for large a. |
| * |
| * I patched up the window for 1+a-b near zero. [GJ] |
| */ |
| static |
| int |
| hyperg_U_series(const double a, const double b, const double x, gsl_sf_result * result) |
| { |
| const double EPS = 2.0 * GSL_DBL_EPSILON; /* EPS = D1MACH(3) */ |
| const double SQRT_EPS = M_SQRT2 * GSL_SQRT_DBL_EPSILON; |
| |
| if(fabs(1.0 + a - b) < SQRT_EPS) { |
| /* Original Comment: ALGORITHM IS BAD WHEN 1+A-B IS NEAR ZERO FOR SMALL X |
| */ |
| /* We can however do the following: |
| * U(a,b,x) = U(a,a+1,x) when 1+a-b=0 |
| * and U(a,a+1,x) = x^(-a). |
| */ |
| double lnr = -a * log(x); |
| int stat_e = gsl_sf_exp_e(lnr, result); |
| result->err += 2.0 * SQRT_EPS * fabs(result->val); |
| return stat_e; |
| } |
| else { |
| double aintb = ( b < 0.0 ? ceil(b-0.5) : floor(b+0.5) ); |
| double beps = b - aintb; |
| int N = aintb; |
| |
| double lnx = log(x); |
| double xeps = exp(-beps*lnx); |
| |
| /* Evaluate finite sum. |
| */ |
| gsl_sf_result sum; |
| int stat_sum = hyperg_U_finite_sum(N, a, b, x, xeps, &sum); |
| |
| |
| /* Evaluate infinite sum. */ |
| |
| int istrt = ( N < 1 ? 1-N : 0 ); |
| double xi = istrt; |
| |
| gsl_sf_result gamr; |
| gsl_sf_result powx; |
| int stat_gamr = gsl_sf_gammainv_e(1.0+a-b, &gamr); |
| int stat_powx = gsl_sf_pow_int_e(x, istrt, &powx); |
| double sarg = beps*M_PI; |
| double sfact = ( sarg != 0.0 ? sarg/sin(sarg) : 1.0 ); |
| double factor_val = sfact * ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * gamr.val * powx.val; |
| double factor_err = fabs(gamr.val) * powx.err + fabs(powx.val) * gamr.err |
| + 2.0 * GSL_DBL_EPSILON * fabs(factor_val); |
| |
| gsl_sf_result pochai; |
| gsl_sf_result gamri1; |
| gsl_sf_result gamrni; |
| int stat_pochai = gsl_sf_poch_e(a, xi, &pochai); |
| int stat_gamri1 = gsl_sf_gammainv_e(xi + 1.0, &gamri1); |
| int stat_gamrni = gsl_sf_gammainv_e(aintb + xi, &gamrni); |
| int stat_gam123 = GSL_ERROR_SELECT_3(stat_gamr, stat_gamri1, stat_gamrni); |
| int stat_gamall = GSL_ERROR_SELECT_4(stat_sum, stat_gam123, stat_pochai, stat_powx); |
| |
| gsl_sf_result pochaxibeps; |
| gsl_sf_result gamrxi1beps; |
| int stat_pochaxibeps = gsl_sf_poch_e(a, xi-beps, &pochaxibeps); |
| int stat_gamrxi1beps = gsl_sf_gammainv_e(xi + 1.0 - beps, &gamrxi1beps); |
| |
| int stat_all = GSL_ERROR_SELECT_3(stat_gamall, stat_pochaxibeps, stat_gamrxi1beps); |
| |
| double b0_val = factor_val * pochaxibeps.val * gamrni.val * gamrxi1beps.val; |
| double b0_err = fabs(factor_val * pochaxibeps.val * gamrni.val) * gamrxi1beps.err |
| + fabs(factor_val * pochaxibeps.val * gamrxi1beps.val) * gamrni.err |
| + fabs(factor_val * gamrni.val * gamrxi1beps.val) * pochaxibeps.err |
| + fabs(pochaxibeps.val * gamrni.val * gamrxi1beps.val) * factor_err |
| + 2.0 * GSL_DBL_EPSILON * fabs(b0_val); |
| |
| if(fabs(xeps-1.0) < 0.5) { |
| /* |
| C X**(-BEPS) IS CLOSE TO 1.0D0, SO WE MUST BE |
| C CAREFUL IN EVALUATING THE DIFFERENCES. |
| */ |
| int i; |
| gsl_sf_result pch1ai; |
| gsl_sf_result pch1i; |
| gsl_sf_result poch1bxibeps; |
| int stat_pch1ai = gsl_sf_pochrel_e(a + xi, -beps, &pch1ai); |
| int stat_pch1i = gsl_sf_pochrel_e(xi + 1.0 - beps, beps, &pch1i); |
| int stat_poch1bxibeps = gsl_sf_pochrel_e(b+xi, -beps, &poch1bxibeps); |
| double c0_t1_val = beps*pch1ai.val*pch1i.val; |
| double c0_t1_err = fabs(beps) * fabs(pch1ai.val) * pch1i.err |
| + fabs(beps) * fabs(pch1i.val) * pch1ai.err |
| + 2.0 * GSL_DBL_EPSILON * fabs(c0_t1_val); |
| double c0_t2_val = -poch1bxibeps.val + pch1ai.val - pch1i.val + c0_t1_val; |
| double c0_t2_err = poch1bxibeps.err + pch1ai.err + pch1i.err + c0_t1_err |
| + 2.0 * GSL_DBL_EPSILON * fabs(c0_t2_val); |
| double c0_val = factor_val * pochai.val * gamrni.val * gamri1.val * c0_t2_val; |
| double c0_err = fabs(factor_val * pochai.val * gamrni.val * gamri1.val) * c0_t2_err |
| + fabs(factor_val * pochai.val * gamrni.val * c0_t2_val) * gamri1.err |
| + fabs(factor_val * pochai.val * gamri1.val * c0_t2_val) * gamrni.err |
| + fabs(factor_val * gamrni.val * gamri1.val * c0_t2_val) * pochai.err |
| + fabs(pochai.val * gamrni.val * gamri1.val * c0_t2_val) * factor_err |
| + 2.0 * GSL_DBL_EPSILON * fabs(c0_val); |
| /* |
| C XEPS1 = (1.0 - X**(-BEPS))/BEPS = (X**(-BEPS) - 1.0)/(-BEPS) |
| */ |
| gsl_sf_result dexprl; |
| int stat_dexprl = gsl_sf_exprel_e(-beps*lnx, &dexprl); |
| double xeps1_val = lnx * dexprl.val; |
| double xeps1_err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(beps*lnx)) * fabs(dexprl.val) |
| + fabs(lnx) * dexprl.err |
| + 2.0 * GSL_DBL_EPSILON * fabs(xeps1_val); |
| double dchu_val = sum.val + c0_val + xeps1_val*b0_val; |
| double dchu_err = sum.err + c0_err |
| + fabs(xeps1_val)*b0_err + xeps1_err * fabs(b0_val) |
| + fabs(b0_val*lnx)*dexprl.err |
| + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(c0_val) + fabs(xeps1_val*b0_val)); |
| double xn = N; |
| double t_val; |
| double t_err; |
| |
| stat_all = GSL_ERROR_SELECT_5(stat_all, stat_dexprl, stat_poch1bxibeps, stat_pch1i, stat_pch1ai); |
| |
| for(i=1; i<2000; i++) { |
| const double xi = istrt + i; |
| const double xi1 = istrt + i - 1; |
| const double tmp = (a-1.0)*(xn+2.0*xi-1.0) + xi*(xi-beps); |
| const double b0_multiplier = (a+xi1-beps)*x/((xn+xi1)*(xi-beps)); |
| const double c0_multiplier_1 = (a+xi1)*x/((b+xi1)*xi); |
| const double c0_multiplier_2 = tmp / (xi*(b+xi1)*(a+xi1-beps)); |
| b0_val *= b0_multiplier; |
| b0_err += fabs(b0_multiplier) * b0_err + fabs(b0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON; |
| c0_val = c0_multiplier_1 * c0_val - c0_multiplier_2 * b0_val; |
| c0_err = fabs(c0_multiplier_1) * c0_err |
| + fabs(c0_multiplier_2) * b0_err |
| + fabs(c0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON |
| + fabs(b0_val * c0_multiplier_2) * 16.0 * 2.0 * GSL_DBL_EPSILON; |
| t_val = c0_val + xeps1_val*b0_val; |
| t_err = c0_err + fabs(xeps1_val)*b0_err; |
| t_err += fabs(b0_val*lnx) * dexprl.err; |
| t_err += fabs(b0_val)*xeps1_err; |
| dchu_val += t_val; |
| dchu_err += t_err; |
| if(fabs(t_val) < EPS*fabs(dchu_val)) break; |
| } |
| |
| result->val = dchu_val; |
| result->err = 2.0 * dchu_err; |
| result->err += 2.0 * fabs(t_val); |
| result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val); |
| result->err *= 2.0; /* FIXME: fudge factor */ |
| |
| if(i >= 2000) { |
| GSL_ERROR ("error", GSL_EMAXITER); |
| } |
| else { |
| return stat_all; |
| } |
| } |
| else { |
| /* |
| C X**(-BEPS) IS VERY DIFFERENT FROM 1.0, SO THE |
| C STRAIGHTFORWARD FORMULATION IS STABLE. |
| */ |
| int i; |
| double dchu_val; |
| double dchu_err; |
| double t_val; |
| double t_err; |
| gsl_sf_result dgamrbxi; |
| int stat_dgamrbxi = gsl_sf_gammainv_e(b+xi, &dgamrbxi); |
| double a0_val = factor_val * pochai.val * dgamrbxi.val * gamri1.val / beps; |
| double a0_err = fabs(factor_val * pochai.val * dgamrbxi.val / beps) * gamri1.err |
| + fabs(factor_val * pochai.val * gamri1.val / beps) * dgamrbxi.err |
| + fabs(factor_val * dgamrbxi.val * gamri1.val / beps) * pochai.err |
| + fabs(pochai.val * dgamrbxi.val * gamri1.val / beps) * factor_err |
| + 2.0 * GSL_DBL_EPSILON * fabs(a0_val); |
| stat_all = GSL_ERROR_SELECT_2(stat_all, stat_dgamrbxi); |
| |
| b0_val = xeps * b0_val / beps; |
| b0_err = fabs(xeps / beps) * b0_err + 4.0 * GSL_DBL_EPSILON * fabs(b0_val); |
| dchu_val = sum.val + a0_val - b0_val; |
| dchu_err = sum.err + a0_err + b0_err |
| + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(a0_val) + fabs(b0_val)); |
| |
| for(i=1; i<2000; i++) { |
| double xi = istrt + i; |
| double xi1 = istrt + i - 1; |
| double a0_multiplier = (a+xi1)*x/((b+xi1)*xi); |
| double b0_multiplier = (a+xi1-beps)*x/((aintb+xi1)*(xi-beps)); |
| a0_val *= a0_multiplier; |
| a0_err += fabs(a0_multiplier) * a0_err; |
| b0_val *= b0_multiplier; |
| b0_err += fabs(b0_multiplier) * b0_err; |
| t_val = a0_val - b0_val; |
| t_err = a0_err + b0_err; |
| dchu_val += t_val; |
| dchu_err += t_err; |
| if(fabs(t_val) < EPS*fabs(dchu_val)) break; |
| } |
| |
| result->val = dchu_val; |
| result->err = 2.0 * dchu_err; |
| result->err += 2.0 * fabs(t_val); |
| result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val); |
| result->err *= 2.0; /* FIXME: fudge factor */ |
| |
| if(i >= 2000) { |
| GSL_ERROR ("error", GSL_EMAXITER); |
| } |
| else { |
| return stat_all; |
| } |
| } |
| } |
| } |
| |
| |
| /* Assumes b > 0 and x > 0. |
| */ |
| static |
| int |
| hyperg_U_small_ab(const double a, const double b, const double x, gsl_sf_result * result) |
| { |
| if(a == -1.0) { |
| /* U(-1,c+1,x) = Laguerre[c,0,x] = -b + x |
| */ |
| result->val = -b + x; |
| result->err = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x)); |
| result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
| return GSL_SUCCESS; |
| } |
| else if(a == 0.0) { |
| /* U(0,c+1,x) = Laguerre[c,0,x] = 1 |
| */ |
| result->val = 1.0; |
| result->err = 0.0; |
| return GSL_SUCCESS; |
| } |
| else if(ASYMP_EVAL_OK(a,b,x)) { |
| double p = pow(x, -a); |
| gsl_sf_result asymp; |
| int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp); |
| result->val = asymp.val * p; |
| result->err = asymp.err * p; |
| result->err += fabs(asymp.val) * GSL_DBL_EPSILON * fabs(a) * p; |
| result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
| return stat_asymp; |
| } |
| else { |
| return hyperg_U_series(a, b, x, result); |
| } |
| } |
| |
| |
| /* Assumes b > 0 and x > 0. |
| */ |
| static |
| int |
| hyperg_U_small_a_bgt0(const double a, const double b, const double x, |
| gsl_sf_result * result, |
| double * ln_multiplier |
| ) |
| { |
| if(a == 0.0) { |
| result->val = 1.0; |
| result->err = 1.0; |
| *ln_multiplier = 0.0; |
| return GSL_SUCCESS; |
| } |
| else if( (b > 5000.0 && x < 0.90 * fabs(b)) |
| || (b > 500.0 && x < 0.50 * fabs(b)) |
| ) { |
| int stat = gsl_sf_hyperg_U_large_b_e(a, b, x, result, ln_multiplier); |
| if(stat == GSL_EOVRFLW) |
| return GSL_SUCCESS; |
| else |
| return stat; |
| } |
| else if(b > 15.0) { |
| /* Recurse up from b near 1. |
| */ |
| double eps = b - floor(b); |
| double b0 = 1.0 + eps; |
| gsl_sf_result r_Ubm1; |
| gsl_sf_result r_Ub; |
| int stat_0 = hyperg_U_small_ab(a, b0, x, &r_Ubm1); |
| int stat_1 = hyperg_U_small_ab(a, b0+1.0, x, &r_Ub); |
| double Ubm1 = r_Ubm1.val; |
| double Ub = r_Ub.val; |
| double Ubp1; |
| double bp; |
| |
| for(bp = b0+1.0; bp<b-0.1; bp += 1.0) { |
| Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x; |
| Ubm1 = Ub; |
| Ub = Ubp1; |
| } |
| result->val = Ub; |
| result->err = (fabs(r_Ubm1.err/r_Ubm1.val) + fabs(r_Ub.err/r_Ub.val)) * fabs(Ub); |
| result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-b0)+1.0) * fabs(Ub); |
| *ln_multiplier = 0.0; |
| return GSL_ERROR_SELECT_2(stat_0, stat_1); |
| } |
| else { |
| *ln_multiplier = 0.0; |
| return hyperg_U_small_ab(a, b, x, result); |
| } |
| } |
| |
| |
| /* We use this to keep track of large |
| * dynamic ranges in the recursions. |
| * This can be important because sometimes |
| * we want to calculate a very large and |
| * a very small number and the answer is |
| * the product, of order 1. This happens, |
| * for instance, when we apply a Kummer |
| * transform to make b positive and |
| * both x and b are large. |
| */ |
| #define RESCALE_2(u0,u1,factor,count) \ |
| do { \ |
| double au0 = fabs(u0); \ |
| if(au0 > factor) { \ |
| u0 /= factor; \ |
| u1 /= factor; \ |
| count++; \ |
| } \ |
| else if(au0 < 1.0/factor) { \ |
| u0 *= factor; \ |
| u1 *= factor; \ |
| count--; \ |
| } \ |
| } while (0) |
| |
| |
| /* Specialization to b >= 1, for integer parameters. |
| * Assumes x > 0. |
| */ |
| static |
| int |
| hyperg_U_int_bge1(const int a, const int b, const double x, |
| gsl_sf_result_e10 * result) |
| { |
| if(a == 0) { |
| result->val = 1.0; |
| result->err = 0.0; |
| result->e10 = 0; |
| return GSL_SUCCESS; |
| } |
| else if(a == -1) { |
| result->val = -b + x; |
| result->err = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x)); |
| result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val); |
| result->e10 = 0; |
| return GSL_SUCCESS; |
| } |
| else if(b == a + 1) { |
| /* U(a,a+1,x) = x^(-a) |
| */ |
| return gsl_sf_exp_e10_e(-a*log(x), result); |
| } |
| else if(ASYMP_EVAL_OK(a,b,x)) { |
| const double ln_pre_val = -a*log(x); |
| const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val); |
| gsl_sf_result asymp; |
| int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp); |
| int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err, |
| asymp.val, asymp.err, |
| result); |
| return GSL_ERROR_SELECT_2(stat_e, stat_asymp); |
| } |
| else if(SERIES_EVAL_OK(a,b,x)) { |
| gsl_sf_result ser; |
| const int stat_ser = hyperg_U_series(a, b, x, &ser); |
| result->val = ser.val; |
| result->err = ser.err; |
| result->e10 = 0; |
| return stat_ser; |
| } |
| else if(a < 0) { |
| /* Recurse backward from a = -1,0. |
| */ |
| int scale_count = 0; |
| const double scale_factor = GSL_SQRT_DBL_MAX; |
| gsl_sf_result lnm; |
| gsl_sf_result y; |
| double lnscale; |
| double Uap1 = 1.0; /* U(0,b,x) */ |
| double Ua = -b + x; /* U(-1,b,x) */ |
| double Uam1; |
| int ap; |
| |
| for(ap=-1; ap>a; ap--) { |
| Uam1 = ap*(b-ap-1.0)*Uap1 + (x+2.0*ap-b)*Ua; |
| Uap1 = Ua; |
| Ua = Uam1; |
| RESCALE_2(Ua,Uap1,scale_factor,scale_count); |
| } |
| |
| lnscale = log(scale_factor); |
| lnm.val = scale_count*lnscale; |
| lnm.err = 2.0 * GSL_DBL_EPSILON * fabs(lnm.val); |
| y.val = Ua; |
| y.err = 4.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(Ua); |
| return gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result); |
| } |
| else if(b >= 2.0*a + x) { |
| /* Recurse forward from a = 0,1. |
| */ |
| int scale_count = 0; |
| const double scale_factor = GSL_SQRT_DBL_MAX; |
| gsl_sf_result r_Ua; |
| gsl_sf_result lnm; |
| gsl_sf_result y; |
| double lnscale; |
| double lm; |
| int stat_1 = hyperg_U_small_a_bgt0(1.0, b, x, &r_Ua, &lm); /* U(1,b,x) */ |
| int stat_e; |
| double Uam1 = 1.0; /* U(0,b,x) */ |
| double Ua = r_Ua.val; |
| double Uap1; |
| int ap; |
| |
| Uam1 *= exp(-lm); |
| |
| for(ap=1; ap<a; ap++) { |
| Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b)); |
| Uam1 = Ua; |
| Ua = Uap1; |
| RESCALE_2(Ua,Uam1,scale_factor,scale_count); |
| } |
| |
| lnscale = log(scale_factor); |
| lnm.val = lm + scale_count * lnscale; |
| lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm) + fabs(scale_count*lnscale)); |
| y.val = Ua; |
| y.err = fabs(r_Ua.err/r_Ua.val) * fabs(Ua); |
| y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ua); |
| stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result); |
| return GSL_ERROR_SELECT_2(stat_e, stat_1); |
| } |
| else { |
| if(b <= x) { |
| /* Recurse backward either to the b=a+1 line |
| * or to a=0, whichever we hit. |
| */ |
| const double scale_factor = GSL_SQRT_DBL_MAX; |
| int scale_count = 0; |
| int stat_CF1; |
| double ru; |
| int CF1_count; |
| int a_target; |
| double lnU_target; |
| double Ua; |
| double Uap1; |
| double Uam1; |
| int ap; |
| |
| if(b < a + 1) { |
| a_target = b-1; |
| lnU_target = -a_target*log(x); |
| } |
| else { |
| a_target = 0; |
| lnU_target = 0.0; |
| } |
| |
| stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count); |
| |
| Ua = 1.0; |
| Uap1 = ru/a * Ua; |
| for(ap=a; ap>a_target; ap--) { |
| Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1); |
| Uap1 = Ua; |
| Ua = Uam1; |
| RESCALE_2(Ua,Uap1,scale_factor,scale_count); |
| } |
| |
| if(Ua == 0.0) { |
| result->val = 0.0; |
| result->err = 0.0; |
| result->e10 = 0; |
| GSL_ERROR ("error", GSL_EZERODIV); |
| } |
| else { |
| double lnscl = -scale_count*log(scale_factor); |
| double lnpre_val = lnU_target + lnscl; |
| double lnpre_err = 2.0 * GSL_DBL_EPSILON * (fabs(lnU_target) + fabs(lnscl)); |
| double oUa_err = 2.0 * (fabs(a_target-a) + CF1_count + 1.0) * GSL_DBL_EPSILON * fabs(1.0/Ua); |
| int stat_e = gsl_sf_exp_mult_err_e10_e(lnpre_val, lnpre_err, |
| 1.0/Ua, oUa_err, |
| result); |
| return GSL_ERROR_SELECT_2(stat_e, stat_CF1); |
| } |
| } |
| else { |
| /* Recurse backward to near the b=2a+x line, then |
| * determine normalization by either direct evaluation |
| * or by a forward recursion. The direct evaluation |
| * is needed when x is small (which is precisely |
| * when it is easy to do). |
| */ |
| const double scale_factor = GSL_SQRT_DBL_MAX; |
| int scale_count_for = 0; |
| int scale_count_bck = 0; |
| int a0 = 1; |
| int a1 = a0 + ceil(0.5*(b-x) - a0); |
| double Ua1_bck_val; |
| double Ua1_bck_err; |
| double Ua1_for_val; |
| double Ua1_for_err; |
| int stat_for; |
| int stat_bck; |
| gsl_sf_result lm_for; |
| |
| { |
| /* Recurse back to determine U(a1,b), sans normalization. |
| */ |
| double ru; |
| int CF1_count; |
| int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count); |
| double Ua = 1.0; |
| double Uap1 = ru/a * Ua; |
| double Uam1; |
| int ap; |
| for(ap=a; ap>a1; ap--) { |
| Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1); |
| Uap1 = Ua; |
| Ua = Uam1; |
| RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck); |
| } |
| Ua1_bck_val = Ua; |
| Ua1_bck_err = 2.0 * GSL_DBL_EPSILON * (fabs(a1-a)+CF1_count+1.0) * fabs(Ua); |
| stat_bck = stat_CF1; |
| } |
| |
| if(b == 2*a1 && a1 > 1) { |
| /* This can happen when x is small, which is |
| * precisely when we need to be careful with |
| * this evaluation. |
| */ |
| hyperg_lnU_beq2a((double)a1, x, &lm_for); |
| Ua1_for_val = 1.0; |
| Ua1_for_err = 0.0; |
| stat_for = GSL_SUCCESS; |
| } |
| else if(b == 2*a1 - 1 && a1 > 1) { |
| /* Similar to the above. Happens when x is small. |
| * Use |
| * U(a,2a-1) = (x U(a,2a) - U(a-1,2(a-1))) / (2a - 2) |
| */ |
| gsl_sf_result lnU00, lnU12; |
| gsl_sf_result U00, U12; |
| hyperg_lnU_beq2a(a1-1.0, x, &lnU00); |
| hyperg_lnU_beq2a(a1, x, &lnU12); |
| if(lnU00.val > lnU12.val) { |
| lm_for.val = lnU00.val; |
| lm_for.err = lnU00.err; |
| U00.val = 1.0; |
| U00.err = 0.0; |
| gsl_sf_exp_err_e(lnU12.val - lm_for.val, lnU12.err + lm_for.err, &U12); |
| } |
| else { |
| lm_for.val = lnU12.val; |
| lm_for.err = lnU12.err; |
| U12.val = 1.0; |
| U12.err = 0.0; |
| gsl_sf_exp_err_e(lnU00.val - lm_for.val, lnU00.err + lm_for.err, &U00); |
| } |
| Ua1_for_val = (x * U12.val - U00.val) / (2.0*a1 - 2.0); |
| Ua1_for_err = (fabs(x)*U12.err + U00.err) / fabs(2.0*a1 - 2.0); |
| Ua1_for_err += 2.0 * GSL_DBL_EPSILON * fabs(Ua1_for_val); |
| stat_for = GSL_SUCCESS; |
| } |
| else { |
| /* Recurse forward to determine U(a1,b) with |
| * absolute normalization. |
| */ |
| gsl_sf_result r_Ua; |
| double Uam1 = 1.0; /* U(a0-1,b,x) = U(0,b,x) */ |
| double Ua; |
| double Uap1; |
| int ap; |
| double lm_for_local; |
| stat_for = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_for_local); /* U(1,b,x) */ |
| Ua = r_Ua.val; |
| Uam1 *= exp(-lm_for_local); |
| lm_for.val = lm_for_local; |
| lm_for.err = 0.0; |
| |
| for(ap=a0; ap<a1; ap++) { |
| Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b)); |
| Uam1 = Ua; |
| Ua = Uap1; |
| RESCALE_2(Ua,Uam1,scale_factor,scale_count_for); |
| } |
| Ua1_for_val = Ua; |
| Ua1_for_err = fabs(Ua) * fabs(r_Ua.err/r_Ua.val); |
| Ua1_for_err += 2.0 * GSL_DBL_EPSILON * (fabs(a1-a0)+1.0) * fabs(Ua1_for_val); |
| } |
| |
| /* Now do the matching to produce the final result. |
| */ |
| if(Ua1_bck_val == 0.0) { |
| result->val = 0.0; |
| result->err = 0.0; |
| result->e10 = 0; |
| GSL_ERROR ("error", GSL_EZERODIV); |
| } |
| else if(Ua1_for_val == 0.0) { |
| /* Should never happen. */ |
| UNDERFLOW_ERROR_E10(result); |
| } |
| else { |
| double lns = (scale_count_for - scale_count_bck)*log(scale_factor); |
| double ln_for_val = log(fabs(Ua1_for_val)); |
| double ln_for_err = GSL_DBL_EPSILON + fabs(Ua1_for_err/Ua1_for_val); |
| double ln_bck_val = log(fabs(Ua1_bck_val)); |
| double ln_bck_err = GSL_DBL_EPSILON + fabs(Ua1_bck_err/Ua1_bck_val); |
| double lnr_val = lm_for.val + ln_for_val - ln_bck_val + lns; |
| double lnr_err = lm_for.err + ln_for_err + ln_bck_err |
| + 2.0 * GSL_DBL_EPSILON * (fabs(lm_for.val) + fabs(ln_for_val) + fabs(ln_bck_val) + fabs(lns)); |
| double sgn = GSL_SIGN(Ua1_for_val) * GSL_SIGN(Ua1_bck_val); |
| int stat_e = gsl_sf_exp_err_e10_e(lnr_val, lnr_err, result); |
| result->val *= sgn; |
| return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for); |
| } |
| } |
| } |
| } |
| |
| |
| /* Handle b >= 1 for generic a,b values. |
| */ |
| static |
| int |
| hyperg_U_bge1(const double a, const double b, const double x, |
| gsl_sf_result_e10 * result) |
| { |
| const double rinta = floor(a+0.5); |
| const int a_neg_integer = (a < 0.0 && fabs(a - rinta) < INT_THRESHOLD); |
| |
| if(a == 0.0) { |
| result->val = 1.0; |
| result->err = 0.0; |
| result->e10 = 0; |
| return GSL_SUCCESS; |
| } |
| else if(a_neg_integer && fabs(rinta) < INT_MAX) { |
| /* U(-n,b,x) = (-1)^n n! Laguerre[n,b-1,x] |
| */ |
| const int n = -(int)rinta; |
| const double sgn = (GSL_IS_ODD(n) ? -1.0 : 1.0); |
| gsl_sf_result lnfact; |
| gsl_sf_result L; |
| const int stat_L = gsl_sf_laguerre_n_e(n, b-1.0, x, &L); |
| gsl_sf_lnfact_e(n, &lnfact); |
| { |
| const int stat_e = gsl_sf_exp_mult_err_e10_e(lnfact.val, lnfact.err, |
| sgn*L.val, L.err, |
| result); |
| return GSL_ERROR_SELECT_2(stat_e, stat_L); |
| } |
| } |
| else if(ASYMP_EVAL_OK(a,b,x)) { |
| const double ln_pre_val = -a*log(x); |
| const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val); |
| gsl_sf_result asymp; |
| int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp); |
| int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err, |
| asymp.val, asymp.err, |
| result); |
| return GSL_ERROR_SELECT_2(stat_e, stat_asymp); |
| } |
| else if(fabs(a) <= 1.0) { |
| gsl_sf_result rU; |
| double ln_multiplier; |
| int stat_U = hyperg_U_small_a_bgt0(a, b, x, &rU, &ln_multiplier); |
| int stat_e = gsl_sf_exp_mult_err_e10_e(ln_multiplier, 2.0*GSL_DBL_EPSILON*fabs(ln_multiplier), rU.val, rU.err, result); |
| return GSL_ERROR_SELECT_2(stat_U, stat_e); |
| } |
| else if(SERIES_EVAL_OK(a,b,x)) { |
| gsl_sf_result ser; |
| const int stat_ser = hyperg_U_series(a, b, x, &ser); |
| result->val = ser.val; |
| result->err = ser.err; |
| result->e10 = 0; |
| return stat_ser; |
| } |
| else if(a < 0.0) { |
| /* Recurse backward on a and then upward on b. |
| */ |
| const double scale_factor = GSL_SQRT_DBL_MAX; |
| const double a0 = a - floor(a) - 1.0; |
| const double b0 = b - floor(b) + 1.0; |
| int scale_count = 0; |
| double lm_0, lm_1; |
| double lm_max; |
| gsl_sf_result r_Uap1; |
| gsl_sf_result r_Ua; |
| int stat_0 = hyperg_U_small_a_bgt0(a0+1.0, b0, x, &r_Uap1, &lm_0); |
| int stat_1 = hyperg_U_small_a_bgt0(a0, b0, x, &r_Ua, &lm_1); |
| int stat_e; |
| double Uap1 = r_Uap1.val; |
| double Ua = r_Ua.val; |
| double Uam1; |
| double ap; |
| lm_max = GSL_MAX(lm_0, lm_1); |
| Uap1 *= exp(lm_0-lm_max); |
| Ua *= exp(lm_1-lm_max); |
| |
| /* Downward recursion on a. |
| */ |
| for(ap=a0; ap>a+0.1; ap -= 1.0) { |
| Uam1 = ap*(b0-ap-1.0)*Uap1 + (x+2.0*ap-b0)*Ua; |
| Uap1 = Ua; |
| Ua = Uam1; |
| RESCALE_2(Ua,Uap1,scale_factor,scale_count); |
| } |
| |
| if(b < 2.0) { |
| /* b == b0, so no recursion necessary |
| */ |
| const double lnscale = log(scale_factor); |
| gsl_sf_result lnm; |
| gsl_sf_result y; |
| lnm.val = lm_max + scale_count * lnscale; |
| lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + scale_count * fabs(lnscale)); |
| y.val = Ua; |
| y.err = fabs(r_Uap1.err/r_Uap1.val) * fabs(Ua); |
| y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua); |
| y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua); |
| y.err *= fabs(lm_0-lm_max) + 1.0; |
| y.err *= fabs(lm_1-lm_max) + 1.0; |
| stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result); |
| } |
| else { |
| /* Upward recursion on b. |
| */ |
| const double err_mult = fabs(b-b0) + fabs(a-a0) + 1.0; |
| const double lnscale = log(scale_factor); |
| gsl_sf_result lnm; |
| gsl_sf_result y; |
| |
| double Ubm1 = Ua; /* U(a,b0) */ |
| double Ub = (a*(b0-a-1.0)*Uap1 + (a+x)*Ua)/x; /* U(a,b0+1) */ |
| double Ubp1; |
| double bp; |
| for(bp=b0+1.0; bp<b-0.1; bp += 1.0) { |
| Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x; |
| Ubm1 = Ub; |
| Ub = Ubp1; |
| RESCALE_2(Ub,Ubm1,scale_factor,scale_count); |
| } |
| |
| lnm.val = lm_max + scale_count * lnscale; |
| lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale)); |
| y.val = Ub; |
| y.err = 2.0 * err_mult * fabs(r_Uap1.err/r_Uap1.val) * fabs(Ub); |
| y.err += 2.0 * err_mult * fabs(r_Ua.err/r_Ua.val) * fabs(Ub); |
| y.err += 2.0 * GSL_DBL_EPSILON * err_mult * fabs(Ub); |
| y.err *= fabs(lm_0-lm_max) + 1.0; |
| y.err *= fabs(lm_1-lm_max) + 1.0; |
| stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result); |
| } |
| return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1); |
| } |
| else if(b >= 2*a + x) { |
| /* Recurse forward from a near zero. |
| * Note that we cannot cross the singularity at |
| * the line b=a+1, because the only way we could |
| * be in that little wedge is if a < 1. But we |
| * have already dealt with the small a case. |
| */ |
| int scale_count = 0; |
| const double a0 = a - floor(a); |
| const double scale_factor = GSL_SQRT_DBL_MAX; |
| double lnscale; |
| double lm_0, lm_1, lm_max; |
| gsl_sf_result r_Uam1; |
| gsl_sf_result r_Ua; |
| int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0); |
| int stat_1 = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_1); |
| int stat_e; |
| gsl_sf_result lnm; |
| gsl_sf_result y; |
| double Uam1 = r_Uam1.val; |
| double Ua = r_Ua.val; |
| double Uap1; |
| double ap; |
| lm_max = GSL_MAX(lm_0, lm_1); |
| Uam1 *= exp(lm_0-lm_max); |
| Ua *= exp(lm_1-lm_max); |
| |
| for(ap=a0; ap<a-0.1; ap += 1.0) { |
| Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b)); |
| Uam1 = Ua; |
| Ua = Uap1; |
| RESCALE_2(Ua,Uam1,scale_factor,scale_count); |
| } |
| |
| lnscale = log(scale_factor); |
| lnm.val = lm_max + scale_count * lnscale; |
| lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale)); |
| y.val = Ua; |
| y.err = fabs(r_Uam1.err/r_Uam1.val) * fabs(Ua); |
| y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua); |
| y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua); |
| y.err *= fabs(lm_0-lm_max) + 1.0; |
| y.err *= fabs(lm_1-lm_max) + 1.0; |
| stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result); |
| return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1); |
| } |
| else { |
| if(b <= x) { |
| /* Recurse backward to a near zero. |
| */ |
| const double a0 = a - floor(a); |
| const double scale_factor = GSL_SQRT_DBL_MAX; |
| int scale_count = 0; |
| gsl_sf_result lnm; |
| gsl_sf_result y; |
| double lnscale; |
| double lm_0; |
| double Uap1; |
| double Ua; |
| double Uam1; |
| gsl_sf_result U0; |
| double ap; |
| double ru; |
| double r; |
| int CF1_count; |
| int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count); |
| int stat_U0; |
| int stat_e; |
| r = ru/a; |
| Ua = GSL_SQRT_DBL_MIN; |
| Uap1 = r * Ua; |
| for(ap=a; ap>a0+0.1; ap -= 1.0) { |
| Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1); |
| Uap1 = Ua; |
| Ua = Uam1; |
| RESCALE_2(Ua,Uap1,scale_factor,scale_count); |
| } |
| |
| stat_U0 = hyperg_U_small_a_bgt0(a0, b, x, &U0, &lm_0); |
| |
| lnscale = log(scale_factor); |
| lnm.val = lm_0 - scale_count * lnscale; |
| lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_0) + fabs(scale_count * lnscale)); |
| y.val = GSL_SQRT_DBL_MIN*(U0.val/Ua); |
| y.err = GSL_SQRT_DBL_MIN*(U0.err/fabs(Ua)); |
| y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a0-a) + CF1_count + 1.0) * fabs(y.val); |
| stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result); |
| return GSL_ERROR_SELECT_3(stat_e, stat_U0, stat_CF1); |
| } |
| else { |
| /* Recurse backward to near the b=2a+x line, then |
| * forward from a near zero to get the normalization. |
| */ |
| int scale_count_for = 0; |
| int scale_count_bck = 0; |
| const double scale_factor = GSL_SQRT_DBL_MAX; |
| const double eps = a - floor(a); |
| const double a0 = ( eps == 0.0 ? 1.0 : eps ); |
| const double a1 = a0 + ceil(0.5*(b-x) - a0); |
| gsl_sf_result lnm; |
| gsl_sf_result y; |
| double lm_for; |
| double lnscale; |
| double Ua1_bck; |
| double Ua1_for; |
| int stat_for; |
| int stat_bck; |
| int stat_e; |
| int CF1_count; |
| |
| { |
| /* Recurse back to determine U(a1,b), sans normalization. |
| */ |
| double Uap1; |
| double Ua; |
| double Uam1; |
| double ap; |
| double ru; |
| double r; |
| int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count); |
| r = ru/a; |
| Ua = GSL_SQRT_DBL_MIN; |
| Uap1 = r * Ua; |
| for(ap=a; ap>a1+0.1; ap -= 1.0) { |
| Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1); |
| Uap1 = Ua; |
| Ua = Uam1; |
| RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck); |
| } |
| Ua1_bck = Ua; |
| stat_bck = stat_CF1; |
| } |
| { |
| /* Recurse forward to determine U(a1,b) with |
| * absolute normalization. |
| */ |
| gsl_sf_result r_Uam1; |
| gsl_sf_result r_Ua; |
| double lm_0, lm_1; |
| int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0); |
| int stat_1 = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_1); |
| double Uam1 = r_Uam1.val; |
| double Ua = r_Ua.val; |
| double Uap1; |
| double ap; |
| |
| lm_for = GSL_MAX(lm_0, lm_1); |
| Uam1 *= exp(lm_0 - lm_for); |
| Ua *= exp(lm_1 - lm_for); |
| |
| for(ap=a0; ap<a1-0.1; ap += 1.0) { |
| Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b)); |
| Uam1 = Ua; |
| Ua = Uap1; |
| RESCALE_2(Ua,Uam1,scale_factor,scale_count_for); |
| } |
| Ua1_for = Ua; |
| stat_for = GSL_ERROR_SELECT_2(stat_0, stat_1); |
| } |
| |
| lnscale = log(scale_factor); |
| lnm.val = lm_for + (scale_count_for - scale_count_bck)*lnscale; |
| lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_for) + fabs(scale_count_for - scale_count_bck)*fabs(lnscale)); |
| y.val = GSL_SQRT_DBL_MIN*Ua1_for/Ua1_bck; |
| y.err = 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + CF1_count + 1.0) * fabs(y.val); |
| stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result); |
| return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for); |
| } |
| } |
| } |
| |
| |
| /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/ |
| |
| |
| int |
| gsl_sf_hyperg_U_int_e10_e(const int a, const int b, const double x, |
| gsl_sf_result_e10 * result) |
| { |
| /* CHECK_POINTER(result) */ |
| |
| if(x <= 0.0) { |
| DOMAIN_ERROR_E10(result); |
| } |
| else { |
| if(b >= 1) { |
| return hyperg_U_int_bge1(a, b, x, result); |
| } |
| else { |
| /* Use the reflection formula |
| * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x) |
| */ |
| gsl_sf_result_e10 U; |
| double ln_x = log(x); |
| int ap = 1 + a - b; |
| int bp = 2 - b; |
| int stat_e; |
| int stat_U = hyperg_U_int_bge1(ap, bp, x, &U); |
| double ln_pre_val = (1.0-b)*ln_x; |
| double ln_pre_err = 2.0 * GSL_DBL_EPSILON * (fabs(b)+1.0) * fabs(ln_x); |
| ln_pre_err += 2.0 * GSL_DBL_EPSILON * fabs(1.0-b); /* error in log(x) */ |
| stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err, |
| U.val, U.err, |
| result); |
| return GSL_ERROR_SELECT_2(stat_e, stat_U); |
| } |
| } |
| } |
| |
| |
| int |
| gsl_sf_hyperg_U_e10_e(const double a, const double b, const double x, |
| gsl_sf_result_e10 * result) |
| { |
| const double rinta = floor(a + 0.5); |
| const double rintb = floor(b + 0.5); |
| const int a_integer = ( fabs(a - rinta) < INT_THRESHOLD ); |
| const int b_integer = ( fabs(b - rintb) < INT_THRESHOLD ); |
| |
| /* CHECK_POINTER(result) */ |
| |
| if(x <= 0.0) { |
| DOMAIN_ERROR_E10(result); |
| } |
| else if(a == 0.0) { |
| result->val = 1.0; |
| result->err = 0.0; |
| result->e10 = 0; |
| return GSL_SUCCESS; |
| } |
| else if(a_integer && b_integer) { |
| return gsl_sf_hyperg_U_int_e10_e(rinta, rintb, x, result); |
| } |
| else { |
| if(b >= 1.0) { |
| /* Use b >= 1 function. |
| */ |
| return hyperg_U_bge1(a, b, x, result); |
| } |
| else { |
| /* Use the reflection formula |
| * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x) |
| */ |
| const double lnx = log(x); |
| const double ln_pre_val = (1.0-b)*lnx; |
| const double ln_pre_err = fabs(lnx) * 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(b)); |
| const double ap = 1.0 + a - b; |
| const double bp = 2.0 - b; |
| gsl_sf_result_e10 U; |
| int stat_U = hyperg_U_bge1(ap, bp, x, &U); |
| int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err, |
| U.val, U.err, |
| result); |
| return GSL_ERROR_SELECT_2(stat_e, stat_U); |
| } |
| } |
| } |
| |
| |
| int |
| gsl_sf_hyperg_U_int_e(const int a, const int b, const double x, gsl_sf_result * result) |
| { |
| gsl_sf_result_e10 re; |
| int stat_U = gsl_sf_hyperg_U_int_e10_e(a, b, x, &re); |
| int stat_c = gsl_sf_result_smash_e(&re, result); |
| return GSL_ERROR_SELECT_2(stat_c, stat_U); |
| } |
| |
| |
| int |
| gsl_sf_hyperg_U_e(const double a, const double b, const double x, gsl_sf_result * result) |
| { |
| gsl_sf_result_e10 re; |
| int stat_U = gsl_sf_hyperg_U_e10_e(a, b, x, &re); |
| int stat_c = gsl_sf_result_smash_e(&re, result); |
| return GSL_ERROR_SELECT_2(stat_c, stat_U); |
| } |
| |
| |
| /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/ |
| |
| #include "eval.h" |
| |
| double gsl_sf_hyperg_U_int(const int a, const int b, const double x) |
| { |
| EVAL_RESULT(gsl_sf_hyperg_U_int_e(a, b, x, &result)); |
| } |
| |
| double gsl_sf_hyperg_U(const double a, const double b, const double x) |
| { |
| EVAL_RESULT(gsl_sf_hyperg_U_e(a, b, x, &result)); |
| } |
