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
-----Original Message----- From: Ondrej Certik [mailto:ondrej@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@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
On 09/24/2012 08:15 AM, Moore, Eric (NIH/NIDDK) [F] wrote:
-----Original Message----- From: Ondrej Certik [mailto:ondrej@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@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@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
On Mon, Sep 24, 2012 at 8:06 AM, Jonathan Helmus <jjhelmus@gmail.com> wrote:
On 09/24/2012 08:15 AM, Moore, Eric (NIH/NIDDK) [F] wrote:
-----Original Message----- From: Ondrej Certik [mailto:ondrej@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@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@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
participants (4)
-
Jonathan Helmus -
Moore, Eric (NIH/NIDDK) [F] -
Ondrej Certik -
Ondřej Čertík