1 /* hyp2f1.c
2 *
3 * Gauss hypergeometric function F
4 * 2 1
5 *
6 *
7 * SYNOPSIS:
8 *
9 * double a, b, c, x, y, hyp2f1();
10 *
11 * y = hyp2f1( a, b, c, x );
12 *
13 *
14 * DESCRIPTION:
15 *
16 *
17 * hyp2f1( a, b, c, x ) = F ( a, b; c; x )
18 * 2 1
19 *
20 * inf.
21 * - a(a+1)...(a+k) b(b+1)...(b+k) k+1
22 * = 1 + > ----------------------------- x .
23 * - c(c+1)...(c+k) (k+1)!
24 * k = 0
25 *
26 * Cases addressed are
27 * Tests and escapes for negative integer a, b, or c
28 * Linear transformation if c - a or c - b negative integer
29 * Special case c = a or c = b
30 * Linear transformation for x near +1
31 * Transformation for x < -0.5
32 * Psi function expansion if x > 0.5 and c - a - b integer
33 * Conditionally, a recurrence on c to make c-a-b > 0
34 *
35 * |x| > 1 is rejected.
36 *
37 * The parameters a, b, c are considered to be integer
38 * valued if they are within 1.0e-14 of the nearest integer
39 * (1.0e-13 for IEEE arithmetic).
40 *
41 * ACCURACY:
42 *
43 *
44 * Relative error (-1 < x < 1):
45 * arithmetic domain # trials peak rms
46 * IEEE -1,7 230000 1.2e-11 5.2e-14
47 *
48 * Several special cases also tested with a, b, c in
49 * the range -7 to 7.
50 *
51 * ERROR MESSAGES:
52 *
53 * A "partial loss of precision" message is printed if
54 * the internally estimated relative error exceeds 1^-12.
55 * A "singularity" message is printed on overflow or
56 * in cases not addressed (such as x < -1).
57 */
58
59 /* hyp2f1 */
60
61
62 /*
63 Cephes Math Library Release 2.2: November, 1992
64 Copyright 1984, 1987, 1992 by Stephen L. Moshier
65 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
66 */
67
68
69 #include "mconf.h"
70 #include "cephes.h"
71
72 #ifdef DEC
73 #define EPS 1.0e-14
74 #define EPS2 1.0e-11
75 #endif
76
77 #ifdef IBMPC
78 #define EPS 1.0e-13
79 #define EPS2 1.0e-10
80 #endif
81
82 #ifdef MIEEE
83 #define EPS 1.0e-13
84 #define EPS2 1.0e-10
85 #endif
86
87 #ifdef UNK
88 #define EPS 1.0e-13
89 #define EPS2 1.0e-10
90 #endif
91
92 #define ETHRESH 1.0e-12
93
94 static double hyt2f1(double, double, double, double, double *);
95 static double hys2f1(double, double, double, double, double *);
96
97 extern double MAXNUM, MACHEP;
98
hyp2f1(a,b,c,x)99 double hyp2f1( a, b, c, x )
100 double a, b, c, x;
101 {
102 double d, d1, d2, e;
103 double p, q, r, s, y, ax;
104 double ia, ib, ic, id, err;
105 int flag, i, aid;
106
107 err = 0.0;
108 ax = fabs(x);
109 s = 1.0 - x;
110 flag = 0;
111 ia = round(a); /* nearest integer to a */
112 ib = round(b);
113
114 if( a <= 0 )
115 {
116 if( fabs(a-ia) < EPS ) /* a is a negative integer */
117 flag |= 1;
118 }
119
120 if( b <= 0 )
121 {
122 if( fabs(b-ib) < EPS ) /* b is a negative integer */
123 flag |= 2;
124 }
125
126 if( ax < 1.0 )
127 {
128 if( fabs(b-c) < EPS ) /* b = c */
129 {
130 y = pow( s, -a ); /* s to the -a power */
131 goto hypdon;
132 }
133 if( fabs(a-c) < EPS ) /* a = c */
134 {
135 y = pow( s, -b ); /* s to the -b power */
136 goto hypdon;
137 }
138 }
139
140
141
142 if( c <= 0.0 )
143 {
144 ic = round(c); /* nearest integer to c */
145 if( fabs(c-ic) < EPS ) /* c is a negative integer */
146 {
147 /* check if termination before explosion */
148 if( (flag & 1) && (ia > ic) )
149 goto hypok;
150 if( (flag & 2) && (ib > ic) )
151 goto hypok;
152 goto hypdiv;
153 }
154 }
155
156 if( flag ) /* function is a polynomial */
157 goto hypok;
158
159 if( ax > 1.0 ) /* series diverges */
160 goto hypdiv;
161
162 p = c - a;
163 ia = round(p); /* nearest integer to c-a */
164 if( (ia <= 0.0) && (fabs(p-ia) < EPS) ) /* negative int c - a */
165 flag |= 4;
166
167 r = c - b;
168 ib = round(r); /* nearest integer to c-b */
169 if( (ib <= 0.0) && (fabs(r-ib) < EPS) ) /* negative int c - b */
170 flag |= 8;
171
172 d = c - a - b;
173 id = round(d); /* nearest integer to d */
174 q = fabs(d-id);
175
176 /* Thanks to Christian Burger <BURGER@DMRHRZ11.HRZ.Uni-Marburg.DE>
177 * for reporting a bug here. */
178 if( fabs(ax-1.0) < EPS ) /* |x| == 1.0 */
179 {
180 if( x > 0.0 )
181 {
182 if( flag & 12 ) /* negative int c-a or c-b */
183 {
184 if( d >= 0.0 )
185 goto hypf;
186 else
187 goto hypdiv;
188 }
189 if( d <= 0.0 )
190 goto hypdiv;
191 y = true_gamma(c)*true_gamma(d)/(true_gamma(p)*true_gamma(r));
192 goto hypdon;
193 }
194
195 if( d <= -1.0 )
196 goto hypdiv;
197
198 }
199
200 /* Conditionally make d > 0 by recurrence on c
201 * AMS55 #15.2.27
202 */
203 if( d < 0.0 )
204 {
205 /* Try the power series first */
206 y = hyt2f1( a, b, c, x, &err );
207 if( err < ETHRESH )
208 goto hypdon;
209 /* Apply the recurrence if power series fails */
210 err = 0.0;
211 aid = 2 - id;
212 e = c + aid;
213 d2 = hyp2f1(a,b,e,x);
214 d1 = hyp2f1(a,b,e+1.0,x);
215 q = a + b + 1.0;
216 for( i=0; i<aid; i++ )
217 {
218 r = e - 1.0;
219 y = (e*(r-(2.0*e-q)*x)*d2 + (e-a)*(e-b)*x*d1)/(e*r*s);
220 e = r;
221 d1 = d2;
222 d2 = y;
223 }
224 goto hypdon;
225 }
226
227
228 if( flag & 12 )
229 goto hypf; /* negative integer c-a or c-b */
230
231 hypok:
232 y = hyt2f1( a, b, c, x, &err );
233
234
235 hypdon:
236 if( err > ETHRESH )
237 {
238 mtherr( "hyp2f1", PLOSS );
239 /* printf( "Estimated err = %.2e\n", err ); */
240 }
241 return(y);
242
243 /* The transformation for c-a or c-b negative integer
244 * AMS55 #15.3.3
245 */
246 hypf:
247 y = pow( s, d ) * hys2f1( c-a, c-b, c, x, &err );
248 goto hypdon;
249
250 /* The alarm exit */
251 hypdiv:
252 mtherr( "hyp2f1", OVERFLOW );
253 return( MAXNUM );
254 }
255
256
257
258
259
260
261 /* Apply transformations for |x| near 1
262 * then call the power series
263 */
hyt2f1(a,b,c,x,loss)264 static double hyt2f1( a, b, c, x, loss )
265 double a, b, c, x;
266 double *loss;
267 {
268 double p, q, r, s, t, y, d, err, err1;
269 double ax, id, d1, d2, e, y1;
270 int i, aid;
271
272 err = 0.0;
273 s = 1.0 - x;
274 if( x < -0.5 )
275 {
276 if( b > a )
277 y = pow( s, -a ) * hys2f1( a, c-b, c, -x/s, &err );
278
279 else
280 y = pow( s, -b ) * hys2f1( c-a, b, c, -x/s, &err );
281
282 goto done;
283 }
284
285 d = c - a - b;
286 id = round(d); /* nearest integer to d */
287
288 if( x > 0.9 )
289 {
290 if( fabs(d-id) > EPS ) /* test for integer c-a-b */
291 {
292 /* Try the power series first */
293 y = hys2f1( a, b, c, x, &err );
294 if( err < ETHRESH )
295 goto done;
296 /* If power series fails, then apply AMS55 #15.3.6 */
297 q = hys2f1( a, b, 1.0-d, s, &err );
298 q *= true_gamma(d) /(true_gamma(c-a) * true_gamma(c-b));
299 r = pow(s,d) * hys2f1( c-a, c-b, d+1.0, s, &err1 );
300 r *= true_gamma(-d)/(true_gamma(a) * true_gamma(b));
301 y = q + r;
302
303 q = fabs(q); /* estimate cancellation error */
304 r = fabs(r);
305 if( q > r )
306 r = q;
307 err += err1 + (MACHEP*r)/y;
308
309 y *= true_gamma(c);
310 goto done;
311 }
312 else
313 {
314 /* Psi function expansion, AMS55 #15.3.10, #15.3.11, #15.3.12 */
315 if( id >= 0.0 )
316 {
317 e = d;
318 d1 = d;
319 d2 = 0.0;
320 aid = id;
321 }
322 else
323 {
324 e = -d;
325 d1 = 0.0;
326 d2 = d;
327 aid = -id;
328 }
329
330 ax = log(s);
331
332 /* sum for t = 0 */
333 y = psi(1.0) + psi(1.0+e) - psi(a+d1) - psi(b+d1) - ax;
334 y /= true_gamma(e+1.0);
335
336 p = (a+d1) * (b+d1) * s / true_gamma(e+2.0); /* Poch for t=1 */
337 t = 1.0;
338 do
339 {
340 r = psi(1.0+t) + psi(1.0+t+e) - psi(a+t+d1)
341 - psi(b+t+d1) - ax;
342 q = p * r;
343 y += q;
344 p *= s * (a+t+d1) / (t+1.0);
345 p *= (b+t+d1) / (t+1.0+e);
346 t += 1.0;
347 }
348 while( fabs(q/y) > EPS );
349
350
351 if( id == 0.0 )
352 {
353 y *= true_gamma(c)/(true_gamma(a)*true_gamma(b));
354 goto psidon;
355 }
356
357 y1 = 1.0;
358
359 if( aid == 1 )
360 goto nosum;
361
362 t = 0.0;
363 p = 1.0;
364 for( i=1; i<aid; i++ )
365 {
366 r = 1.0-e+t;
367 p *= s * (a+t+d2) * (b+t+d2) / r;
368 t += 1.0;
369 p /= t;
370 y1 += p;
371 }
372 nosum:
373 p = true_gamma(c);
374 y1 *= true_gamma(e) * p / (true_gamma(a+d1) * true_gamma(b+d1));
375
376 y *= p / (true_gamma(a+d2) * true_gamma(b+d2));
377 if( (aid & 1) != 0 )
378 y = -y;
379
380 q = pow( s, id ); /* s to the id power */
381 if( id > 0.0 )
382 y *= q;
383 else
384 y1 *= q;
385
386 y += y1;
387 psidon:
388 goto done;
389 }
390
391 }
392
393 /* Use defining power series if no special cases */
394 y = hys2f1( a, b, c, x, &err );
395
396 done:
397 *loss = err;
398 return(y);
399 }
400
401
402
403
404
405 /* Defining power series expansion of Gauss hypergeometric function */
406
hys2f1(a,b,c,x,loss)407 static double hys2f1( a, b, c, x, loss )
408 double a, b, c, x;
409 double *loss; /* estimates loss of significance */
410 {
411 double f, g, h, k, m, s, u, umax;
412 int i;
413
414 i = 0;
415 umax = 0.0;
416 f = a;
417 g = b;
418 h = c;
419 s = 1.0;
420 u = 1.0;
421 k = 0.0;
422 do
423 {
424 if( fabs(h) < EPS )
425 {
426 *loss = 1.0;
427 return( MAXNUM );
428 }
429 m = k + 1.0;
430 u = u * ((f+k) * (g+k) * x / ((h+k) * m));
431 s += u;
432 k = fabs(u); /* remember largest term summed */
433 if( k > umax )
434 umax = k;
435 k = m;
436 if( ++i > 10000 ) /* should never happen */
437 {
438 *loss = 1.0;
439 return(s);
440 }
441 }
442 while( fabs(u/s) > MACHEP );
443
444 /* return estimated relative error */
445 *loss = (MACHEP*umax)/fabs(s) + (MACHEP*i);
446
447 return(s);
448 }
449