Re: [Bug-gsl] Re: [Help-gsl] Hypergeometric function

[Imad-Eddine Srairi]

Hi,

I looked again at the code of function "hyperg_2F1_reflect" found in 
file hyperg_2F1.c .

As I had noticed and reported in my first e-mail regarding this topic, 
this function handles separately the case where c - a - b = m is an 
integer (which is the case for the bug we are interested in: a = -0.5 ; 
b = 1.5 ; c = 1.0) and computes:

F1 + sgn * F2

where sgn depends on parity of the integer m (sgn = -1 if m is odd, +1 
otherwise).

In the test case we are looking into, F1 = 0.
I am now also pretty sure that F2 is the second part (i.e. the infinite 
sum) of 15.3.11 in Abramowitz&Stegun [1] and of 15.8.10 in the DLMF [2]

There are only two differences. The author of the code has chosen to 
distribute the minus sign: instead of

ln(1-x)-psi(k+1)-psi(k+m+1)+psi(a+k+m)+psi(b+k+m)

(where psi is the digamma function [3]) s/he computes

psi(k+1)+psi(k+m+1)-psi(a+k+m)-psi(b+k+m)-ln(1-x)

I can see nothing wrong with this.

The other difference lies within the factor before the infinite sum.
Both A&S and the DLMF write (x - 1)^m whereas the corresponding code is

double ln_pre2_val = lng_c.val + d1*ln_omx - lng_ad2.val - lng_bd2.val;

with

const double ln_omx = log(1.0 - x);
(omx probably stands for "one minus x")

This is the origin of "sgn" in "F1 + sgn * F2": we need to multiply by 
(-1)^m since (x - 1)^m = (-1)^m * (1 - x)^m
But this time again, I can see nothing wrong with this.

So as far as I can see, the code is correct or at least corresponds to 
the formulas found in A & S as well as in the DLMF.
But yet the result is wrong...
I am quite puzzled...

imad

[1] http://people.math.sfu.ca/~cbm/aands/page_559.htm
[2] http://dlmf.nist.gov/15.8.E10
[3] http://people.math.sfu.ca/~cbm/aands/page_258.htm

On 02/11/2011 06:00 PM, Brian Gough wrote:
> At Fri, 04 Feb 2011 17:19:32 +0100,
> Imad-Eddine Srairi wrote:
> > /* Do the reflection described in [Moshier, p. 334].
> >    * Assumes a,b,c != neg integer.
> >    */
> >
> > I have no access to this book, but I am currently trying to match the
> > code against the special cases described in [2].
>
> Thanks for tracking that down.  I have the book so I will look at it
> and see if there is any more information but Moshier doesn't usually
> have much explanation of the algorithms.