hyp2f1f.c 7.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442
  1. /* hyp2f1f.c
  2. *
  3. * Gauss hypergeometric function F
  4. * 2 1
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * float a, b, c, x, y, hyp2f1f();
  10. *
  11. * y = hyp2f1f( 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-6 of the nearest integer.
  39. *
  40. * ACCURACY:
  41. *
  42. * Relative error (-1 < x < 1):
  43. * arithmetic domain # trials peak rms
  44. * IEEE 0,3 30000 5.8e-4 4.3e-6
  45. */
  46. /* hyp2f1 */
  47. /*
  48. Cephes Math Library Release 2.2: November, 1992
  49. Copyright 1984, 1987, 1992 by Stephen L. Moshier
  50. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  51. */
  52. #include <math.h>
  53. #define EPS 1.0e-5
  54. #define EPS2 1.0e-5
  55. #define ETHRESH 1.0e-5
  56. extern float MAXNUMF, MACHEPF;
  57. #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
  58. #ifdef ANSIC
  59. float powf(float, float);
  60. static float hys2f1f(float, float, float, float, float *);
  61. static float hyt2f1f(float, float, float, float, float *);
  62. float gammaf(float), logf(float), expf(float), psif(float);
  63. float floorf(float);
  64. #else
  65. float powf(), gammaf(), logf(), expf(), psif();
  66. float floorf();
  67. static float hyt2f1f(), hys2f1f();
  68. #endif
  69. #define roundf(x) (floorf((x)+(float )0.5))
  70. float hyp2f1f( float aa, float bb, float cc, float xx )
  71. {
  72. float a, b, c, x;
  73. float d, d1, d2, e;
  74. float p, q, r, s, y, ax;
  75. float ia, ib, ic, id, err;
  76. int flag, i, aid;
  77. a = aa;
  78. b = bb;
  79. c = cc;
  80. x = xx;
  81. err = 0.0;
  82. ax = fabsf(x);
  83. s = 1.0 - x;
  84. flag = 0;
  85. ia = roundf(a); /* nearest integer to a */
  86. ib = roundf(b);
  87. if( a <= 0 )
  88. {
  89. if( fabsf(a-ia) < EPS ) /* a is a negative integer */
  90. flag |= 1;
  91. }
  92. if( b <= 0 )
  93. {
  94. if( fabsf(b-ib) < EPS ) /* b is a negative integer */
  95. flag |= 2;
  96. }
  97. if( ax < 1.0 )
  98. {
  99. if( fabsf(b-c) < EPS ) /* b = c */
  100. {
  101. y = powf( s, -a ); /* s to the -a power */
  102. goto hypdon;
  103. }
  104. if( fabsf(a-c) < EPS ) /* a = c */
  105. {
  106. y = powf( s, -b ); /* s to the -b power */
  107. goto hypdon;
  108. }
  109. }
  110. if( c <= 0.0 )
  111. {
  112. ic = roundf(c); /* nearest integer to c */
  113. if( fabsf(c-ic) < EPS ) /* c is a negative integer */
  114. {
  115. /* check if termination before explosion */
  116. if( (flag & 1) && (ia > ic) )
  117. goto hypok;
  118. if( (flag & 2) && (ib > ic) )
  119. goto hypok;
  120. goto hypdiv;
  121. }
  122. }
  123. if( flag ) /* function is a polynomial */
  124. goto hypok;
  125. if( ax > 1.0 ) /* series diverges */
  126. goto hypdiv;
  127. p = c - a;
  128. ia = roundf(p);
  129. if( (ia <= 0.0) && (fabsf(p-ia) < EPS) ) /* negative int c - a */
  130. flag |= 4;
  131. r = c - b;
  132. ib = roundf(r); /* nearest integer to r */
  133. if( (ib <= 0.0) && (fabsf(r-ib) < EPS) ) /* negative int c - b */
  134. flag |= 8;
  135. d = c - a - b;
  136. id = roundf(d); /* nearest integer to d */
  137. q = fabsf(d-id);
  138. if( fabsf(ax-1.0) < EPS ) /* |x| == 1.0 */
  139. {
  140. if( x > 0.0 )
  141. {
  142. if( flag & 12 ) /* negative int c-a or c-b */
  143. {
  144. if( d >= 0.0 )
  145. goto hypf;
  146. else
  147. goto hypdiv;
  148. }
  149. if( d <= 0.0 )
  150. goto hypdiv;
  151. y = gammaf(c)*gammaf(d)/(gammaf(p)*gammaf(r));
  152. goto hypdon;
  153. }
  154. if( d <= -1.0 )
  155. goto hypdiv;
  156. }
  157. /* Conditionally make d > 0 by recurrence on c
  158. * AMS55 #15.2.27
  159. */
  160. if( d < 0.0 )
  161. {
  162. /* Try the power series first */
  163. y = hyt2f1f( a, b, c, x, &err );
  164. if( err < ETHRESH )
  165. goto hypdon;
  166. /* Apply the recurrence if power series fails */
  167. err = 0.0;
  168. aid = 2 - id;
  169. e = c + aid;
  170. d2 = hyp2f1f(a,b,e,x);
  171. d1 = hyp2f1f(a,b,e+1.0,x);
  172. q = a + b + 1.0;
  173. for( i=0; i<aid; i++ )
  174. {
  175. r = e - 1.0;
  176. y = (e*(r-(2.0*e-q)*x)*d2 + (e-a)*(e-b)*x*d1)/(e*r*s);
  177. e = r;
  178. d1 = d2;
  179. d2 = y;
  180. }
  181. goto hypdon;
  182. }
  183. if( flag & 12 )
  184. goto hypf; /* negative integer c-a or c-b */
  185. hypok:
  186. y = hyt2f1f( a, b, c, x, &err );
  187. hypdon:
  188. if( err > ETHRESH )
  189. {
  190. mtherr( "hyp2f1", PLOSS );
  191. /* printf( "Estimated err = %.2e\n", err );*/
  192. }
  193. return(y);
  194. /* The transformation for c-a or c-b negative integer
  195. * AMS55 #15.3.3
  196. */
  197. hypf:
  198. y = powf( s, d ) * hys2f1f( c-a, c-b, c, x, &err );
  199. goto hypdon;
  200. /* The alarm exit */
  201. hypdiv:
  202. mtherr( "hyp2f1f", OVERFLOW );
  203. return( MAXNUMF );
  204. }
  205. /* Apply transformations for |x| near 1
  206. * then call the power series
  207. */
  208. static float hyt2f1f( float aa, float bb, float cc, float xx, float *loss )
  209. {
  210. float a, b, c, x;
  211. float p, q, r, s, t, y, d, err, err1;
  212. float ax, id, d1, d2, e, y1;
  213. int i, aid;
  214. a = aa;
  215. b = bb;
  216. c = cc;
  217. x = xx;
  218. err = 0.0;
  219. s = 1.0 - x;
  220. if( x < -0.5 )
  221. {
  222. if( b > a )
  223. y = powf( s, -a ) * hys2f1f( a, c-b, c, -x/s, &err );
  224. else
  225. y = powf( s, -b ) * hys2f1f( c-a, b, c, -x/s, &err );
  226. goto done;
  227. }
  228. d = c - a - b;
  229. id = roundf(d); /* nearest integer to d */
  230. if( x > 0.8 )
  231. {
  232. if( fabsf(d-id) > EPS2 ) /* test for integer c-a-b */
  233. {
  234. /* Try the power series first */
  235. y = hys2f1f( a, b, c, x, &err );
  236. if( err < ETHRESH )
  237. goto done;
  238. /* If power series fails, then apply AMS55 #15.3.6 */
  239. q = hys2f1f( a, b, 1.0-d, s, &err );
  240. q *= gammaf(d) /(gammaf(c-a) * gammaf(c-b));
  241. r = powf(s,d) * hys2f1f( c-a, c-b, d+1.0, s, &err1 );
  242. r *= gammaf(-d)/(gammaf(a) * gammaf(b));
  243. y = q + r;
  244. q = fabsf(q); /* estimate cancellation error */
  245. r = fabsf(r);
  246. if( q > r )
  247. r = q;
  248. err += err1 + (MACHEPF*r)/y;
  249. y *= gammaf(c);
  250. goto done;
  251. }
  252. else
  253. {
  254. /* Psi function expansion, AMS55 #15.3.10, #15.3.11, #15.3.12 */
  255. if( id >= 0.0 )
  256. {
  257. e = d;
  258. d1 = d;
  259. d2 = 0.0;
  260. aid = id;
  261. }
  262. else
  263. {
  264. e = -d;
  265. d1 = 0.0;
  266. d2 = d;
  267. aid = -id;
  268. }
  269. ax = logf(s);
  270. /* sum for t = 0 */
  271. y = psif(1.0) + psif(1.0+e) - psif(a+d1) - psif(b+d1) - ax;
  272. y /= gammaf(e+1.0);
  273. p = (a+d1) * (b+d1) * s / gammaf(e+2.0); /* Poch for t=1 */
  274. t = 1.0;
  275. do
  276. {
  277. r = psif(1.0+t) + psif(1.0+t+e) - psif(a+t+d1)
  278. - psif(b+t+d1) - ax;
  279. q = p * r;
  280. y += q;
  281. p *= s * (a+t+d1) / (t+1.0);
  282. p *= (b+t+d1) / (t+1.0+e);
  283. t += 1.0;
  284. }
  285. while( fabsf(q/y) > EPS );
  286. if( id == 0.0 )
  287. {
  288. y *= gammaf(c)/(gammaf(a)*gammaf(b));
  289. goto psidon;
  290. }
  291. y1 = 1.0;
  292. if( aid == 1 )
  293. goto nosum;
  294. t = 0.0;
  295. p = 1.0;
  296. for( i=1; i<aid; i++ )
  297. {
  298. r = 1.0-e+t;
  299. p *= s * (a+t+d2) * (b+t+d2) / r;
  300. t += 1.0;
  301. p /= t;
  302. y1 += p;
  303. }
  304. nosum:
  305. p = gammaf(c);
  306. y1 *= gammaf(e) * p / (gammaf(a+d1) * gammaf(b+d1));
  307. y *= p / (gammaf(a+d2) * gammaf(b+d2));
  308. if( (aid & 1) != 0 )
  309. y = -y;
  310. q = powf( s, id ); /* s to the id power */
  311. if( id > 0.0 )
  312. y *= q;
  313. else
  314. y1 *= q;
  315. y += y1;
  316. psidon:
  317. goto done;
  318. }
  319. }
  320. /* Use defining power series if no special cases */
  321. y = hys2f1f( a, b, c, x, &err );
  322. done:
  323. *loss = err;
  324. return(y);
  325. }
  326. /* Defining power series expansion of Gauss hypergeometric function */
  327. static float hys2f1f( float aa, float bb, float cc, float xx, float *loss )
  328. {
  329. int i;
  330. float a, b, c, x;
  331. float f, g, h, k, m, s, u, umax;
  332. a = aa;
  333. b = bb;
  334. c = cc;
  335. x = xx;
  336. i = 0;
  337. umax = 0.0;
  338. f = a;
  339. g = b;
  340. h = c;
  341. k = 0.0;
  342. s = 1.0;
  343. u = 1.0;
  344. do
  345. {
  346. if( fabsf(h) < EPS )
  347. return( MAXNUMF );
  348. m = k + 1.0;
  349. u = u * ((f+k) * (g+k) * x / ((h+k) * m));
  350. s += u;
  351. k = fabsf(u); /* remember largest term summed */
  352. if( k > umax )
  353. umax = k;
  354. k = m;
  355. if( ++i > 10000 ) /* should never happen */
  356. {
  357. *loss = 1.0;
  358. return(s);
  359. }
  360. }
  361. while( fabsf(u/s) > MACHEPF );
  362. /* return estimated relative error */
  363. *loss = (MACHEPF*umax)/fabsf(s) + (MACHEPF*i);
  364. return(s);
  365. }