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

[Ondřej Čertík]

On Mon, Sep 24, 2012 at 8:06 AM, Jonathan Helmus <jjhelmus at gmail.com> wrote:
>
>
>
> 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()).

I think the Fortran implementation is not optimal. One probably has to
use log_gamma()
and logs when doing intermediate calculations to avoid this issue and
only exponentiate the final answer, since we know
it is 0.019816636488030055..., which is a nice number, not too small,
not too big.

Ondrej