From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 2413 invoked by alias); 9 Jan 2008 13:50:56 -0000 Received: (qmail 27392 invoked by uid 22791); 9 Jan 2008 13:38:41 -0000 X-Spam-Check-By: sourceware.org Date: Wed, 09 Jan 2008 13:50:00 -0000 From: Richard Mathar To: gsl-discuss@sourceware.org Subject: Re^2: Hypergeometric 2F1() specfunc: cases x=1 Message-ID: <20080109133751.GA22080@strw.leidenuniv.nl> MIME-Version: 1.0 Content-Type: multipart/mixed; boundary="X1bOJ3K7DJ5YkBrT" Content-Disposition: inline User-Agent: Mutt/1.5.17 (2007-11-01) X-LeidenObservatory-MailScanner-Information: Processed through MailScanner at Leiden Observatory X-LeidenObservatory-MailScanner: Virusscan: Found to be clean X-LeidenObservatory-MailScanner-SpamCheck: not spam, SpamAssassin (not cached, score=-2.599, required 4.5, autolearn=not spam, BAYES_00 -2.60) X-LeidenObservatory-MailScanner-From: mathar@mail.strw.leidenuniv.nl Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2008-q1/txt/msg00010.txt.bz2 --X1bOJ3K7DJ5YkBrT Content-Type: text/plain; charset=us-ascii Content-Disposition: inline Content-length: 613 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 --X1bOJ3K7DJ5YkBrT Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="hyperg_2F1.c.diff" Content-length: 1433 --- 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); } --X1bOJ3K7DJ5YkBrT Content-Type: text/plain; charset=us-ascii Content-Disposition: attachment; filename="test_hyperg.c.diff" Content-length: 2854 --- 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 */ --X1bOJ3K7DJ5YkBrT--