#17328 closed defect (fixed)
incomplete gamma function bugs for certain arguments
Reported by: | kcrisman | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-6.4 |
Component: | calculus | Keywords: | incomplete gamma function |
Cc: | rws | Merged in: | |
Authors: | Peter Bruin | Reviewers: | Karl-Dieter Crisman |
Report Upstream: | N/A | Work issues: | |
Branch: | 3cae7c8 (Commits, GitHub, GitLab) | Commit: | |
Dependencies: | Stopgaps: |
Description
See this sage-support thread for details.
sage: N(integral(1/log(x)^2,(x,2,3))) 0.536566859259958 # quite wrong sage: integral(1/(ln(x))^2, x,2,3) gamma(-1, -log(3)) - gamma(-1, -log(2))
We get the antiderivative and answer from Maxima, which evaluates this correctly numerically.
(%i1) display2d : false $ (%i2) foo : 1/log(x)^2 $ (%i3) integrate (foo, x, 2, 3); (%o3) gamma_incomplete(-1,-log(3))-gamma_incomplete(-1,-log(2)) (%i4) %, numer; (%o4) 1.273097216447114 (%i5) integrate (foo, x); (%o5) gamma_incomplete(-1,-log(x)) (%i6) ev (%, x=3) - ev (%, x=2); (%o6) gamma_incomplete(-1,-log(3))-gamma_incomplete(-1,-log(2)) (%i7) %, numer; (%o7) 1.273097216447114
See in particular this ugly plot.
sage: plot(lambda t: numerical_integral(1/ln(x)^2,2,t)[0],2,3)+plot(lambda t: gamma(-1, -log(t)).real(),2,3,color='red')
And even that we have to call real()
to remove the numerical noise is not good...
Change History (13)
comment:1 Changed 7 years ago by
- Cc rws added
comment:2 Changed 7 years ago by
comment:3 Changed 7 years ago by
I can't seem to reproduce this in Sage 6.4:
sage: numerical_integral(1/log(x)^2,2,3) (1.273097216447114, 1.4134218422857824e-14) sage: N(integral(1/log(x)^2,(x,2,3))) 1.27309721644711
In fact, the incomplete gamma function is evaluated using PARI, and the bug appears to have been fixed by the latest PARI upgrade (see #15767):
gp-2.5 > incgam(-1,-1) %1 = -1.5817775437849803404714887501837998073 gp-2.5 > incgam(-1,-log(3)) - incgam(-1,-log(2)) %2 = 0.54869101686030432559293241097969276965 gp-2.7 > incgam(-1,-1) %1 = -0.82316401210310847989376653702102822882 + 3.1415926535897932384626433832795028842*I gp-2.7 > incgam(-1,-log(3)) - incgam(-1,-log(2)) %2 = 1.2730972164471138219094623429485715030 + 0.E-37*I
comment:4 follow-ups: ↓ 5 ↓ 7 Changed 7 years ago by
Interesting and confirmed. Of course, possibly due to something completely different, we now have
sage: plot(lambda t: gamma(-1, -log(t)),2,3,color='red') AttributeError: type object 'float' has no attribute 'precision'
as opposed to a wacky plot, no plot.
comment:5 in reply to: ↑ 4 Changed 7 years ago by
Replying to kcrisman:
Interesting and confirmed. Of course, possibly due to something completely different, we now have
sage: plot(lambda t: gamma(-1, -log(t)),2,3,color='red') AttributeError: type object 'float' has no attribute 'precision'as opposed to a wacky plot, no plot.
That is because of a patch I made at #7099; it assumed that the parent had a precision()
method (previously the precision was ignored when converting the arguments into a ComplexField
). I will fix this as soon as I have time.
We should also add a doctest for the original bug.
comment:6 Changed 7 years ago by
- Branch set to u/pbruin/17328-incgam
- Commit set to 3cae7c8469d6c2f5dd3d15a872965f528f637c5c
- Priority changed from critical to major
- Status changed from new to needs_review
comment:7 in reply to: ↑ 4 Changed 7 years ago by
- Keywords incomplete gamma function added
- Summary changed from incomplete gamma in integral is wrong to incomplete gamma function bugs for certain arguments
Replying to kcrisman:
Interesting and confirmed. Of course, possibly due to something completely different, we now have
sage: plot(lambda t: gamma(-1, -log(t)),2,3,color='red') AttributeError: type object 'float' has no attribute 'precision'as opposed to a wacky plot, no plot.
It would be good if you (or someone else) could try plotting this with the above patch and check if it looks good to you.
comment:8 follow-up: ↓ 10 Changed 7 years ago by
Seems to give correct answers in various precisions.
plot(lambda t: gamma(-1, -log(t)).real() - numerical_integral(1/ln(x)^2,2,t)[0],2,3,color='red')
is a very nice straight line just below zero now. And the imaginary part seems to be consistently pi for t>1. And the fix is a good one.
We should also add a doctest for the original bug.
Did you want to add
N(integral(1/log(x)^2,(x,2,3)))
as well then, or is the example you have with the actual antideriv and floats sufficient?
comment:9 Changed 7 years ago by
Passes relevant tests, so other than that question we're good to go.
comment:10 in reply to: ↑ 8 Changed 7 years ago by
Replying to kcrisman:
Seems to give correct answers in various precisions.
plot(lambda t: gamma(-1, -log(t)).real() - numerical_integral(1/ln(x)^2,2,t)[0],2,3,color='red')is a very nice straight line just below zero now. And the imaginary part seems to be consistently pi for t>1. And the fix is a good one.
OK, thanks!
We should also add a doctest for the original bug.
Did you want to add
N(integral(1/log(x)^2,(x,2,3)))as well then, or is the example you have with the actual antideriv and floats sufficient?
I don't think it is necessary to add another doctest for the integral; this would just test in addition that Maxima computes the correct antiderivative. Of course the above command was how the bug was originally found, but I don't really see why we should doctest exactly that command. What I meant to say is "add a doctest to show that the PARI bug has been fixed".
comment:11 Changed 7 years ago by
- Reviewers set to Karl-Dieter Crisman
- Status changed from needs_review to positive_review
comment:12 Changed 7 years ago by
- Branch changed from u/pbruin/17328-incgam to 3cae7c8469d6c2f5dd3d15a872965f528f637c5c
- Resolution set to fixed
- Status changed from positive_review to closed
comment:13 Changed 7 years ago by
- Commit 3cae7c8469d6c2f5dd3d15a872965f528f637c5c deleted
Possibly - even likely - related is #11164.