1 /* specfunc/hyperg_2F1.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2004 Gerard Jungman
4 * Copyright (C) 2009 Brian Gough
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
19 */
20
21 /* Author: G. Jungman */
22
23 #include <config.h>
24 #include <gsl/gsl_math.h>
25 #include <gsl/gsl_errno.h>
26 #include <gsl/gsl_sf_exp.h>
27 #include <gsl/gsl_sf_pow_int.h>
28 #include <gsl/gsl_sf_gamma.h>
29 #include <gsl/gsl_sf_psi.h>
30 #include <gsl/gsl_sf_hyperg.h>
31
32 #include "error.h"
33
34 #define locEPS (1000.0*GSL_DBL_EPSILON)
35
36
37 /* Assumes c != negative integer.
38 */
39 static int
hyperg_2F1_series(const double a,const double b,const double c,const double x,gsl_sf_result * result)40 hyperg_2F1_series(const double a, const double b, const double c,
41 const double x,
42 gsl_sf_result * result
43 )
44 {
45 double sum_pos = 1.0;
46 double sum_neg = 0.0;
47 double del_pos = 1.0;
48 double del_neg = 0.0;
49 double del = 1.0;
50 double del_prev;
51 double k = 0.0;
52 int i = 0;
53
54 if(fabs(c) < GSL_DBL_EPSILON) {
55 result->val = 0.0; /* FIXME: ?? */
56 result->err = 1.0;
57 GSL_ERROR ("error", GSL_EDOM);
58 }
59
60 do {
61 if(++i > 30000) {
62 result->val = sum_pos - sum_neg;
63 result->err = del_pos + del_neg;
64 result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
65 result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
66 GSL_ERROR ("error", GSL_EMAXITER);
67 }
68 del_prev = del;
69 del *= (a+k)*(b+k) * x / ((c+k) * (k+1.0)); /* Gauss series */
70
71 if(del > 0.0) {
72 del_pos = del;
73 sum_pos += del;
74 }
75 else if(del == 0.0) {
76 /* Exact termination (a or b was a negative integer).
77 */
78 del_pos = 0.0;
79 del_neg = 0.0;
80 break;
81 }
82 else {
83 del_neg = -del;
84 sum_neg -= del;
85 }
86
87 /*
88 * This stopping criteria is taken from the thesis
89 * "Computation of Hypergeometic Functions" by J. Pearson, pg. 31
90 * (http://people.maths.ox.ac.uk/porterm/research/pearson_final.pdf)
91 * and fixes bug #45926
92 */
93 if (fabs(del_prev / (sum_pos - sum_neg)) < GSL_DBL_EPSILON &&
94 fabs(del / (sum_pos - sum_neg)) < GSL_DBL_EPSILON)
95 break;
96
97 k += 1.0;
98 } while(fabs((del_pos + del_neg)/(sum_pos-sum_neg)) > GSL_DBL_EPSILON);
99
100 result->val = sum_pos - sum_neg;
101 result->err = del_pos + del_neg;
102 result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
103 result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);
104
105 return GSL_SUCCESS;
106 }
107
108
109 /* a = aR + i aI, b = aR - i aI */
110 static
111 int
hyperg_2F1_conj_series(const double aR,const double aI,const double c,double x,gsl_sf_result * result)112 hyperg_2F1_conj_series(const double aR, const double aI, const double c,
113 double x,
114 gsl_sf_result * result)
115 {
116 if(c == 0.0) {
117 result->val = 0.0; /* FIXME: should be Inf */
118 result->err = 0.0;
119 GSL_ERROR ("error", GSL_EDOM);
120 }
121 else {
122 double sum_pos = 1.0;
123 double sum_neg = 0.0;
124 double del_pos = 1.0;
125 double del_neg = 0.0;
126 double del = 1.0;
127 double k = 0.0;
128 do {
129 del *= ((aR+k)*(aR+k) + aI*aI)/((k+1.0)*(c+k)) * x;
130
131 if(del >= 0.0) {
132 del_pos = del;
133 sum_pos += del;
134 }
135 else {
136 del_neg = -del;
137 sum_neg -= del;
138 }
139
140 if(k > 30000) {
141 result->val = sum_pos - sum_neg;
142 result->err = del_pos + del_neg;
143 result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
144 result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k)+1.0) * fabs(result->val);
145 GSL_ERROR ("error", GSL_EMAXITER);
146 }
147
148 k += 1.0;
149 } while(fabs((del_pos + del_neg)/(sum_pos - sum_neg)) > GSL_DBL_EPSILON);
150
151 result->val = sum_pos - sum_neg;
152 result->err = del_pos + del_neg;
153 result->err += 2.0 * GSL_DBL_EPSILON * (sum_pos + sum_neg);
154 result->err += 2.0 * GSL_DBL_EPSILON * (2.0*sqrt(k) + 1.0) * fabs(result->val);
155
156 return GSL_SUCCESS;
157 }
158 }
159
160
161 /* Luke's rational approximation. The most accesible
162 * discussion is in [Kolbig, CPC 23, 51 (1981)].
163 * The convergence is supposedly guaranteed for x < 0.
164 * You have to read Luke's books to see this and other
165 * results. Unfortunately, the stability is not so
166 * clear to me, although it seems very efficient when
167 * it works.
168 */
169 static
170 int
hyperg_2F1_luke(const double a,const double b,const double c,const double xin,gsl_sf_result * result)171 hyperg_2F1_luke(const double a, const double b, const double c,
172 const double xin,
173 gsl_sf_result * result)
174 {
175 int stat_iter;
176 const double RECUR_BIG = 1.0e+50;
177 const int nmax = 20000;
178 int n = 3;
179 const double x = -xin;
180 const double x3 = x*x*x;
181 const double t0 = a*b/c;
182 const double t1 = (a+1.0)*(b+1.0)/(2.0*c);
183 const double t2 = (a+2.0)*(b+2.0)/(2.0*(c+1.0));
184 double F = 1.0;
185 double prec;
186
187 double Bnm3 = 1.0; /* B0 */
188 double Bnm2 = 1.0 + t1 * x; /* B1 */
189 double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x); /* B2 */
190
191 double Anm3 = 1.0; /* A0 */
192 double Anm2 = Bnm2 - t0 * x; /* A1 */
193 double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x; /* A2 */
194
195 while(1) {
196 double npam1 = n + a - 1;
197 double npbm1 = n + b - 1;
198 double npcm1 = n + c - 1;
199 double npam2 = n + a - 2;
200 double npbm2 = n + b - 2;
201 double npcm2 = n + c - 2;
202 double tnm1 = 2*n - 1;
203 double tnm3 = 2*n - 3;
204 double tnm5 = 2*n - 5;
205 double n2 = n*n;
206 double F1 = (3.0*n2 + (a+b-6)*n + 2 - a*b - 2*(a+b)) / (2*tnm3*npcm1);
207 double F2 = -(3.0*n2 - (a+b+6)*n + 2 - a*b)*npam1*npbm1/(4*tnm1*tnm3*npcm2*npcm1);
208 double F3 = (npam2*npam1*npbm2*npbm1*(n-a-2)*(n-b-2)) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
209 double E = -npam1*npbm1*(n-c-1) / (2*tnm3*npcm2*npcm1);
210
211 double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
212 double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
213 double r = An/Bn;
214
215 prec = fabs((F - r)/F);
216 F = r;
217
218 if(prec < GSL_DBL_EPSILON || n > nmax) break;
219
220 if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
221 An /= RECUR_BIG;
222 Bn /= RECUR_BIG;
223 Anm1 /= RECUR_BIG;
224 Bnm1 /= RECUR_BIG;
225 Anm2 /= RECUR_BIG;
226 Bnm2 /= RECUR_BIG;
227 Anm3 /= RECUR_BIG;
228 Bnm3 /= RECUR_BIG;
229 }
230 else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
231 An *= RECUR_BIG;
232 Bn *= RECUR_BIG;
233 Anm1 *= RECUR_BIG;
234 Bnm1 *= RECUR_BIG;
235 Anm2 *= RECUR_BIG;
236 Bnm2 *= RECUR_BIG;
237 Anm3 *= RECUR_BIG;
238 Bnm3 *= RECUR_BIG;
239 }
240
241 n++;
242 Bnm3 = Bnm2;
243 Bnm2 = Bnm1;
244 Bnm1 = Bn;
245 Anm3 = Anm2;
246 Anm2 = Anm1;
247 Anm1 = An;
248 }
249
250 result->val = F;
251 result->err = 2.0 * fabs(prec * F);
252 result->err += 2.0 * GSL_DBL_EPSILON * (n+1.0) * fabs(F);
253
254 /* FIXME: just a hack: there's a lot of shit going on here */
255 result->err *= 8.0 * (fabs(a) + fabs(b) + 1.0);
256
257 stat_iter = (n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );
258
259 return stat_iter;
260 }
261
262
263 /* Luke's rational approximation for the
264 * case a = aR + i aI, b = aR - i aI.
265 */
266 static
267 int
hyperg_2F1_conj_luke(const double aR,const double aI,const double c,const double xin,gsl_sf_result * result)268 hyperg_2F1_conj_luke(const double aR, const double aI, const double c,
269 const double xin,
270 gsl_sf_result * result)
271 {
272 int stat_iter;
273 const double RECUR_BIG = 1.0e+50;
274 const int nmax = 10000;
275 int n = 3;
276 const double x = -xin;
277 const double x3 = x*x*x;
278 const double atimesb = aR*aR + aI*aI;
279 const double apb = 2.0*aR;
280 const double t0 = atimesb/c;
281 const double t1 = (atimesb + apb + 1.0)/(2.0*c);
282 const double t2 = (atimesb + 2.0*apb + 4.0)/(2.0*(c+1.0));
283 double F = 1.0;
284 double prec;
285
286 double Bnm3 = 1.0; /* B0 */
287 double Bnm2 = 1.0 + t1 * x; /* B1 */
288 double Bnm1 = 1.0 + t2 * x * (1.0 + t1/3.0 * x); /* B2 */
289
290 double Anm3 = 1.0; /* A0 */
291 double Anm2 = Bnm2 - t0 * x; /* A1 */
292 double Anm1 = Bnm1 - t0*(1.0 + t2*x)*x + t0 * t1 * (c/(c+1.0)) * x*x; /* A2 */
293
294 while(1) {
295 double nm1 = n - 1;
296 double nm2 = n - 2;
297 double npam1_npbm1 = atimesb + nm1*apb + nm1*nm1;
298 double npam2_npbm2 = atimesb + nm2*apb + nm2*nm2;
299 double npcm1 = nm1 + c;
300 double npcm2 = nm2 + c;
301 double tnm1 = 2*n - 1;
302 double tnm3 = 2*n - 3;
303 double tnm5 = 2*n - 5;
304 double n2 = n*n;
305 double F1 = (3.0*n2 + (apb-6)*n + 2 - atimesb - 2*apb) / (2*tnm3*npcm1);
306 double F2 = -(3.0*n2 - (apb+6)*n + 2 - atimesb)*npam1_npbm1/(4*tnm1*tnm3*npcm2*npcm1);
307 double F3 = (npam2_npbm2*npam1_npbm1*(nm2*nm2 - nm2*apb + atimesb)) / (8*tnm3*tnm3*tnm5*(n+c-3)*npcm2*npcm1);
308 double E = -npam1_npbm1*(n-c-1) / (2*tnm3*npcm2*npcm1);
309
310 double An = (1.0+F1*x)*Anm1 + (E + F2*x)*x*Anm2 + F3*x3*Anm3;
311 double Bn = (1.0+F1*x)*Bnm1 + (E + F2*x)*x*Bnm2 + F3*x3*Bnm3;
312 double r = An/Bn;
313
314 prec = fabs(F - r)/fabs(F);
315 F = r;
316
317 if(prec < GSL_DBL_EPSILON || n > nmax) break;
318
319 if(fabs(An) > RECUR_BIG || fabs(Bn) > RECUR_BIG) {
320 An /= RECUR_BIG;
321 Bn /= RECUR_BIG;
322 Anm1 /= RECUR_BIG;
323 Bnm1 /= RECUR_BIG;
324 Anm2 /= RECUR_BIG;
325 Bnm2 /= RECUR_BIG;
326 Anm3 /= RECUR_BIG;
327 Bnm3 /= RECUR_BIG;
328 }
329 else if(fabs(An) < 1.0/RECUR_BIG || fabs(Bn) < 1.0/RECUR_BIG) {
330 An *= RECUR_BIG;
331 Bn *= RECUR_BIG;
332 Anm1 *= RECUR_BIG;
333 Bnm1 *= RECUR_BIG;
334 Anm2 *= RECUR_BIG;
335 Bnm2 *= RECUR_BIG;
336 Anm3 *= RECUR_BIG;
337 Bnm3 *= RECUR_BIG;
338 }
339
340 n++;
341 Bnm3 = Bnm2;
342 Bnm2 = Bnm1;
343 Bnm1 = Bn;
344 Anm3 = Anm2;
345 Anm2 = Anm1;
346 Anm1 = An;
347 }
348
349 result->val = F;
350 result->err = 2.0 * fabs(prec * F);
351 result->err += 2.0 * GSL_DBL_EPSILON * (n+1.0) * fabs(F);
352
353 /* FIXME: see above */
354 result->err *= 8.0 * (fabs(aR) + fabs(aI) + 1.0);
355
356 stat_iter = (n >= nmax ? GSL_EMAXITER : GSL_SUCCESS );
357
358 return stat_iter;
359 }
360
361
362 /* Do the reflection described in [Moshier, p. 334].
363 * Assumes a,b,c != neg integer.
364 */
365 static
366 int
hyperg_2F1_reflect(const double a,const double b,const double c,const double x,gsl_sf_result * result)367 hyperg_2F1_reflect(const double a, const double b, const double c,
368 const double x, gsl_sf_result * result)
369 {
370 const double d = c - a - b;
371 const int intd = floor(d+0.5);
372 const int d_integer = ( fabs(d - intd) < locEPS );
373
374 if(d_integer) {
375 const double ln_omx = log(1.0 - x);
376 const double ad = fabs(d);
377 int stat_F2 = GSL_SUCCESS;
378 double sgn_2;
379 gsl_sf_result F1;
380 gsl_sf_result F2;
381 double d1, d2;
382 gsl_sf_result lng_c;
383 gsl_sf_result lng_ad2;
384 gsl_sf_result lng_bd2;
385 int stat_c;
386 int stat_ad2;
387 int stat_bd2;
388
389 if(d >= 0.0) {
390 d1 = d;
391 d2 = 0.0;
392 }
393 else {
394 d1 = 0.0;
395 d2 = d;
396 }
397
398 stat_ad2 = gsl_sf_lngamma_e(a+d2, &lng_ad2);
399 stat_bd2 = gsl_sf_lngamma_e(b+d2, &lng_bd2);
400 stat_c = gsl_sf_lngamma_e(c, &lng_c);
401
402 /* Evaluate F1.
403 */
404 if(ad < GSL_DBL_EPSILON) {
405 /* d = 0 */
406 F1.val = 0.0;
407 F1.err = 0.0;
408 }
409 else {
410 gsl_sf_result lng_ad;
411 gsl_sf_result lng_ad1;
412 gsl_sf_result lng_bd1;
413 int stat_ad = gsl_sf_lngamma_e(ad, &lng_ad);
414 int stat_ad1 = gsl_sf_lngamma_e(a+d1, &lng_ad1);
415 int stat_bd1 = gsl_sf_lngamma_e(b+d1, &lng_bd1);
416
417 if(stat_ad1 == GSL_SUCCESS && stat_bd1 == GSL_SUCCESS && stat_ad == GSL_SUCCESS) {
418 /* Gamma functions in the denominator are ok.
419 * Proceed with evaluation.
420 */
421 int i;
422 double sum1 = 1.0;
423 double term = 1.0;
424 double ln_pre1_val = lng_ad.val + lng_c.val + d2*ln_omx - lng_ad1.val - lng_bd1.val;
425 double ln_pre1_err = lng_ad.err + lng_c.err + lng_ad1.err + lng_bd1.err + GSL_DBL_EPSILON * fabs(ln_pre1_val);
426 int stat_e;
427
428 /* Do F1 sum.
429 */
430 for(i=1; i<ad; i++) {
431 int j = i-1;
432 term *= (a + d2 + j) * (b + d2 + j) / (1.0 + d2 + j) / i * (1.0-x);
433 sum1 += term;
434 }
435
436 stat_e = gsl_sf_exp_mult_err_e(ln_pre1_val, ln_pre1_err,
437 sum1, GSL_DBL_EPSILON*fabs(sum1),
438 &F1);
439 if(stat_e == GSL_EOVRFLW) {
440 OVERFLOW_ERROR(result);
441 }
442 }
443 else {
444 /* Gamma functions in the denominator were not ok.
445 * So the F1 term is zero.
446 */
447 F1.val = 0.0;
448 F1.err = 0.0;
449 }
450 } /* end F1 evaluation */
451
452
453 /* Evaluate F2.
454 */
455 if(stat_ad2 == GSL_SUCCESS && stat_bd2 == GSL_SUCCESS) {
456 /* Gamma functions in the denominator are ok.
457 * Proceed with evaluation.
458 */
459 const int maxiter = 2000;
460 double psi_1 = -M_EULER;
461 gsl_sf_result psi_1pd;
462 gsl_sf_result psi_apd1;
463 gsl_sf_result psi_bpd1;
464 int stat_1pd = gsl_sf_psi_e(1.0 + ad, &psi_1pd);
465 int stat_apd1 = gsl_sf_psi_e(a + d1, &psi_apd1);
466 int stat_bpd1 = gsl_sf_psi_e(b + d1, &psi_bpd1);
467 int stat_dall = GSL_ERROR_SELECT_3(stat_1pd, stat_apd1, stat_bpd1);
468
469 double psi_val = psi_1 + psi_1pd.val - psi_apd1.val - psi_bpd1.val - ln_omx;
470 double psi_err = psi_1pd.err + psi_apd1.err + psi_bpd1.err + GSL_DBL_EPSILON*fabs(psi_val);
471 double fact = 1.0;
472 double sum2_val = psi_val;
473 double sum2_err = psi_err;
474 double ln_pre2_val = lng_c.val + d1*ln_omx - lng_ad2.val - lng_bd2.val;
475 double ln_pre2_err = lng_c.err + lng_ad2.err + lng_bd2.err + GSL_DBL_EPSILON * fabs(ln_pre2_val);
476 int stat_e;
477
478 int j;
479
480 /* Do F2 sum.
481 */
482 for(j=1; j<maxiter; j++) {
483 /* values for psi functions use recurrence; Abramowitz+Stegun 6.3.5 */
484 double term1 = 1.0/(double)j + 1.0/(ad+j);
485 double term2 = 1.0/(a+d1+j-1.0) + 1.0/(b+d1+j-1.0);
486 double delta = 0.0;
487 psi_val += term1 - term2;
488 psi_err += GSL_DBL_EPSILON * (fabs(term1) + fabs(term2));
489 fact *= (a+d1+j-1.0)*(b+d1+j-1.0)/((ad+j)*j) * (1.0-x);
490 delta = fact * psi_val;
491 sum2_val += delta;
492 sum2_err += fabs(fact * psi_err) + GSL_DBL_EPSILON*fabs(delta);
493 if(fabs(delta) < GSL_DBL_EPSILON * fabs(sum2_val)) break;
494 }
495
496 if(j == maxiter) stat_F2 = GSL_EMAXITER;
497
498 if(sum2_val == 0.0) {
499 F2.val = 0.0;
500 F2.err = 0.0;
501 }
502 else {
503 stat_e = gsl_sf_exp_mult_err_e(ln_pre2_val, ln_pre2_err,
504 sum2_val, sum2_err,
505 &F2);
506 if(stat_e == GSL_EOVRFLW) {
507 result->val = 0.0;
508 result->err = 0.0;
509 GSL_ERROR ("error", GSL_EOVRFLW);
510 }
511 }
512 stat_F2 = GSL_ERROR_SELECT_2(stat_F2, stat_dall);
513 }
514 else {
515 /* Gamma functions in the denominator not ok.
516 * So the F2 term is zero.
517 */
518 F2.val = 0.0;
519 F2.err = 0.0;
520 } /* end F2 evaluation */
521
522 sgn_2 = ( GSL_IS_ODD(intd) ? -1.0 : 1.0 );
523 result->val = F1.val + sgn_2 * F2.val;
524 result->err = F1.err + F2. err;
525 result->err += 2.0 * GSL_DBL_EPSILON * (fabs(F1.val) + fabs(F2.val));
526 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
527 return stat_F2;
528 }
529 else {
530 /* d not an integer */
531
532 gsl_sf_result pre1, pre2;
533 double sgn1, sgn2;
534 gsl_sf_result F1, F2;
535 int status_F1, status_F2;
536
537 /* These gamma functions appear in the denominator, so we
538 * catch their harmless domain errors and set the terms to zero.
539 */
540 gsl_sf_result ln_g1ca, ln_g1cb, ln_g2a, ln_g2b;
541 double sgn_g1ca, sgn_g1cb, sgn_g2a, sgn_g2b;
542 int stat_1ca = gsl_sf_lngamma_sgn_e(c-a, &ln_g1ca, &sgn_g1ca);
543 int stat_1cb = gsl_sf_lngamma_sgn_e(c-b, &ln_g1cb, &sgn_g1cb);
544 int stat_2a = gsl_sf_lngamma_sgn_e(a, &ln_g2a, &sgn_g2a);
545 int stat_2b = gsl_sf_lngamma_sgn_e(b, &ln_g2b, &sgn_g2b);
546 int ok1 = (stat_1ca == GSL_SUCCESS && stat_1cb == GSL_SUCCESS);
547 int ok2 = (stat_2a == GSL_SUCCESS && stat_2b == GSL_SUCCESS);
548
549 gsl_sf_result ln_gc, ln_gd, ln_gmd;
550 double sgn_gc, sgn_gd, sgn_gmd;
551 gsl_sf_lngamma_sgn_e( c, &ln_gc, &sgn_gc);
552 gsl_sf_lngamma_sgn_e( d, &ln_gd, &sgn_gd);
553 gsl_sf_lngamma_sgn_e(-d, &ln_gmd, &sgn_gmd);
554
555 sgn1 = sgn_gc * sgn_gd * sgn_g1ca * sgn_g1cb;
556 sgn2 = sgn_gc * sgn_gmd * sgn_g2a * sgn_g2b;
557
558 if(ok1 && ok2) {
559 double ln_pre1_val = ln_gc.val + ln_gd.val - ln_g1ca.val - ln_g1cb.val;
560 double ln_pre2_val = ln_gc.val + ln_gmd.val - ln_g2a.val - ln_g2b.val + d*log(1.0-x);
561 double ln_pre1_err = ln_gc.err + ln_gd.err + ln_g1ca.err + ln_g1cb.err;
562 double ln_pre2_err = ln_gc.err + ln_gmd.err + ln_g2a.err + ln_g2b.err;
563 if(ln_pre1_val < GSL_LOG_DBL_MAX && ln_pre2_val < GSL_LOG_DBL_MAX) {
564 gsl_sf_exp_err_e(ln_pre1_val, ln_pre1_err, &pre1);
565 gsl_sf_exp_err_e(ln_pre2_val, ln_pre2_err, &pre2);
566 pre1.val *= sgn1;
567 pre2.val *= sgn2;
568 }
569 else {
570 OVERFLOW_ERROR(result);
571 }
572 }
573 else if(ok1 && !ok2) {
574 double ln_pre1_val = ln_gc.val + ln_gd.val - ln_g1ca.val - ln_g1cb.val;
575 double ln_pre1_err = ln_gc.err + ln_gd.err + ln_g1ca.err + ln_g1cb.err;
576 if(ln_pre1_val < GSL_LOG_DBL_MAX) {
577 gsl_sf_exp_err_e(ln_pre1_val, ln_pre1_err, &pre1);
578 pre1.val *= sgn1;
579 pre2.val = 0.0;
580 pre2.err = 0.0;
581 }
582 else {
583 OVERFLOW_ERROR(result);
584 }
585 }
586 else if(!ok1 && ok2) {
587 double ln_pre2_val = ln_gc.val + ln_gmd.val - ln_g2a.val - ln_g2b.val + d*log(1.0-x);
588 double ln_pre2_err = ln_gc.err + ln_gmd.err + ln_g2a.err + ln_g2b.err;
589 if(ln_pre2_val < GSL_LOG_DBL_MAX) {
590 pre1.val = 0.0;
591 pre1.err = 0.0;
592 gsl_sf_exp_err_e(ln_pre2_val, ln_pre2_err, &pre2);
593 pre2.val *= sgn2;
594 }
595 else {
596 OVERFLOW_ERROR(result);
597 }
598 }
599 else {
600 pre1.val = 0.0;
601 pre2.val = 0.0;
602 UNDERFLOW_ERROR(result);
603 }
604
605 status_F1 = hyperg_2F1_series( a, b, 1.0-d, 1.0-x, &F1);
606 status_F2 = hyperg_2F1_series(c-a, c-b, 1.0+d, 1.0-x, &F2);
607
608 result->val = pre1.val*F1.val + pre2.val*F2.val;
609 result->err = fabs(pre1.val*F1.err) + fabs(pre2.val*F2.err);
610 result->err += fabs(pre1.err*F1.val) + fabs(pre2.err*F2.val);
611 result->err += 2.0 * GSL_DBL_EPSILON * (fabs(pre1.val*F1.val) + fabs(pre2.val*F2.val));
612 result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
613
614 if (status_F1)
615 return status_F1;
616
617 if (status_F2)
618 return status_F2;
619
620 return GSL_SUCCESS;
621 }
622 }
623
624
pow_omx(const double x,const double p,gsl_sf_result * result)625 static int pow_omx(const double x, const double p, gsl_sf_result * result)
626 {
627 double ln_omx;
628 double ln_result;
629 if(fabs(x) < GSL_ROOT5_DBL_EPSILON) {
630 ln_omx = -x*(1.0 + x*(1.0/2.0 + x*(1.0/3.0 + x/4.0 + x*x/5.0)));
631 }
632 else {
633 ln_omx = log(1.0-x);
634 }
635 ln_result = p * ln_omx;
636 return gsl_sf_exp_err_e(ln_result, GSL_DBL_EPSILON * fabs(ln_result), result);
637 }
638
639
640 /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
641
642 int
gsl_sf_hyperg_2F1_e(double a,double b,const double c,const double x,gsl_sf_result * result)643 gsl_sf_hyperg_2F1_e(double a, double b, const double c,
644 const double x,
645 gsl_sf_result * result)
646 {
647 const double d = c - a - b;
648 const double rinta = floor(a + 0.5);
649 const double rintb = floor(b + 0.5);
650 const double rintc = floor(c + 0.5);
651 const int a_neg_integer = ( a < 0.0 && fabs(a - rinta) < locEPS );
652 const int b_neg_integer = ( b < 0.0 && fabs(b - rintb) < locEPS );
653 const int c_neg_integer = ( c < 0.0 && fabs(c - rintc) < locEPS );
654
655 result->val = 0.0;
656 result->err = 0.0;
657
658 /* Handle x == 1.0 RJM */
659
660 if (fabs (x - 1.0) < locEPS && (c - a - b) > 0 && c != 0 && !c_neg_integer) {
661 gsl_sf_result lngamc, lngamcab, lngamca, lngamcb;
662 double lngamc_sgn, lngamca_sgn, lngamcb_sgn;
663 int status;
664 int stat1 = gsl_sf_lngamma_sgn_e (c, &lngamc, &lngamc_sgn);
665 int stat2 = gsl_sf_lngamma_e (c - a - b, &lngamcab);
666 int stat3 = gsl_sf_lngamma_sgn_e (c - a, &lngamca, &lngamca_sgn);
667 int stat4 = gsl_sf_lngamma_sgn_e (c - b, &lngamcb, &lngamcb_sgn);
668
669 if (stat1 != GSL_SUCCESS || stat2 != GSL_SUCCESS
670 || stat3 != GSL_SUCCESS || stat4 != GSL_SUCCESS)
671 {
672 DOMAIN_ERROR (result);
673 }
674
675 status =
676 gsl_sf_exp_err_e (lngamc.val + lngamcab.val - lngamca.val - lngamcb.val,
677 lngamc.err + lngamcab.err + lngamca.err + lngamcb.err,
678 result);
679
680 result->val *= lngamc_sgn / (lngamca_sgn * lngamcb_sgn);
681 return status;
682 }
683
684 if(x < -1.0 || 1.0 <= x) {
685 DOMAIN_ERROR(result);
686 }
687
688 if(c_neg_integer) {
689 /* If c is a negative integer, then either a or b must be a
690 negative integer of smaller magnitude than c to ensure
691 cancellation of the series. */
692 if(! (a_neg_integer && a > c + 0.1) && ! (b_neg_integer && b > c + 0.1)) {
693 DOMAIN_ERROR(result);
694 }
695 }
696
697 if(fabs(c-b) < locEPS || fabs(c-a) < locEPS) {
698 return pow_omx(x, d, result); /* (1-x)^(c-a-b) */
699 }
700
701 if(a >= 0.0 && b >= 0.0 && c >=0.0 && x >= 0.0 && x < 0.995) {
702 /* Series has all positive definite
703 * terms and x is not close to 1.
704 */
705 return hyperg_2F1_series(a, b, c, x, result);
706 }
707
708 if(fabs(a) < 10.0 && fabs(b) < 10.0) {
709 /* a and b are not too large, so we attempt
710 * variations on the series summation.
711 */
712 if(a_neg_integer) {
713 return hyperg_2F1_series(rinta, b, c, x, result);
714 }
715 if(b_neg_integer) {
716 return hyperg_2F1_series(a, rintb, c, x, result);
717 }
718
719 if(x < -0.25) {
720 return hyperg_2F1_luke(a, b, c, x, result);
721 }
722 else if(x < 0.5) {
723 return hyperg_2F1_series(a, b, c, x, result);
724 }
725 else {
726 if(fabs(c) > 10.0) {
727 return hyperg_2F1_series(a, b, c, x, result);
728 }
729 else {
730 return hyperg_2F1_reflect(a, b, c, x, result);
731 }
732 }
733 }
734 else {
735 /* Either a or b or both large.
736 * Introduce some new variables ap,bp so that bp is
737 * the larger in magnitude.
738 */
739 double ap, bp;
740 if(fabs(a) > fabs(b)) {
741 bp = a;
742 ap = b;
743 }
744 else {
745 bp = b;
746 ap = a;
747 }
748
749 if(x < 0.0) {
750 /* What the hell, maybe Luke will converge.
751 */
752 return hyperg_2F1_luke(a, b, c, x, result);
753 }
754
755 if(GSL_MAX_DBL(fabs(ap),1.0)*fabs(bp)*fabs(x) < 2.0*fabs(c)) {
756 /* If c is large enough or x is small enough,
757 * we can attempt the series anyway.
758 */
759 return hyperg_2F1_series(a, b, c, x, result);
760 }
761
762 if(fabs(bp*bp*x*x) < 0.001*fabs(bp) && fabs(ap) < 10.0) {
763 /* The famous but nearly worthless "large b" asymptotic.
764 */
765 int stat = gsl_sf_hyperg_1F1_e(ap, c, bp*x, result);
766 result->err = 0.001 * fabs(result->val);
767 return stat;
768 }
769
770 /* We give up. */
771 result->val = 0.0;
772 result->err = 0.0;
773 GSL_ERROR ("error", GSL_EUNIMPL);
774 }
775 }
776
777
778 int
gsl_sf_hyperg_2F1_conj_e(const double aR,const double aI,const double c,const double x,gsl_sf_result * result)779 gsl_sf_hyperg_2F1_conj_e(const double aR, const double aI, const double c,
780 const double x,
781 gsl_sf_result * result)
782 {
783 const double ax = fabs(x);
784 const double rintc = floor(c + 0.5);
785 const int c_neg_integer = ( c < 0.0 && fabs(c - rintc) < locEPS );
786
787 result->val = 0.0;
788 result->err = 0.0;
789
790 if(ax >= 1.0 || c_neg_integer || c == 0.0) {
791 DOMAIN_ERROR(result);
792 }
793
794 if( (ax < 0.25 && fabs(aR) < 20.0 && fabs(aI) < 20.0)
795 || (c > 0.0 && x > 0.0)
796 ) {
797 return hyperg_2F1_conj_series(aR, aI, c, x, result);
798 }
799 else if(fabs(aR) < 10.0 && fabs(aI) < 10.0) {
800 if(x < -0.25) {
801 return hyperg_2F1_conj_luke(aR, aI, c, x, result);
802 }
803 else {
804 return hyperg_2F1_conj_series(aR, aI, c, x, result);
805 }
806 }
807 else {
808 if(x < 0.0) {
809 /* What the hell, maybe Luke will converge.
810 */
811 return hyperg_2F1_conj_luke(aR, aI, c, x, result);
812 }
813
814 /* Give up. */
815 result->val = 0.0;
816 result->err = 0.0;
817 GSL_ERROR ("error", GSL_EUNIMPL);
818 }
819 }
820
821
822 int
gsl_sf_hyperg_2F1_renorm_e(const double a,const double b,const double c,const double x,gsl_sf_result * result)823 gsl_sf_hyperg_2F1_renorm_e(const double a, const double b, const double c,
824 const double x,
825 gsl_sf_result * result
826 )
827 {
828 const double rinta = floor(a + 0.5);
829 const double rintb = floor(b + 0.5);
830 const double rintc = floor(c + 0.5);
831 const int a_neg_integer = ( a < 0.0 && fabs(a - rinta) < locEPS );
832 const int b_neg_integer = ( b < 0.0 && fabs(b - rintb) < locEPS );
833 const int c_neg_integer = ( c < 0.0 && fabs(c - rintc) < locEPS );
834
835 if(c_neg_integer) {
836 if((a_neg_integer && a > c+0.1) || (b_neg_integer && b > c+0.1)) {
837 /* 2F1 terminates early */
838 result->val = 0.0;
839 result->err = 0.0;
840 return GSL_SUCCESS;
841 }
842 else {
843 /* 2F1 does not terminate early enough, so something survives */
844 /* [Abramowitz+Stegun, 15.1.2] */
845 gsl_sf_result g1, g2, g3, g4, g5;
846 double s1, s2, s3, s4, s5;
847 int stat = 0;
848 stat += gsl_sf_lngamma_sgn_e(a-c+1, &g1, &s1);
849 stat += gsl_sf_lngamma_sgn_e(b-c+1, &g2, &s2);
850 stat += gsl_sf_lngamma_sgn_e(a, &g3, &s3);
851 stat += gsl_sf_lngamma_sgn_e(b, &g4, &s4);
852 stat += gsl_sf_lngamma_sgn_e(-c+2, &g5, &s5);
853 if(stat != 0) {
854 DOMAIN_ERROR(result);
855 }
856 else {
857 gsl_sf_result F;
858 int stat_F = gsl_sf_hyperg_2F1_e(a-c+1, b-c+1, -c+2, x, &F);
859 double ln_pre_val = g1.val + g2.val - g3.val - g4.val - g5.val;
860 double ln_pre_err = g1.err + g2.err + g3.err + g4.err + g5.err;
861 double sg = s1 * s2 * s3 * s4 * s5;
862 int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
863 sg * F.val, F.err,
864 result);
865 return GSL_ERROR_SELECT_2(stat_e, stat_F);
866 }
867 }
868 }
869 else {
870 /* generic c */
871 gsl_sf_result F;
872 gsl_sf_result lng;
873 double sgn;
874 int stat_g = gsl_sf_lngamma_sgn_e(c, &lng, &sgn);
875 int stat_F = gsl_sf_hyperg_2F1_e(a, b, c, x, &F);
876 int stat_e = gsl_sf_exp_mult_err_e(-lng.val, lng.err,
877 sgn*F.val, F.err,
878 result);
879 return GSL_ERROR_SELECT_3(stat_e, stat_F, stat_g);
880 }
881 }
882
883
884 int
gsl_sf_hyperg_2F1_conj_renorm_e(const double aR,const double aI,const double c,const double x,gsl_sf_result * result)885 gsl_sf_hyperg_2F1_conj_renorm_e(const double aR, const double aI, const double c,
886 const double x,
887 gsl_sf_result * result
888 )
889 {
890 const double rintc = floor(c + 0.5);
891 const double rinta = floor(aR + 0.5);
892 const int a_neg_integer = ( aR < 0.0 && fabs(aR-rinta) < locEPS && aI == 0.0);
893 const int c_neg_integer = ( c < 0.0 && fabs(c - rintc) < locEPS );
894
895 if(c_neg_integer) {
896 if(a_neg_integer && aR > c+0.1) {
897 /* 2F1 terminates early */
898 result->val = 0.0;
899 result->err = 0.0;
900 return GSL_SUCCESS;
901 }
902 else {
903 /* 2F1 does not terminate early enough, so something survives */
904 /* [Abramowitz+Stegun, 15.1.2] */
905 gsl_sf_result g1, g2;
906 gsl_sf_result g3;
907 gsl_sf_result a1, a2;
908 int stat = 0;
909 stat += gsl_sf_lngamma_complex_e(aR-c+1, aI, &g1, &a1);
910 stat += gsl_sf_lngamma_complex_e(aR, aI, &g2, &a2);
911 stat += gsl_sf_lngamma_e(-c+2.0, &g3);
912 if(stat != 0) {
913 DOMAIN_ERROR(result);
914 }
915 else {
916 gsl_sf_result F;
917 int stat_F = gsl_sf_hyperg_2F1_conj_e(aR-c+1, aI, -c+2, x, &F);
918 double ln_pre_val = 2.0*(g1.val - g2.val) - g3.val;
919 double ln_pre_err = 2.0 * (g1.err + g2.err) + g3.err;
920 int stat_e = gsl_sf_exp_mult_err_e(ln_pre_val, ln_pre_err,
921 F.val, F.err,
922 result);
923 return GSL_ERROR_SELECT_2(stat_e, stat_F);
924 }
925 }
926 }
927 else {
928 /* generic c */
929 gsl_sf_result F;
930 gsl_sf_result lng;
931 double sgn;
932 int stat_g = gsl_sf_lngamma_sgn_e(c, &lng, &sgn);
933 int stat_F = gsl_sf_hyperg_2F1_conj_e(aR, aI, c, x, &F);
934 int stat_e = gsl_sf_exp_mult_err_e(-lng.val, lng.err,
935 sgn*F.val, F.err,
936 result);
937 return GSL_ERROR_SELECT_3(stat_e, stat_F, stat_g);
938 }
939 }
940
941
942 /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
943
944 #include "eval.h"
945
gsl_sf_hyperg_2F1(double a,double b,double c,double x)946 double gsl_sf_hyperg_2F1(double a, double b, double c, double x)
947 {
948 EVAL_RESULT(gsl_sf_hyperg_2F1_e(a, b, c, x, &result));
949 }
950
gsl_sf_hyperg_2F1_conj(double aR,double aI,double c,double x)951 double gsl_sf_hyperg_2F1_conj(double aR, double aI, double c, double x)
952 {
953 EVAL_RESULT(gsl_sf_hyperg_2F1_conj_e(aR, aI, c, x, &result));
954 }
955
gsl_sf_hyperg_2F1_renorm(double a,double b,double c,double x)956 double gsl_sf_hyperg_2F1_renorm(double a, double b, double c, double x)
957 {
958 EVAL_RESULT(gsl_sf_hyperg_2F1_renorm_e(a, b, c, x, &result));
959 }
960
gsl_sf_hyperg_2F1_conj_renorm(double aR,double aI,double c,double x)961 double gsl_sf_hyperg_2F1_conj_renorm(double aR, double aI, double c, double x)
962 {
963 EVAL_RESULT(gsl_sf_hyperg_2F1_conj_renorm_e(aR, aI, c, x, &result));
964 }
965