#18210: numerical bug in incomplete gamma function
-------------------------------------------------+-------------------------
Reporter: VivianePons | Owner:
Type: defect | Status: new
Priority: major | Milestone: sage-6.7
Component: symbolics | Resolution:
Keywords: sd67 | Merged in:
Authors: | Reviewers:
Report Upstream: Reported upstream. Developers | Work issues:
acknowledge bug. | Commit:
Branch: | Stopgaps:
Dependencies: |
-------------------------------------------------+-------------------------
Comment (by buck):
Discusion migrated from #16697
----
by buck:
... this patch doesn't fix the numerical instability in incomplete gamma
as promised elsewhere.
{{{
sage: gamma_inc_lower(25, 14.5).n()
-6.63538851954338e22
sage: gamma_inc_lower(25, 14.5).n(algorithm='mpmath')
-6.63538851954338e22
sage: gamma_inc_lower(25, 14.5).n(algorithm='pari')
-6.63538851954338e22
}}}
by buck:
Of note, gamma (even incomplete) is strictly positive in the reals; it's
an integral of a strictly positive function.
https://en.wikipedia.org/wiki/Incomplete_gamma_function#Definition
by buck:
Wolfram gives the value of the same as 4.7e21, which seems correct.
http://www.wolframalpha.com/input/?i=Gamma%5B25%2C+0%2C+14.5%5D
by buck:
Actually, if I ask pari directly, it gets the correct value, so maybe the
algorithm argument is broken?
{{{
sage: pari('gamma(25) - incgam(25, 14.5)')
4.71197173824555 E21
sage: pari('incgamc(25, 14.5)')
4.71197173824555 E21
}}}
by buck:
Same story if I ask mpmath directly. I can't imagine what's wrong.
{{{
sage: call(mpmath.gammainc, 25, 0, 14.5)
4.71197173824555e21
}}}
by rws:
Replying to buck:
> Same story if I ask mpmath directly. I can't imagine what's wrong.
'algorithm' isn't propagated to evalf so we always get the buggy pari
values. As we would have to switch the default anyway as soon as a later
pari version fixes the bug, I'll now simply disable pari. We'll have to
live with mpmath's 53 bits for now.
by rws:
Replying to buck:
> Actually, if I ask pari directly, it gets the correct value, so maybe
the algorithm argument is broken?
> {{{
> sage: pari('gamma(25) - incgam(25, 14.5)')
> 4.71197173824555 E21
> sage: pari('incgamc(25, 14.5)')
> 4.71197173824555 E21
> }}}
This works because it uses the gp interface (the libpari might be broken
on our side). Ah, that's two different bugs. However, even if it is in our
interface to libpari and we fix this, or if we use the gp interface, there
is still the unfixed pari bug at
http://pari.math.u-bordeaux.fr/cgi-bin/bugreport.cgi?bug=1689PARI/GP
by rws:
Actually the 53-bit restriction applies only to larger values, so all is
not lost.
--
Ticket URL: <http://trac.sagemath.org/ticket/18210#comment:21>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-trac.
For more options, visit https://groups.google.com/d/optout.
