Hello,
I use the function 'hypgeo(a,b,c,z)' given in section 6.12 to evaluate the hypergeometric function with complex arguments.
It works pretty well for a large range of values of the arguments, but fails for others.

For instance, with the values
a = cmplx(3.02619478964082,0.366612579197369)
b = cmplx(1.97380521035918,0.366612579197369)
c = cmplx(3.00000000000000,0.0000000000E+000)
z = cmplx(1.21104634938969,0.0000000000E+000)
I get the following answer :
Fortran : (1.1694339666388947,1.3493422205836327)
Python : (108.68967630069800,155.17225602885699)
with both f77 and f90 versions of the Numerical Recipes code.

Mathematica gives the same answer as Python, the error therefore comes from the Fortran.

I did not completely figure out what is made in NR to evaluate the function, but looking at the code I guess the error comes from the insufficient precision used for the integration performed. I tried to increase the precision with the following steps:

changing EPS=1.e-15 and hmin=1.e-15 in subroutine 'hypgeo'
--> no change
EPS=1.e-20.
--> evaluation fails : STOP stepsize smaller than minimum in odeint
hmin=1.e-20 when calling odeint
-->STOP stepsize underflow in bsstep
compilation with real*8 option
--> STOP too many steps in odeint
MAXSTP=1e7 in odeint
--> evaluation takes a very long time, but no change at the end


What else can I do to get the correct evaluation of the hypergeometric function ?
Thanks very much for your help.

Hi Mazkime,

The answer to this one is a little tricky. The hypergeometric function F has a branch cut in the complex plane that starts at z=1 and is usually taken to run along the real axis to plus infinity. This means that the value of F for z real and bigger than 1 has to be defined by convention: either as the limit approaching the real axis from the upper half plane or from below. The standard is to work with F for which arg(-z) < pi, which you can check means approaching the real axis from below for z > 1.

Unfortunately in NR we implemented this incorrectly: the values are being returned for approaching the real axis from above. You can verify this by adding a tiny positive imaginary piece to z and confirming that Mathematica and NR now agree.

The fix is to change line 13 of the function hypgeo from
else z0=Complex(0.0,imag(z) >= 0.0 ? 0.5 : -0.5);
to
else z0=Complex(0.0,imag(z) > 0.0 ? 0.5 : -0.5);

Thanks for letting us know about this.

Saul Teukolsky

Hi Saul,
thanks very much for the answer and the very clear explanations.
After correction of the routine, everything works fine.

Cheers,
Mazkime

Hello again,
I encountered a new problem while using the hypergeometric function. I am using the f77 version of NR.
An example is given with the following values :

a = cmplx( 22.3604742608975 ,-46.6320376130814)
b = cmplx(-19.3604742608975 ,-46.6320376130814)
c = cmplx( 2.00000000000000 , 0.000000000000000E+000)
z = cmplx( 0.680320117079294, 0.000000000000000E+000)


I get the following answers :

fortran : (-7.0194405E+22,-3.0385708E+22)
python : (-2.9433980E-03,4.2219213E-03)

Mathematica gives the same result as python, and fortran clearly seems to diverge.
The previous fix given in post #2 does not change the fortran result, as real(z) < 1.

How can I fix this issue ?
Many thanks for your help.
Mazkime

Dear Mazkime,

As the text in the book explains, the routine hypgeo is meant to illustrate a particular technique and is not a universal routine for computing the hypergeometric function.

The problem you are finding is because for real values of z less than 1, the routine initializes the ode integrator by computing the hypergeometric series for z=0.5. But for "large" values of a and b (eg of order 20 as in your example), the terms in the series get very large and cancel before converging to a small answer. The result is that with only 16 significant digits, no accuracy is left in the final answer. (This is the same behavior you get when evaluating sin x by summing the series when x is too big.)

The fix is to set the starting value z0 to a smaller value, say 0.01 (or even 0). Your example will then work. But coming up with good values of z0 for all values of a,b,c is very complicated.

Thanks for writing.
Saul Teukolsky

Saul,
thank you for your answer. The problem is now very clear to me.
Indeed, selecting a different value of z0 fixes the issue for this specific case, but choosing a good value of z0 in the general case is very difficult.
Another solution is to compile the code with the -r16 option (real*16 for increased precision), which drastically improves the evaluation of the hypergeometric function with the initial value of z0 (the error is only 4% in that case).

Regards,
Mazkime