Continuing the discussion of
http://sourceware.org/ml/gsl-discuss/2007-q4/msg00038.html :
Two files with the "diff -u" format are attached
which implement Gaussian Hypergeometric Functions 2F1() for argument x=1.
The implementation in hyperg_2F1.c has been changed to handle cases
with negative values of the intermediate Gamma-functions correctly.
All 4 tests added to test_hyperg.c are passed on my gcc 4.1.2 (i386 RedHat).
On an unrelated issue, the testing for some 2F0() has been improved by
insertion of more accurate results in test_hyperg.c.
--
Richard J Mathar
http://www.strw.leidenuniv.nl/~mathar
--- gsl-1.10/specfunc/hyperg_2F1.c.org 2008-01-09 12:52:30.000000000 +0100
+++ gsl-1.10/specfunc/hyperg_2F1.c 2008-01-09 14:00:53.000000000 +0100
@@ -636,6 +636,39 @@
result->val = 0.0;
result->err = 0.0;
+ /* start proposed handling of x=1.0 RJM */
+ if ( fabs(x - 1.0) < locEPS && c-a-b > 0. && c != 0. && ! c_neg_integer )
+ {
+ gsl_sf_result lnGamc, lnGamcab, lnGamca, lnGamcb ;
+ double lnGamcSi, lnGamcaSi, lnGamcbSi ;
+ int stat = gsl_sf_lngamma_sgn_e(c, & lnGamc, & lnGamcSi);
+ if ( stat != GSL_SUCCESS )
+ {
+ DOMAIN_ERROR(result);
+ }
+ stat = gsl_sf_lngamma_e(c-a-b, & lnGamcab);
+ if ( stat != GSL_SUCCESS )
+ {
+ DOMAIN_ERROR(result);
+ }
+ stat = gsl_sf_lngamma_sgn_e(c-a, & lnGamca, & lnGamcaSi);
+ if ( stat != GSL_SUCCESS )
+ {
+ DOMAIN_ERROR(result);
+ }
+ stat = gsl_sf_lngamma_sgn_e(c-b, & lnGamcb, & lnGamcbSi);
+ if ( stat != GSL_SUCCESS )
+ {
+ DOMAIN_ERROR(result);
+ }
+ stat = gsl_sf_exp_err_e(lnGamc.val+lnGamcab.val-lnGamca.val
+
-lnGamcb.val,lnGamc.err+lnGamcab.err+lnGamca.err+lnGamcb.err,result) ;
+
+ result->val *= lnGamcSi/(lnGamcaSi*lnGamcbSi) ;
+ return stat ;
+ }
+/*end proposed handling of x=1.0 RJM */
+
if(x < -1.0 || 1.0 <= x) {
DOMAIN_ERROR(result);
}
--- gsl-1.10/specfunc/test_hyperg.c.org 2008-01-09 13:03:18.000000000 +0100
+++ gsl-1.10/specfunc/test_hyperg.c 2008-01-09 14:22:13.000000000 +0100
@@ -459,6 +459,23 @@
TEST_SF(s, gsl_sf_hyperg_2F1_e, (1.5, 0.5, 2.0, -1.0+1.0/65536.0, &r),
0.762762124908845424449, TEST_TOL0, GSL_SUCCESS);
TEST_SF(s, gsl_sf_hyperg_2F1_e, (1.5, 0.5, 2.0, -1.0+1.0/1048576.0, &r),
0.762759911089064738044, TEST_TOL0, GSL_SUCCESS);
+ /* added special handling with argument = 1.0 , Richard J. Mathar, 2008-01-09
+ * Digits := 30 :
+ * interface(prettyprint=0) ;
+ * a := 1.5; b := 0.5 ; c := 3.0 ;
+ * hypergeom([a,b],[c],1.0) ;
+ * a := 1.5; b := -4.2 ; c := 3.0 ;
+ * hypergeom([a,b],[c],1.0) ;
+ * a := -7.4; b := 0.7 ; c := -1.5 ;
+ * hypergeom([a,b],[c],1.0) ;
+ * a := 0.1; b := -2.7 ; c := -1.5 ;
+ * hypergeom([a,b],[c],1.0) ;
+ */
+ TEST_SF(s, gsl_sf_hyperg_2F1_e, (1.5, 0.5, 3.0, 1.0, &r),
1.6976527263135502482014268 , TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_hyperg_2F1_e, (1.5, -4.2, 3.0, 1.0, &r),
.15583601560025710649555254 , TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_hyperg_2F1_e, (-7.4, 0.7, -1.5, 1.0, &r),
-.34478866959246584996859 , TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_hyperg_2F1_e, (0.1, -2.7, -1.5, 1.0, &r),
1.059766766063610122925 , TEST_TOL0, GSL_SUCCESS);
+
/* 2F1 conj */
@@ -469,14 +486,14 @@
TEST_SF(s, gsl_sf_hyperg_2F1_conj_e, (8, 8, 5, -0.5, &r),
0.00023301916814505902809, TEST_TOL3, GSL_SUCCESS);
TEST_SF(s, gsl_sf_hyperg_2F1_conj_e, (25, 25, 1, -0.5, &r),
5.1696944096e-06, TEST_SQRT_TOL0, GSL_SUCCESS);
- /* FIXME: the "true" values here may not be so good */
- /*
- TEST_SF(s, gsl_sf_hyperg_2F0_e, (0.01, 1.0, -0.02, &r), 0.999803886708565
, TEST_TOL0, GSL_SUCCESS);
- TEST_SF(s, gsl_sf_hyperg_2F0_e, (0.1, 0.5, -0.02, &r), 0.999015947934831
, TEST_TOL0, GSL_SUCCESS);
- TEST_SF(s, gsl_sf_hyperg_2F0_e, (1, 1, -0.02, &r), 0.980755496569062 ,
TEST_TOL0, GSL_SUCCESS);
- TEST_SF(s, gsl_sf_hyperg_2F0_e, (8, 8, -0.02, &r), 0.3299059284994299 ,
TEST_TOL0, GSL_SUCCESS);
- TEST_SF(s, gsl_sf_hyperg_2F0_e, (50, 50, -0.02, &r), 2.688995263773233e-13,
TEST_TOL0, GSL_SUCCESS);
- */
+
+ /* updated correct values, testing enabled, Richard J. Mathar, 2008-01-09 */
+
+ TEST_SF(s, gsl_sf_hyperg_2F0_e, (0.01, 1.0, -0.02, &r),
.99980388665511730901180717 , TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_hyperg_2F0_e, (0.1, 0.5, -0.02, &r),
.99901595171179281891589794 , TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_hyperg_2F0_e, (1, 1, -0.02, &r),
.98075549650574351826538049000 , TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_hyperg_2F0_e, (8, 8, -0.02, &r),
.32990592849626965538692141 , TEST_TOL0, GSL_SUCCESS);
+ TEST_SF(s, gsl_sf_hyperg_2F0_e, (50, 50, -0.02, &r),
.2688995263772964415245902e-12 , TEST_TOL0, GSL_SUCCESS);
/* 2F1 renorm */
