Re^2: Hypergeometric 2F1() specfunc: cases x=1

Re^2: Hypergeometric 2F1() specfunc: cases x=1

[Richard Mathar]

[-- Attachment #1: Type: text/plain, Size: 613 bytes --]

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

[-- Attachment #2: hyperg_2F1.c.diff --]
[-- Type: text/plain, Size: 1433 bytes --]

--- 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);
   }

[-- Attachment #3: test_hyperg.c.diff --]
[-- Type: text/plain, Size: 2854 bytes --]

--- 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 */

Re: Re^2: Hypergeometric 2F1() specfunc: cases x=1

[Brian Gough]

At Wed, 9 Jan 2008 14:37:51 +0100,
Richard Mathar wrote:
> 
> 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.
> 

Thanks, I've applied that patch.

-- 
Brian Gough