[SciPy-User] hyp1f1(0.5, 1.5, -2000) fails

[Jonathan Helmus]

On 09/24/2012 08:15 AM, Moore, Eric (NIH/NIDDK) [F] wrote:
>> -----Original Message-----
>> From: Ondrej Certik [mailto:ondrej at certik.cz]
>> Sent: Sunday, September 23, 2012 7:32 PM
>> To: SciPy Users List
>> Subject: [SciPy-User] hyp1f1(0.5, 1.5, -2000) fails
>>
>> Hi,
>>
>> I noticed, that hyp1f1(0.5, 1.5, -2000) fails:
>>
>> In [5]: scipy.special.hyp1f1(0.5, 1.5, 0)
>> Out[5]: 1.0
>>
>> In [6]: scipy.special.hyp1f1(0.5, 1.5, -20)
>> Out[6]: 0.19816636482997368
>>
>> In [7]: scipy.special.hyp1f1(0.5, 1.5, -200)
>> Out[7]: 0.062665706865775023
>>
>> In [8]: scipy.special.hyp1f1(0.5, 1.5, -2000)
>> Warning: invalid value encountered in hyp1f1
>> Out[8]: nan
>>
>>
>> The values [5], [6] and [7] are correct. The value [8] should be:
>>
>> 0.019816636488030055...
>>
>> Ondrej
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
> Its blowing up around -709:
>
> In [60]: s.hyp1f1(0.5, 1.5, -709.7827128933)
> Out[60]: 0.03326459435722777
>
> In [61]: s.hyp1f1(0.5, 1.5, -709.7827128934)
> Out[61]: inf
>
> Eric
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

exp(709.78) is right at the edge of causing a overflow in double 
precision, 709.79 causes an overflow.


In [2]: np.exp(np.array([709.78], dtype='float64'))
Out[2]: array([  1.79282279e+308])

In [3]: np.exp(np.array([709.79], dtype='float64'))
/home/jhelmus/bin/ipython:1: RuntimeWarning: overflow encountered in exp
   #!/home/jhelmus/bin/epd-7.3-1-rh5-x86_64/bin/python
Out[3]: array([ inf])

Line 5695 in specfun.f is calculating an exponential (and I believe is 
causing this overflow).  The complex version of this subroutine (CCHG) 
also suffers from this problem.  I do not know enough about this type of 
function to suggest a meaningful fix other than suggesting a check that 
abs(x3) <= np.log(np.finfo('d').max()).

- Jonathan Helmus