1 /* specfunc/hyperg_U.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 /* Author: G. Jungman */
21
22 #include "gsl__config.h"
23 #include "gsl_math.h"
24 #include "gsl_errno.h"
25 #include "gsl_sf_exp.h"
26 #include "gsl_sf_gamma.h"
27 #include "gsl_sf_bessel.h"
28 #include "gsl_sf_laguerre.h"
29 #include "gsl_sf_pow_int.h"
30 #include "gsl_sf_hyperg.h"
31
32 #include "gsl_specfunc__error.h"
33 #include "gsl_specfunc__hyperg.h"
34
35 #define INT_THRESHOLD (1000.0*GSL_DBL_EPSILON)
36
37 #define SERIES_EVAL_OK(a,b,x) ((fabs(a) < 5 && b < 5 && x < 2.0) || (fabs(a) < 10 && b < 10 && x < 1.0))
38
39 #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))
40
41
42 /* Log[U(a,2a,x)]
43 * [Abramowitz+stegun, 13.6.21]
44 * Assumes x > 0, a > 1/2.
45 */
46 static
47 int
hyperg_lnU_beq2a(const double a,const double x,gsl_sf_result * result)48 hyperg_lnU_beq2a(const double a, const double x, gsl_sf_result * result)
49 {
50 const double lx = log(x);
51 const double nu = a - 0.5;
52 const double lnpre = 0.5*(x - M_LNPI) - nu*lx;
53 gsl_sf_result lnK;
54 gsl_sf_bessel_lnKnu_e(nu, 0.5*x, &lnK);
55 result->val = lnpre + lnK.val;
56 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(0.5*x) + 0.5*M_LNPI + fabs(nu*lx));
57 result->err += lnK.err;
58 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
59 return GSL_SUCCESS;
60 }
61
62
63 /* Evaluate u_{N+1}/u_N by Steed's continued fraction method.
64 *
65 * u_N := Gamma[a+N]/Gamma[a] U(a + N, b, x)
66 *
67 * u_{N+1}/u_N = (a+N) U(a+N+1,b,x)/U(a+N,b,x)
68 */
69 static
70 int
hyperg_U_CF1(const double a,const double b,const int N,const double x,double * result,int * count)71 hyperg_U_CF1(const double a, const double b, const int N, const double x,
72 double * result, int * count)
73 {
74 const double RECUR_BIG = GSL_SQRT_DBL_MAX;
75 const int maxiter = 20000;
76 int n = 1;
77 double Anm2 = 1.0;
78 double Bnm2 = 0.0;
79 double Anm1 = 0.0;
80 double Bnm1 = 1.0;
81 double a1 = -(a + N);
82 double b1 = (b - 2.0*a - x - 2.0*(N+1));
83 double An = b1*Anm1 + a1*Anm2;
84 double Bn = b1*Bnm1 + a1*Bnm2;
85 double an, bn;
86 double fn = An/Bn;
87
88 while(n < maxiter) {
89 double old_fn;
90 double del;
91 n++;
92 Anm2 = Anm1;
93 Bnm2 = Bnm1;
94 Anm1 = An;
95 Bnm1 = Bn;
96 an = -(a + N + n - b)*(a + N + n - 1.0);
97 bn = (b - 2.0*a - x - 2.0*(N+n));
98 An = bn*Anm1 + an*Anm2;
99 Bn = bn*Bnm1 + an*Bnm2;
100
101 if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
102 An /= RECUR_BIG;
103 Bn /= RECUR_BIG;
104 Anm1 /= RECUR_BIG;
105 Bnm1 /= RECUR_BIG;
106 Anm2 /= RECUR_BIG;
107 Bnm2 /= RECUR_BIG;
108 }
109
110 old_fn = fn;
111 fn = An/Bn;
112 del = old_fn/fn;
113
114 if(fabs(del - 1.0) < 10.0*GSL_DBL_EPSILON) break;
115 }
116
117 *result = fn;
118 *count = n;
119
120 if(n == maxiter)
121 GSL_ERROR ("error", GSL_EMAXITER);
122 else
123 return GSL_SUCCESS;
124 }
125
126
127 /* Large x asymptotic for x^a U(a,b,x)
128 * Based on SLATEC D9CHU() [W. Fullerton]
129 *
130 * Uses a rational approximation due to Luke.
131 * See [Luke, Algorithms for the Computation of Special Functions, p. 252]
132 * [Luke, Utilitas Math. (1977)]
133 *
134 * z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z)
135 *
136 * This assumes that a is not a negative integer and
137 * that 1+a-b is not a negative integer. If one of them
138 * is, then the 2F0 actually terminates, the above
139 * relation is an equality, and the sum should be
140 * evaluated directly [see below].
141 */
142 static
143 int
d9chu(const double a,const double b,const double x,gsl_sf_result * result)144 d9chu(const double a, const double b, const double x, gsl_sf_result * result)
145 {
146 const double EPS = 8.0 * GSL_DBL_EPSILON; /* EPS = 4.0D0*D1MACH(4) */
147 const int maxiter = 500;
148 double aa[4], bb[4];
149 int i;
150
151 double bp = 1.0 + a - b;
152 double ab = a*bp;
153 double ct2 = 2.0 * (x - ab);
154 double sab = a + bp;
155
156 double ct3 = sab + 1.0 + ab;
157 double anbn = ct3 + sab + 3.0;
158 double ct1 = 1.0 + 2.0*x/anbn;
159
160 bb[0] = 1.0;
161 aa[0] = 1.0;
162
163 bb[1] = 1.0 + 2.0*x/ct3;
164 aa[1] = 1.0 + ct2/ct3;
165
166 bb[2] = 1.0 + 6.0*ct1*x/ct3;
167 aa[2] = 1.0 + 6.0*ab/anbn + 3.0*ct1*ct2/ct3;
168
169 for(i=4; i<maxiter; i++) {
170 int j;
171 double c2;
172 double d1z;
173 double g1, g2, g3;
174 double x2i1 = 2*i - 3;
175 ct1 = x2i1/(x2i1-2.0);
176 anbn += x2i1 + sab;
177 ct2 = (x2i1 - 1.0)/anbn;
178 c2 = x2i1*ct2 - 1.0;
179 d1z = 2.0*x2i1*x/anbn;
180
181 ct3 = sab*ct2;
182 g1 = d1z + ct1*(c2+ct3);
183 g2 = d1z - c2;
184 g3 = ct1*(1.0 - ct3 - 2.0*ct2);
185
186 bb[3] = g1*bb[2] + g2*bb[1] + g3*bb[0];
187 aa[3] = g1*aa[2] + g2*aa[1] + g3*aa[0];
188
189 if(fabs(aa[3]*bb[0]-aa[0]*bb[3]) < EPS*fabs(bb[3]*bb[0])) break;
190
191 for(j=0; j<3; j++) {
192 aa[j] = aa[j+1];
193 bb[j] = bb[j+1];
194 }
195 }
196
197 result->val = aa[3]/bb[3];
198 result->err = 8.0 * GSL_DBL_EPSILON * fabs(result->val);
199
200 if(i == maxiter) {
201 GSL_ERROR ("error", GSL_EMAXITER);
202 }
203 else {
204 return GSL_SUCCESS;
205 }
206 }
207
208
209 /* Evaluate asymptotic for z^a U(a,b,z) ~ 2F0(a,1+a-b,-1/z)
210 * We check for termination of the 2F0 as a special case.
211 * Assumes x > 0.
212 * Also assumes a,b are not too large compared to x.
213 */
214 static
215 int
hyperg_zaU_asymp(const double a,const double b,const double x,gsl_sf_result * result)216 hyperg_zaU_asymp(const double a, const double b, const double x, gsl_sf_result *result)
217 {
218 const double ap = a;
219 const double bp = 1.0 + a - b;
220 const double rintap = floor(ap + 0.5);
221 const double rintbp = floor(bp + 0.5);
222 const int ap_neg_int = ( ap < 0.0 && fabs(ap - rintap) < INT_THRESHOLD );
223 const int bp_neg_int = ( bp < 0.0 && fabs(bp - rintbp) < INT_THRESHOLD );
224
225 if(ap_neg_int || bp_neg_int) {
226 /* Evaluate 2F0 polynomial.
227 */
228 double mxi = -1.0/x;
229 double nmax = -(int)(GSL_MIN(ap,bp) - 0.1);
230 double tn = 1.0;
231 double sum = 1.0;
232 double n = 1.0;
233 double sum_err = 0.0;
234 while(n <= nmax) {
235 double apn = (ap+n-1.0);
236 double bpn = (bp+n-1.0);
237 tn *= ((apn/n)*mxi)*bpn;
238 sum += tn;
239 sum_err += 2.0 * GSL_DBL_EPSILON * fabs(tn);
240 n += 1.0;
241 }
242 result->val = sum;
243 result->err = sum_err;
244 result->err += 2.0 * GSL_DBL_EPSILON * (fabs(nmax)+1.0) * fabs(sum);
245 return GSL_SUCCESS;
246 }
247 else {
248 return d9chu(a,b,x,result);
249 }
250 }
251
252
253 /* Evaluate finite sum which appears below.
254 */
255 static
256 int
hyperg_U_finite_sum(int N,double a,double b,double x,double xeps,gsl_sf_result * result)257 hyperg_U_finite_sum(int N, double a, double b, double x, double xeps,
258 gsl_sf_result * result)
259 {
260 int i;
261 double sum_val;
262 double sum_err;
263
264 if(N <= 0) {
265 double t_val = 1.0;
266 double t_err = 0.0;
267 gsl_sf_result poch;
268 int stat_poch;
269
270 sum_val = 1.0;
271 sum_err = 0.0;
272 for(i=1; i<= -N; i++) {
273 const double xi1 = i - 1;
274 const double mult = (a+xi1)*x/((b+xi1)*(xi1+1.0));
275 t_val *= mult;
276 t_err += fabs(mult) * t_err + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
277 sum_val += t_val;
278 sum_err += t_err;
279 }
280
281 stat_poch = gsl_sf_poch_e(1.0+a-b, -a, &poch);
282
283 result->val = sum_val * poch.val;
284 result->err = fabs(sum_val) * poch.err + sum_err * fabs(poch.val);
285 result->err += fabs(poch.val) * (fabs(N) + 2.0) * GSL_DBL_EPSILON * fabs(sum_val);
286 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
287 result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
288 return stat_poch;
289 }
290 else {
291 const int M = N - 2;
292 if(M < 0) {
293 result->val = 0.0;
294 result->err = 0.0;
295 return GSL_SUCCESS;
296 }
297 else {
298 gsl_sf_result gbm1;
299 gsl_sf_result gamr;
300 int stat_gbm1;
301 double t_val = 1.0;
302 double t_err = 0.0;
303
304 sum_val = 1.0;
305 sum_err = 0.0;
306 for(i=1; i<=M; i++) {
307 const double mult = (a-b+i)*x/((1.0-b+i)*i);
308 t_val *= mult;
309 t_err += t_err * fabs(mult) + fabs(t_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
310 sum_val += t_val;
311 sum_err += t_err;
312 }
313
314 stat_gbm1 = gsl_sf_gamma_e(b-1.0, &gbm1);
315 (void) gsl_sf_gammainv_e(a, &gamr);
316
317 if(stat_gbm1 == GSL_SUCCESS) {
318 gsl_sf_result powx1N;
319 int stat_p = gsl_sf_pow_int_e(x, 1-N, &powx1N);
320 double pe_val = powx1N.val * xeps;
321 double pe_err = powx1N.err * fabs(xeps) + 2.0 * GSL_DBL_EPSILON * fabs(pe_val);
322 double coeff_val = gbm1.val * gamr.val * pe_val;
323 double coeff_err = gbm1.err * fabs(gamr.val * pe_val)
324 + gamr.err * fabs(gbm1.val * pe_val)
325 + fabs(gbm1.val * gamr.val) * pe_err
326 + 2.0 * GSL_DBL_EPSILON * fabs(coeff_val);
327
328 result->val = sum_val * coeff_val;
329 result->err = fabs(sum_val) * coeff_err + sum_err * fabs(coeff_val);
330 result->err += 2.0 * GSL_DBL_EPSILON * (M+2.0) * fabs(result->val);
331 result->err *= 2.0; /* FIXME: fudge factor... why is the error estimate too small? */
332 return stat_p;
333 }
334 else {
335 result->val = 0.0;
336 result->err = 0.0;
337 return stat_gbm1;
338 }
339 }
340 }
341 }
342
343
344 /* Based on SLATEC DCHU() [W. Fullerton]
345 * Assumes x > 0.
346 * This is just a series summation method, and
347 * it is not good for large a.
348 *
349 * I patched up the window for 1+a-b near zero. [GJ]
350 */
351 static
352 int
hyperg_U_series(const double a,const double b,const double x,gsl_sf_result * result)353 hyperg_U_series(const double a, const double b, const double x, gsl_sf_result * result)
354 {
355 const double EPS = 2.0 * GSL_DBL_EPSILON; /* EPS = D1MACH(3) */
356 const double SQRT_EPS = M_SQRT2 * GSL_SQRT_DBL_EPSILON;
357
358 if(fabs(1.0 + a - b) < SQRT_EPS) {
359 /* Original Comment: ALGORITHM IS BAD WHEN 1+A-B IS NEAR ZERO FOR SMALL X
360 */
361 /* We can however do the following:
362 * U(a,b,x) = U(a,a+1,x) when 1+a-b=0
363 * and U(a,a+1,x) = x^(-a).
364 */
365 double lnr = -a * log(x);
366 int stat_e = gsl_sf_exp_e(lnr, result);
367 result->err += 2.0 * SQRT_EPS * fabs(result->val);
368 return stat_e;
369 }
370 else {
371 double aintb = ( b < 0.0 ? ceil(b-0.5) : floor(b+0.5) );
372 double beps = b - aintb;
373 int N = aintb;
374
375 double lnx = log(x);
376 double xeps = exp(-beps*lnx);
377
378 /* Evaluate finite sum.
379 */
380 gsl_sf_result sum;
381 int stat_sum = hyperg_U_finite_sum(N, a, b, x, xeps, &sum);
382
383
384 /* Evaluate infinite sum. */
385
386 int istrt = ( N < 1 ? 1-N : 0 );
387 double xi = istrt;
388
389 gsl_sf_result gamr;
390 gsl_sf_result powx;
391 int stat_gamr = gsl_sf_gammainv_e(1.0+a-b, &gamr);
392 int stat_powx = gsl_sf_pow_int_e(x, istrt, &powx);
393 double sarg = beps*M_PI;
394 double sfact = ( sarg != 0.0 ? sarg/sin(sarg) : 1.0 );
395 double factor_val = sfact * ( GSL_IS_ODD(N) ? -1.0 : 1.0 ) * gamr.val * powx.val;
396 double factor_err = fabs(gamr.val) * powx.err + fabs(powx.val) * gamr.err
397 + 2.0 * GSL_DBL_EPSILON * fabs(factor_val);
398
399 gsl_sf_result pochai;
400 gsl_sf_result gamri1;
401 gsl_sf_result gamrni;
402 int stat_pochai = gsl_sf_poch_e(a, xi, &pochai);
403 int stat_gamri1 = gsl_sf_gammainv_e(xi + 1.0, &gamri1);
404 int stat_gamrni = gsl_sf_gammainv_e(aintb + xi, &gamrni);
405 int stat_gam123 = GSL_ERROR_SELECT_3(stat_gamr, stat_gamri1, stat_gamrni);
406 int stat_gamall = GSL_ERROR_SELECT_4(stat_sum, stat_gam123, stat_pochai, stat_powx);
407
408 gsl_sf_result pochaxibeps;
409 gsl_sf_result gamrxi1beps;
410 int stat_pochaxibeps = gsl_sf_poch_e(a, xi-beps, &pochaxibeps);
411 int stat_gamrxi1beps = gsl_sf_gammainv_e(xi + 1.0 - beps, &gamrxi1beps);
412
413 int stat_all = GSL_ERROR_SELECT_3(stat_gamall, stat_pochaxibeps, stat_gamrxi1beps);
414
415 double b0_val = factor_val * pochaxibeps.val * gamrni.val * gamrxi1beps.val;
416 double b0_err = fabs(factor_val * pochaxibeps.val * gamrni.val) * gamrxi1beps.err
417 + fabs(factor_val * pochaxibeps.val * gamrxi1beps.val) * gamrni.err
418 + fabs(factor_val * gamrni.val * gamrxi1beps.val) * pochaxibeps.err
419 + fabs(pochaxibeps.val * gamrni.val * gamrxi1beps.val) * factor_err
420 + 2.0 * GSL_DBL_EPSILON * fabs(b0_val);
421
422 if(fabs(xeps-1.0) < 0.5) {
423 /*
424 C X**(-BEPS) IS CLOSE TO 1.0D0, SO WE MUST BE
425 C CAREFUL IN EVALUATING THE DIFFERENCES.
426 */
427 int i;
428 gsl_sf_result pch1ai;
429 gsl_sf_result pch1i;
430 gsl_sf_result poch1bxibeps;
431 int stat_pch1ai = gsl_sf_pochrel_e(a + xi, -beps, &pch1ai);
432 int stat_pch1i = gsl_sf_pochrel_e(xi + 1.0 - beps, beps, &pch1i);
433 int stat_poch1bxibeps = gsl_sf_pochrel_e(b+xi, -beps, &poch1bxibeps);
434 double c0_t1_val = beps*pch1ai.val*pch1i.val;
435 double c0_t1_err = fabs(beps) * fabs(pch1ai.val) * pch1i.err
436 + fabs(beps) * fabs(pch1i.val) * pch1ai.err
437 + 2.0 * GSL_DBL_EPSILON * fabs(c0_t1_val);
438 double c0_t2_val = -poch1bxibeps.val + pch1ai.val - pch1i.val + c0_t1_val;
439 double c0_t2_err = poch1bxibeps.err + pch1ai.err + pch1i.err + c0_t1_err
440 + 2.0 * GSL_DBL_EPSILON * fabs(c0_t2_val);
441 double c0_val = factor_val * pochai.val * gamrni.val * gamri1.val * c0_t2_val;
442 double c0_err = fabs(factor_val * pochai.val * gamrni.val * gamri1.val) * c0_t2_err
443 + fabs(factor_val * pochai.val * gamrni.val * c0_t2_val) * gamri1.err
444 + fabs(factor_val * pochai.val * gamri1.val * c0_t2_val) * gamrni.err
445 + fabs(factor_val * gamrni.val * gamri1.val * c0_t2_val) * pochai.err
446 + fabs(pochai.val * gamrni.val * gamri1.val * c0_t2_val) * factor_err
447 + 2.0 * GSL_DBL_EPSILON * fabs(c0_val);
448 /*
449 C XEPS1 = (1.0 - X**(-BEPS))/BEPS = (X**(-BEPS) - 1.0)/(-BEPS)
450 */
451 gsl_sf_result dexprl;
452 int stat_dexprl = gsl_sf_exprel_e(-beps*lnx, &dexprl);
453 double xeps1_val = lnx * dexprl.val;
454 double xeps1_err = 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(beps*lnx)) * fabs(dexprl.val)
455 + fabs(lnx) * dexprl.err
456 + 2.0 * GSL_DBL_EPSILON * fabs(xeps1_val);
457 double dchu_val = sum.val + c0_val + xeps1_val*b0_val;
458 double dchu_err = sum.err + c0_err
459 + fabs(xeps1_val)*b0_err + xeps1_err * fabs(b0_val)
460 + fabs(b0_val*lnx)*dexprl.err
461 + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(c0_val) + fabs(xeps1_val*b0_val));
462 double xn = N;
463 double t_val;
464 double t_err;
465
466 stat_all = GSL_ERROR_SELECT_5(stat_all, stat_dexprl, stat_poch1bxibeps, stat_pch1i, stat_pch1ai);
467
468 for(i=1; i<2000; i++) {
469 const double xi = istrt + i;
470 const double xi1 = istrt + i - 1;
471 const double tmp = (a-1.0)*(xn+2.0*xi-1.0) + xi*(xi-beps);
472 const double b0_multiplier = (a+xi1-beps)*x/((xn+xi1)*(xi-beps));
473 const double c0_multiplier_1 = (a+xi1)*x/((b+xi1)*xi);
474 const double c0_multiplier_2 = tmp / (xi*(b+xi1)*(a+xi1-beps));
475 b0_val *= b0_multiplier;
476 b0_err += fabs(b0_multiplier) * b0_err + fabs(b0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON;
477 c0_val = c0_multiplier_1 * c0_val - c0_multiplier_2 * b0_val;
478 c0_err = fabs(c0_multiplier_1) * c0_err
479 + fabs(c0_multiplier_2) * b0_err
480 + fabs(c0_val) * 8.0 * 2.0 * GSL_DBL_EPSILON
481 + fabs(b0_val * c0_multiplier_2) * 16.0 * 2.0 * GSL_DBL_EPSILON;
482 t_val = c0_val + xeps1_val*b0_val;
483 t_err = c0_err + fabs(xeps1_val)*b0_err;
484 t_err += fabs(b0_val*lnx) * dexprl.err;
485 t_err += fabs(b0_val)*xeps1_err;
486 dchu_val += t_val;
487 dchu_err += t_err;
488 if(fabs(t_val) < EPS*fabs(dchu_val)) break;
489 }
490
491 result->val = dchu_val;
492 result->err = 2.0 * dchu_err;
493 result->err += 2.0 * fabs(t_val);
494 result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);
495 result->err *= 2.0; /* FIXME: fudge factor */
496
497 if(i >= 2000) {
498 GSL_ERROR ("error", GSL_EMAXITER);
499 }
500 else {
501 return stat_all;
502 }
503 }
504 else {
505 /*
506 C X**(-BEPS) IS VERY DIFFERENT FROM 1.0, SO THE
507 C STRAIGHTFORWARD FORMULATION IS STABLE.
508 */
509 int i;
510 double dchu_val;
511 double dchu_err;
512 double t_val;
513 double t_err;
514 gsl_sf_result dgamrbxi;
515 int stat_dgamrbxi = gsl_sf_gammainv_e(b+xi, &dgamrbxi);
516 double a0_val = factor_val * pochai.val * dgamrbxi.val * gamri1.val / beps;
517 double a0_err = fabs(factor_val * pochai.val * dgamrbxi.val / beps) * gamri1.err
518 + fabs(factor_val * pochai.val * gamri1.val / beps) * dgamrbxi.err
519 + fabs(factor_val * dgamrbxi.val * gamri1.val / beps) * pochai.err
520 + fabs(pochai.val * dgamrbxi.val * gamri1.val / beps) * factor_err
521 + 2.0 * GSL_DBL_EPSILON * fabs(a0_val);
522 stat_all = GSL_ERROR_SELECT_2(stat_all, stat_dgamrbxi);
523
524 b0_val = xeps * b0_val / beps;
525 b0_err = fabs(xeps / beps) * b0_err + 4.0 * GSL_DBL_EPSILON * fabs(b0_val);
526 dchu_val = sum.val + a0_val - b0_val;
527 dchu_err = sum.err + a0_err + b0_err
528 + 2.0 * GSL_DBL_EPSILON * (fabs(sum.val) + fabs(a0_val) + fabs(b0_val));
529
530 for(i=1; i<2000; i++) {
531 double xi = istrt + i;
532 double xi1 = istrt + i - 1;
533 double a0_multiplier = (a+xi1)*x/((b+xi1)*xi);
534 double b0_multiplier = (a+xi1-beps)*x/((aintb+xi1)*(xi-beps));
535 a0_val *= a0_multiplier;
536 a0_err += fabs(a0_multiplier) * a0_err;
537 b0_val *= b0_multiplier;
538 b0_err += fabs(b0_multiplier) * b0_err;
539 t_val = a0_val - b0_val;
540 t_err = a0_err + b0_err;
541 dchu_val += t_val;
542 dchu_err += t_err;
543 if(fabs(t_val) < EPS*fabs(dchu_val)) break;
544 }
545
546 result->val = dchu_val;
547 result->err = 2.0 * dchu_err;
548 result->err += 2.0 * fabs(t_val);
549 result->err += 4.0 * GSL_DBL_EPSILON * (i+2.0) * fabs(dchu_val);
550 result->err *= 2.0; /* FIXME: fudge factor */
551
552 if(i >= 2000) {
553 GSL_ERROR ("error", GSL_EMAXITER);
554 }
555 else {
556 return stat_all;
557 }
558 }
559 }
560 }
561
562
563 /* Assumes b > 0 and x > 0.
564 */
565 static
566 int
hyperg_U_small_ab(const double a,const double b,const double x,gsl_sf_result * result)567 hyperg_U_small_ab(const double a, const double b, const double x, gsl_sf_result * result)
568 {
569 if(a == -1.0) {
570 /* U(-1,c+1,x) = Laguerre[c,0,x] = -b + x
571 */
572 result->val = -b + x;
573 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));
574 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
575 return GSL_SUCCESS;
576 }
577 else if(a == 0.0) {
578 /* U(0,c+1,x) = Laguerre[c,0,x] = 1
579 */
580 result->val = 1.0;
581 result->err = 0.0;
582 return GSL_SUCCESS;
583 }
584 else if(ASYMP_EVAL_OK(a,b,x)) {
585 double p = pow(x, -a);
586 gsl_sf_result asymp;
587 int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
588 result->val = asymp.val * p;
589 result->err = asymp.err * p;
590 result->err += fabs(asymp.val) * GSL_DBL_EPSILON * fabs(a) * p;
591 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
592 return stat_asymp;
593 }
594 else {
595 return hyperg_U_series(a, b, x, result);
596 }
597 }
598
599
600 /* Assumes b > 0 and x > 0.
601 */
602 static
603 int
hyperg_U_small_a_bgt0(const double a,const double b,const double x,gsl_sf_result * result,double * ln_multiplier)604 hyperg_U_small_a_bgt0(const double a, const double b, const double x,
605 gsl_sf_result * result,
606 double * ln_multiplier
607 )
608 {
609 if(a == 0.0) {
610 result->val = 1.0;
611 result->err = 1.0;
612 *ln_multiplier = 0.0;
613 return GSL_SUCCESS;
614 }
615 else if( (b > 5000.0 && x < 0.90 * fabs(b))
616 || (b > 500.0 && x < 0.50 * fabs(b))
617 ) {
618 int stat = gsl_sf_hyperg_U_large_b_e(a, b, x, result, ln_multiplier);
619 if(stat == GSL_EOVRFLW)
620 return GSL_SUCCESS;
621 else
622 return stat;
623 }
624 else if(b > 15.0) {
625 /* Recurse up from b near 1.
626 */
627 double eps = b - floor(b);
628 double b0 = 1.0 + eps;
629 gsl_sf_result r_Ubm1;
630 gsl_sf_result r_Ub;
631 int stat_0 = hyperg_U_small_ab(a, b0, x, &r_Ubm1);
632 int stat_1 = hyperg_U_small_ab(a, b0+1.0, x, &r_Ub);
633 double Ubm1 = r_Ubm1.val;
634 double Ub = r_Ub.val;
635 double Ubp1;
636 double bp;
637
638 for(bp = b0+1.0; bp<b-0.1; bp += 1.0) {
639 Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x;
640 Ubm1 = Ub;
641 Ub = Ubp1;
642 }
643 result->val = Ub;
644 result->err = (fabs(r_Ubm1.err/r_Ubm1.val) + fabs(r_Ub.err/r_Ub.val)) * fabs(Ub);
645 result->err += 2.0 * GSL_DBL_EPSILON * (fabs(b-b0)+1.0) * fabs(Ub);
646 *ln_multiplier = 0.0;
647 return GSL_ERROR_SELECT_2(stat_0, stat_1);
648 }
649 else {
650 *ln_multiplier = 0.0;
651 return hyperg_U_small_ab(a, b, x, result);
652 }
653 }
654
655
656 /* We use this to keep track of large
657 * dynamic ranges in the recursions.
658 * This can be important because sometimes
659 * we want to calculate a very large and
660 * a very small number and the answer is
661 * the product, of order 1. This happens,
662 * for instance, when we apply a Kummer
663 * transform to make b positive and
664 * both x and b are large.
665 */
666 #define RESCALE_2(u0,u1,factor,count) \
667 do { \
668 double au0 = fabs(u0); \
669 if(au0 > factor) { \
670 u0 /= factor; \
671 u1 /= factor; \
672 count++; \
673 } \
674 else if(au0 < 1.0/factor) { \
675 u0 *= factor; \
676 u1 *= factor; \
677 count--; \
678 } \
679 } while (0)
680
681
682 /* Specialization to b >= 1, for integer parameters.
683 * Assumes x > 0.
684 */
685 static
686 int
hyperg_U_int_bge1(const int a,const int b,const double x,gsl_sf_result_e10 * result)687 hyperg_U_int_bge1(const int a, const int b, const double x,
688 gsl_sf_result_e10 * result)
689 {
690 if(a == 0) {
691 result->val = 1.0;
692 result->err = 0.0;
693 result->e10 = 0;
694 return GSL_SUCCESS;
695 }
696 else if(a == -1) {
697 result->val = -b + x;
698 result->err = 2.0 * GSL_DBL_EPSILON * (fabs(b) + fabs(x));
699 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
700 result->e10 = 0;
701 return GSL_SUCCESS;
702 }
703 else if(b == a + 1) {
704 /* U(a,a+1,x) = x^(-a)
705 */
706 return gsl_sf_exp_e10_e(-a*log(x), result);
707 }
708 else if(ASYMP_EVAL_OK(a,b,x)) {
709 const double ln_pre_val = -a*log(x);
710 const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val);
711 gsl_sf_result asymp;
712 int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
713 int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err,
714 asymp.val, asymp.err,
715 result);
716 return GSL_ERROR_SELECT_2(stat_e, stat_asymp);
717 }
718 else if(SERIES_EVAL_OK(a,b,x)) {
719 gsl_sf_result ser;
720 const int stat_ser = hyperg_U_series(a, b, x, &ser);
721 result->val = ser.val;
722 result->err = ser.err;
723 result->e10 = 0;
724 return stat_ser;
725 }
726 else if(a < 0) {
727 /* Recurse backward from a = -1,0.
728 */
729 int scale_count = 0;
730 const double scale_factor = GSL_SQRT_DBL_MAX;
731 gsl_sf_result lnm;
732 gsl_sf_result y;
733 double lnscale;
734 double Uap1 = 1.0; /* U(0,b,x) */
735 double Ua = -b + x; /* U(-1,b,x) */
736 double Uam1;
737 int ap;
738
739 for(ap=-1; ap>a; ap--) {
740 Uam1 = ap*(b-ap-1.0)*Uap1 + (x+2.0*ap-b)*Ua;
741 Uap1 = Ua;
742 Ua = Uam1;
743 RESCALE_2(Ua,Uap1,scale_factor,scale_count);
744 }
745
746 lnscale = log(scale_factor);
747 lnm.val = scale_count*lnscale;
748 lnm.err = 2.0 * GSL_DBL_EPSILON * fabs(lnm.val);
749 y.val = Ua;
750 y.err = 4.0 * GSL_DBL_EPSILON * (fabs(a)+1.0) * fabs(Ua);
751 return gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
752 }
753 else if(b >= 2.0*a + x) {
754 /* Recurse forward from a = 0,1.
755 */
756 int scale_count = 0;
757 const double scale_factor = GSL_SQRT_DBL_MAX;
758 gsl_sf_result r_Ua;
759 gsl_sf_result lnm;
760 gsl_sf_result y;
761 double lnscale;
762 double lm;
763 int stat_1 = hyperg_U_small_a_bgt0(1.0, b, x, &r_Ua, &lm); /* U(1,b,x) */
764 int stat_e;
765 double Uam1 = 1.0; /* U(0,b,x) */
766 double Ua = r_Ua.val;
767 double Uap1;
768 int ap;
769
770 Uam1 *= exp(-lm);
771
772 for(ap=1; ap<a; ap++) {
773 Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
774 Uam1 = Ua;
775 Ua = Uap1;
776 RESCALE_2(Ua,Uam1,scale_factor,scale_count);
777 }
778
779 lnscale = log(scale_factor);
780 lnm.val = lm + scale_count * lnscale;
781 lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm) + fabs(scale_count*lnscale));
782 y.val = Ua;
783 y.err = fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
784 y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a) + 1.0) * fabs(Ua);
785 stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
786 return GSL_ERROR_SELECT_2(stat_e, stat_1);
787 }
788 else {
789 if(b <= x) {
790 /* Recurse backward either to the b=a+1 line
791 * or to a=0, whichever we hit.
792 */
793 const double scale_factor = GSL_SQRT_DBL_MAX;
794 int scale_count = 0;
795 int stat_CF1;
796 double ru;
797 int CF1_count;
798 int a_target;
799 double lnU_target;
800 double Ua;
801 double Uap1;
802 double Uam1;
803 int ap;
804
805 if(b < a + 1) {
806 a_target = b-1;
807 lnU_target = -a_target*log(x);
808 }
809 else {
810 a_target = 0;
811 lnU_target = 0.0;
812 }
813
814 stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
815
816 Ua = 1.0;
817 Uap1 = ru/a * Ua;
818 for(ap=a; ap>a_target; ap--) {
819 Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
820 Uap1 = Ua;
821 Ua = Uam1;
822 RESCALE_2(Ua,Uap1,scale_factor,scale_count);
823 }
824
825 if(Ua == 0.0) {
826 result->val = 0.0;
827 result->err = 0.0;
828 result->e10 = 0;
829 GSL_ERROR ("error", GSL_EZERODIV);
830 }
831 else {
832 double lnscl = -scale_count*log(scale_factor);
833 double lnpre_val = lnU_target + lnscl;
834 double lnpre_err = 2.0 * GSL_DBL_EPSILON * (fabs(lnU_target) + fabs(lnscl));
835 double oUa_err = 2.0 * (fabs(a_target-a) + CF1_count + 1.0) * GSL_DBL_EPSILON * fabs(1.0/Ua);
836 int stat_e = gsl_sf_exp_mult_err_e10_e(lnpre_val, lnpre_err,
837 1.0/Ua, oUa_err,
838 result);
839 return GSL_ERROR_SELECT_2(stat_e, stat_CF1);
840 }
841 }
842 else {
843 /* Recurse backward to near the b=2a+x line, then
844 * determine normalization by either direct evaluation
845 * or by a forward recursion. The direct evaluation
846 * is needed when x is small (which is precisely
847 * when it is easy to do).
848 */
849 const double scale_factor = GSL_SQRT_DBL_MAX;
850 int scale_count_for = 0;
851 int scale_count_bck = 0;
852 int a0 = 1;
853 int a1 = a0 + ceil(0.5*(b-x) - a0);
854 double Ua1_bck_val;
855 double Ua1_bck_err;
856 double Ua1_for_val;
857 double Ua1_for_err;
858 int stat_for;
859 int stat_bck;
860 gsl_sf_result lm_for;
861
862 {
863 /* Recurse back to determine U(a1,b), sans normalization.
864 */
865 double ru;
866 int CF1_count;
867 int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
868 double Ua = 1.0;
869 double Uap1 = ru/a * Ua;
870 double Uam1;
871 int ap;
872 for(ap=a; ap>a1; ap--) {
873 Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
874 Uap1 = Ua;
875 Ua = Uam1;
876 RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck);
877 }
878 Ua1_bck_val = Ua;
879 Ua1_bck_err = 2.0 * GSL_DBL_EPSILON * (fabs(a1-a)+CF1_count+1.0) * fabs(Ua);
880 stat_bck = stat_CF1;
881 }
882
883 if(b == 2*a1 && a1 > 1) {
884 /* This can happen when x is small, which is
885 * precisely when we need to be careful with
886 * this evaluation.
887 */
888 hyperg_lnU_beq2a((double)a1, x, &lm_for);
889 Ua1_for_val = 1.0;
890 Ua1_for_err = 0.0;
891 stat_for = GSL_SUCCESS;
892 }
893 else if(b == 2*a1 - 1 && a1 > 1) {
894 /* Similar to the above. Happens when x is small.
895 * Use
896 * U(a,2a-1) = (x U(a,2a) - U(a-1,2(a-1))) / (2a - 2)
897 */
898 gsl_sf_result lnU00, lnU12;
899 gsl_sf_result U00, U12;
900 hyperg_lnU_beq2a(a1-1.0, x, &lnU00);
901 hyperg_lnU_beq2a(a1, x, &lnU12);
902 if(lnU00.val > lnU12.val) {
903 lm_for.val = lnU00.val;
904 lm_for.err = lnU00.err;
905 U00.val = 1.0;
906 U00.err = 0.0;
907 gsl_sf_exp_err_e(lnU12.val - lm_for.val, lnU12.err + lm_for.err, &U12);
908 }
909 else {
910 lm_for.val = lnU12.val;
911 lm_for.err = lnU12.err;
912 U12.val = 1.0;
913 U12.err = 0.0;
914 gsl_sf_exp_err_e(lnU00.val - lm_for.val, lnU00.err + lm_for.err, &U00);
915 }
916 Ua1_for_val = (x * U12.val - U00.val) / (2.0*a1 - 2.0);
917 Ua1_for_err = (fabs(x)*U12.err + U00.err) / fabs(2.0*a1 - 2.0);
918 Ua1_for_err += 2.0 * GSL_DBL_EPSILON * fabs(Ua1_for_val);
919 stat_for = GSL_SUCCESS;
920 }
921 else {
922 /* Recurse forward to determine U(a1,b) with
923 * absolute normalization.
924 */
925 gsl_sf_result r_Ua;
926 double Uam1 = 1.0; /* U(a0-1,b,x) = U(0,b,x) */
927 double Ua;
928 double Uap1;
929 int ap;
930 double lm_for_local;
931 stat_for = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_for_local); /* U(1,b,x) */
932 Ua = r_Ua.val;
933 Uam1 *= exp(-lm_for_local);
934 lm_for.val = lm_for_local;
935 lm_for.err = 0.0;
936
937 for(ap=a0; ap<a1; ap++) {
938 Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
939 Uam1 = Ua;
940 Ua = Uap1;
941 RESCALE_2(Ua,Uam1,scale_factor,scale_count_for);
942 }
943 Ua1_for_val = Ua;
944 Ua1_for_err = fabs(Ua) * fabs(r_Ua.err/r_Ua.val);
945 Ua1_for_err += 2.0 * GSL_DBL_EPSILON * (fabs(a1-a0)+1.0) * fabs(Ua1_for_val);
946 }
947
948 /* Now do the matching to produce the final result.
949 */
950 if(Ua1_bck_val == 0.0) {
951 result->val = 0.0;
952 result->err = 0.0;
953 result->e10 = 0;
954 GSL_ERROR ("error", GSL_EZERODIV);
955 }
956 else if(Ua1_for_val == 0.0) {
957 /* Should never happen. */
958 UNDERFLOW_ERROR_E10(result);
959 }
960 else {
961 double lns = (scale_count_for - scale_count_bck)*log(scale_factor);
962 double ln_for_val = log(fabs(Ua1_for_val));
963 double ln_for_err = GSL_DBL_EPSILON + fabs(Ua1_for_err/Ua1_for_val);
964 double ln_bck_val = log(fabs(Ua1_bck_val));
965 double ln_bck_err = GSL_DBL_EPSILON + fabs(Ua1_bck_err/Ua1_bck_val);
966 double lnr_val = lm_for.val + ln_for_val - ln_bck_val + lns;
967 double lnr_err = lm_for.err + ln_for_err + ln_bck_err
968 + 2.0 * GSL_DBL_EPSILON * (fabs(lm_for.val) + fabs(ln_for_val) + fabs(ln_bck_val) + fabs(lns));
969 double sgn = GSL_SIGN(Ua1_for_val) * GSL_SIGN(Ua1_bck_val);
970 int stat_e = gsl_sf_exp_err_e10_e(lnr_val, lnr_err, result);
971 result->val *= sgn;
972 return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for);
973 }
974 }
975 }
976 }
977
978
979 /* Handle b >= 1 for generic a,b values.
980 */
981 static
982 int
hyperg_U_bge1(const double a,const double b,const double x,gsl_sf_result_e10 * result)983 hyperg_U_bge1(const double a, const double b, const double x,
984 gsl_sf_result_e10 * result)
985 {
986 const double rinta = floor(a+0.5);
987 const int a_neg_integer = (a < 0.0 && fabs(a - rinta) < INT_THRESHOLD);
988
989 if(a == 0.0) {
990 result->val = 1.0;
991 result->err = 0.0;
992 result->e10 = 0;
993 return GSL_SUCCESS;
994 }
995 else if(a_neg_integer && fabs(rinta) < INT_MAX) {
996 /* U(-n,b,x) = (-1)^n n! Laguerre[n,b-1,x]
997 */
998 const int n = -(int)rinta;
999 const double sgn = (GSL_IS_ODD(n) ? -1.0 : 1.0);
1000 gsl_sf_result lnfact;
1001 gsl_sf_result L;
1002 const int stat_L = gsl_sf_laguerre_n_e(n, b-1.0, x, &L);
1003 gsl_sf_lnfact_e(n, &lnfact);
1004 {
1005 const int stat_e = gsl_sf_exp_mult_err_e10_e(lnfact.val, lnfact.err,
1006 sgn*L.val, L.err,
1007 result);
1008 return GSL_ERROR_SELECT_2(stat_e, stat_L);
1009 }
1010 }
1011 else if(ASYMP_EVAL_OK(a,b,x)) {
1012 const double ln_pre_val = -a*log(x);
1013 const double ln_pre_err = 2.0 * GSL_DBL_EPSILON * fabs(ln_pre_val);
1014 gsl_sf_result asymp;
1015 int stat_asymp = hyperg_zaU_asymp(a, b, x, &asymp);
1016 int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val, ln_pre_err,
1017 asymp.val, asymp.err,
1018 result);
1019 return GSL_ERROR_SELECT_2(stat_e, stat_asymp);
1020 }
1021 else if(fabs(a) <= 1.0) {
1022 gsl_sf_result rU;
1023 double ln_multiplier;
1024 int stat_U = hyperg_U_small_a_bgt0(a, b, x, &rU, &ln_multiplier);
1025 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);
1026 return GSL_ERROR_SELECT_2(stat_U, stat_e);
1027 }
1028 else if(SERIES_EVAL_OK(a,b,x)) {
1029 gsl_sf_result ser;
1030 const int stat_ser = hyperg_U_series(a, b, x, &ser);
1031 result->val = ser.val;
1032 result->err = ser.err;
1033 result->e10 = 0;
1034 return stat_ser;
1035 }
1036 else if(a < 0.0) {
1037 /* Recurse backward on a and then upward on b.
1038 */
1039 const double scale_factor = GSL_SQRT_DBL_MAX;
1040 const double a0 = a - floor(a) - 1.0;
1041 const double b0 = b - floor(b) + 1.0;
1042 int scale_count = 0;
1043 double lm_0, lm_1;
1044 double lm_max;
1045 gsl_sf_result r_Uap1;
1046 gsl_sf_result r_Ua;
1047 int stat_0 = hyperg_U_small_a_bgt0(a0+1.0, b0, x, &r_Uap1, &lm_0);
1048 int stat_1 = hyperg_U_small_a_bgt0(a0, b0, x, &r_Ua, &lm_1);
1049 int stat_e;
1050 double Uap1 = r_Uap1.val;
1051 double Ua = r_Ua.val;
1052 double Uam1;
1053 double ap;
1054 lm_max = GSL_MAX(lm_0, lm_1);
1055 Uap1 *= exp(lm_0-lm_max);
1056 Ua *= exp(lm_1-lm_max);
1057
1058 /* Downward recursion on a.
1059 */
1060 for(ap=a0; ap>a+0.1; ap -= 1.0) {
1061 Uam1 = ap*(b0-ap-1.0)*Uap1 + (x+2.0*ap-b0)*Ua;
1062 Uap1 = Ua;
1063 Ua = Uam1;
1064 RESCALE_2(Ua,Uap1,scale_factor,scale_count);
1065 }
1066
1067 if(b < 2.0) {
1068 /* b == b0, so no recursion necessary
1069 */
1070 const double lnscale = log(scale_factor);
1071 gsl_sf_result lnm;
1072 gsl_sf_result y;
1073 lnm.val = lm_max + scale_count * lnscale;
1074 lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + scale_count * fabs(lnscale));
1075 y.val = Ua;
1076 y.err = fabs(r_Uap1.err/r_Uap1.val) * fabs(Ua);
1077 y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
1078 y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua);
1079 y.err *= fabs(lm_0-lm_max) + 1.0;
1080 y.err *= fabs(lm_1-lm_max) + 1.0;
1081 stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
1082 }
1083 else {
1084 /* Upward recursion on b.
1085 */
1086 const double err_mult = fabs(b-b0) + fabs(a-a0) + 1.0;
1087 const double lnscale = log(scale_factor);
1088 gsl_sf_result lnm;
1089 gsl_sf_result y;
1090
1091 double Ubm1 = Ua; /* U(a,b0) */
1092 double Ub = (a*(b0-a-1.0)*Uap1 + (a+x)*Ua)/x; /* U(a,b0+1) */
1093 double Ubp1;
1094 double bp;
1095 for(bp=b0+1.0; bp<b-0.1; bp += 1.0) {
1096 Ubp1 = ((1.0+a-bp)*Ubm1 + (bp+x-1.0)*Ub)/x;
1097 Ubm1 = Ub;
1098 Ub = Ubp1;
1099 RESCALE_2(Ub,Ubm1,scale_factor,scale_count);
1100 }
1101
1102 lnm.val = lm_max + scale_count * lnscale;
1103 lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale));
1104 y.val = Ub;
1105 y.err = 2.0 * err_mult * fabs(r_Uap1.err/r_Uap1.val) * fabs(Ub);
1106 y.err += 2.0 * err_mult * fabs(r_Ua.err/r_Ua.val) * fabs(Ub);
1107 y.err += 2.0 * GSL_DBL_EPSILON * err_mult * fabs(Ub);
1108 y.err *= fabs(lm_0-lm_max) + 1.0;
1109 y.err *= fabs(lm_1-lm_max) + 1.0;
1110 stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
1111 }
1112 return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1);
1113 }
1114 else if(b >= 2*a + x) {
1115 /* Recurse forward from a near zero.
1116 * Note that we cannot cross the singularity at
1117 * the line b=a+1, because the only way we could
1118 * be in that little wedge is if a < 1. But we
1119 * have already dealt with the small a case.
1120 */
1121 int scale_count = 0;
1122 const double a0 = a - floor(a);
1123 const double scale_factor = GSL_SQRT_DBL_MAX;
1124 double lnscale;
1125 double lm_0, lm_1, lm_max;
1126 gsl_sf_result r_Uam1;
1127 gsl_sf_result r_Ua;
1128 int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0);
1129 int stat_1 = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_1);
1130 int stat_e;
1131 gsl_sf_result lnm;
1132 gsl_sf_result y;
1133 double Uam1 = r_Uam1.val;
1134 double Ua = r_Ua.val;
1135 double Uap1;
1136 double ap;
1137 lm_max = GSL_MAX(lm_0, lm_1);
1138 Uam1 *= exp(lm_0-lm_max);
1139 Ua *= exp(lm_1-lm_max);
1140
1141 for(ap=a0; ap<a-0.1; ap += 1.0) {
1142 Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
1143 Uam1 = Ua;
1144 Ua = Uap1;
1145 RESCALE_2(Ua,Uam1,scale_factor,scale_count);
1146 }
1147
1148 lnscale = log(scale_factor);
1149 lnm.val = lm_max + scale_count * lnscale;
1150 lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_max) + fabs(scale_count * lnscale));
1151 y.val = Ua;
1152 y.err = fabs(r_Uam1.err/r_Uam1.val) * fabs(Ua);
1153 y.err += fabs(r_Ua.err/r_Ua.val) * fabs(Ua);
1154 y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + 1.0) * fabs(Ua);
1155 y.err *= fabs(lm_0-lm_max) + 1.0;
1156 y.err *= fabs(lm_1-lm_max) + 1.0;
1157 stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
1158 return GSL_ERROR_SELECT_3(stat_e, stat_0, stat_1);
1159 }
1160 else {
1161 if(b <= x) {
1162 /* Recurse backward to a near zero.
1163 */
1164 const double a0 = a - floor(a);
1165 const double scale_factor = GSL_SQRT_DBL_MAX;
1166 int scale_count = 0;
1167 gsl_sf_result lnm;
1168 gsl_sf_result y;
1169 double lnscale;
1170 double lm_0;
1171 double Uap1;
1172 double Ua;
1173 double Uam1;
1174 gsl_sf_result U0;
1175 double ap;
1176 double ru;
1177 double r;
1178 int CF1_count;
1179 int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
1180 int stat_U0;
1181 int stat_e;
1182 r = ru/a;
1183 Ua = GSL_SQRT_DBL_MIN;
1184 Uap1 = r * Ua;
1185 for(ap=a; ap>a0+0.1; ap -= 1.0) {
1186 Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
1187 Uap1 = Ua;
1188 Ua = Uam1;
1189 RESCALE_2(Ua,Uap1,scale_factor,scale_count);
1190 }
1191
1192 stat_U0 = hyperg_U_small_a_bgt0(a0, b, x, &U0, &lm_0);
1193
1194 lnscale = log(scale_factor);
1195 lnm.val = lm_0 - scale_count * lnscale;
1196 lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_0) + fabs(scale_count * lnscale));
1197 y.val = GSL_SQRT_DBL_MIN*(U0.val/Ua);
1198 y.err = GSL_SQRT_DBL_MIN*(U0.err/fabs(Ua));
1199 y.err += 2.0 * GSL_DBL_EPSILON * (fabs(a0-a) + CF1_count + 1.0) * fabs(y.val);
1200 stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
1201 return GSL_ERROR_SELECT_3(stat_e, stat_U0, stat_CF1);
1202 }
1203 else {
1204 /* Recurse backward to near the b=2a+x line, then
1205 * forward from a near zero to get the normalization.
1206 */
1207 int scale_count_for = 0;
1208 int scale_count_bck = 0;
1209 const double scale_factor = GSL_SQRT_DBL_MAX;
1210 const double eps = a - floor(a);
1211 const double a0 = ( eps == 0.0 ? 1.0 : eps );
1212 const double a1 = a0 + ceil(0.5*(b-x) - a0);
1213 gsl_sf_result lnm;
1214 gsl_sf_result y;
1215 double lm_for;
1216 double lnscale;
1217 double Ua1_bck;
1218 double Ua1_for;
1219 int stat_for;
1220 int stat_bck;
1221 int stat_e;
1222 int CF1_count;
1223
1224 {
1225 /* Recurse back to determine U(a1,b), sans normalization.
1226 */
1227 double Uap1;
1228 double Ua;
1229 double Uam1;
1230 double ap;
1231 double ru;
1232 double r;
1233 int stat_CF1 = hyperg_U_CF1(a, b, 0, x, &ru, &CF1_count);
1234 r = ru/a;
1235 Ua = GSL_SQRT_DBL_MIN;
1236 Uap1 = r * Ua;
1237 for(ap=a; ap>a1+0.1; ap -= 1.0) {
1238 Uam1 = -((b-2.0*ap-x)*Ua + ap*(1.0+ap-b)*Uap1);
1239 Uap1 = Ua;
1240 Ua = Uam1;
1241 RESCALE_2(Ua,Uap1,scale_factor,scale_count_bck);
1242 }
1243 Ua1_bck = Ua;
1244 stat_bck = stat_CF1;
1245 }
1246 {
1247 /* Recurse forward to determine U(a1,b) with
1248 * absolute normalization.
1249 */
1250 gsl_sf_result r_Uam1;
1251 gsl_sf_result r_Ua;
1252 double lm_0, lm_1;
1253 int stat_0 = hyperg_U_small_a_bgt0(a0-1.0, b, x, &r_Uam1, &lm_0);
1254 int stat_1 = hyperg_U_small_a_bgt0(a0, b, x, &r_Ua, &lm_1);
1255 double Uam1 = r_Uam1.val;
1256 double Ua = r_Ua.val;
1257 double Uap1;
1258 double ap;
1259
1260 lm_for = GSL_MAX(lm_0, lm_1);
1261 Uam1 *= exp(lm_0 - lm_for);
1262 Ua *= exp(lm_1 - lm_for);
1263
1264 for(ap=a0; ap<a1-0.1; ap += 1.0) {
1265 Uap1 = -(Uam1 + (b-2.0*ap-x)*Ua)/(ap*(1.0+ap-b));
1266 Uam1 = Ua;
1267 Ua = Uap1;
1268 RESCALE_2(Ua,Uam1,scale_factor,scale_count_for);
1269 }
1270 Ua1_for = Ua;
1271 stat_for = GSL_ERROR_SELECT_2(stat_0, stat_1);
1272 }
1273
1274 lnscale = log(scale_factor);
1275 lnm.val = lm_for + (scale_count_for - scale_count_bck)*lnscale;
1276 lnm.err = 2.0 * GSL_DBL_EPSILON * (fabs(lm_for) + fabs(scale_count_for - scale_count_bck)*fabs(lnscale));
1277 y.val = GSL_SQRT_DBL_MIN*Ua1_for/Ua1_bck;
1278 y.err = 2.0 * GSL_DBL_EPSILON * (fabs(a-a0) + CF1_count + 1.0) * fabs(y.val);
1279 stat_e = gsl_sf_exp_mult_err_e10_e(lnm.val, lnm.err, y.val, y.err, result);
1280 return GSL_ERROR_SELECT_3(stat_e, stat_bck, stat_for);
1281 }
1282 }
1283 }
1284
1285
1286 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
1287
1288
1289 int
gsl_sf_hyperg_U_int_e10_e(const int a,const int b,const double x,gsl_sf_result_e10 * result)1290 gsl_sf_hyperg_U_int_e10_e(const int a, const int b, const double x,
1291 gsl_sf_result_e10 * result)
1292 {
1293 /* CHECK_POINTER(result) */
1294
1295 if(x <= 0.0) {
1296 DOMAIN_ERROR_E10(result);
1297 }
1298 else {
1299 if(b >= 1) {
1300 return hyperg_U_int_bge1(a, b, x, result);
1301 }
1302 else {
1303 /* Use the reflection formula
1304 * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x)
1305 */
1306 gsl_sf_result_e10 U;
1307 double ln_x = log(x);
1308 int ap = 1 + a - b;
1309 int bp = 2 - b;
1310 int stat_e;
1311 int stat_U = hyperg_U_int_bge1(ap, bp, x, &U);
1312 double ln_pre_val = (1.0-b)*ln_x;
1313 double ln_pre_err = 2.0 * GSL_DBL_EPSILON * (fabs(b)+1.0) * fabs(ln_x);
1314 ln_pre_err += 2.0 * GSL_DBL_EPSILON * fabs(1.0-b); /* error in log(x) */
1315 stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
1316 U.val, U.err,
1317 result);
1318 return GSL_ERROR_SELECT_2(stat_e, stat_U);
1319 }
1320 }
1321 }
1322
1323
1324 int
gsl_sf_hyperg_U_e10_e(const double a,const double b,const double x,gsl_sf_result_e10 * result)1325 gsl_sf_hyperg_U_e10_e(const double a, const double b, const double x,
1326 gsl_sf_result_e10 * result)
1327 {
1328 const double rinta = floor(a + 0.5);
1329 const double rintb = floor(b + 0.5);
1330 const int a_integer = ( fabs(a - rinta) < INT_THRESHOLD );
1331 const int b_integer = ( fabs(b - rintb) < INT_THRESHOLD );
1332
1333 /* CHECK_POINTER(result) */
1334
1335 if(x <= 0.0) {
1336 DOMAIN_ERROR_E10(result);
1337 }
1338 else if(a == 0.0) {
1339 result->val = 1.0;
1340 result->err = 0.0;
1341 result->e10 = 0;
1342 return GSL_SUCCESS;
1343 }
1344 else if(a_integer && b_integer) {
1345 return gsl_sf_hyperg_U_int_e10_e(rinta, rintb, x, result);
1346 }
1347 else {
1348 if(b >= 1.0) {
1349 /* Use b >= 1 function.
1350 */
1351 return hyperg_U_bge1(a, b, x, result);
1352 }
1353 else {
1354 /* Use the reflection formula
1355 * U(a,b,x) = x^(1-b) U(1+a-b,2-b,x)
1356 */
1357 const double lnx = log(x);
1358 const double ln_pre_val = (1.0-b)*lnx;
1359 const double ln_pre_err = fabs(lnx) * 2.0 * GSL_DBL_EPSILON * (1.0 + fabs(b));
1360 const double ap = 1.0 + a - b;
1361 const double bp = 2.0 - b;
1362 gsl_sf_result_e10 U;
1363 int stat_U = hyperg_U_bge1(ap, bp, x, &U);
1364 int stat_e = gsl_sf_exp_mult_err_e10_e(ln_pre_val + U.e10*M_LN10, ln_pre_err,
1365 U.val, U.err,
1366 result);
1367 return GSL_ERROR_SELECT_2(stat_e, stat_U);
1368 }
1369 }
1370 }
1371
1372
1373 int
gsl_sf_hyperg_U_int_e(const int a,const int b,const double x,gsl_sf_result * result)1374 gsl_sf_hyperg_U_int_e(const int a, const int b, const double x, gsl_sf_result * result)
1375 {
1376 gsl_sf_result_e10 re;
1377 int stat_U = gsl_sf_hyperg_U_int_e10_e(a, b, x, &re);
1378 int stat_c = gsl_sf_result_smash_e(&re, result);
1379 return GSL_ERROR_SELECT_2(stat_c, stat_U);
1380 }
1381
1382
1383 int
gsl_sf_hyperg_U_e(const double a,const double b,const double x,gsl_sf_result * result)1384 gsl_sf_hyperg_U_e(const double a, const double b, const double x, gsl_sf_result * result)
1385 {
1386 gsl_sf_result_e10 re;
1387 int stat_U = gsl_sf_hyperg_U_e10_e(a, b, x, &re);
1388 int stat_c = gsl_sf_result_smash_e(&re, result);
1389 return GSL_ERROR_SELECT_2(stat_c, stat_U);
1390 }
1391
1392
1393 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
1394
1395 #include "gsl_specfunc__eval.h"
1396
gsl_sf_hyperg_U_int(const int a,const int b,const double x)1397 double gsl_sf_hyperg_U_int(const int a, const int b, const double x)
1398 {
1399 EVAL_RESULT(gsl_sf_hyperg_U_int_e(a, b, x, &result));
1400 }
1401
gsl_sf_hyperg_U(const double a,const double b,const double x)1402 double gsl_sf_hyperg_U(const double a, const double b, const double x)
1403 {
1404 EVAL_RESULT(gsl_sf_hyperg_U_e(a, b, x, &result));
1405 }
1406