#17328 closed defect (fixed)
incomplete gamma function bugs for certain arguments
Reported by: | Karl-Dieter Crisman | Owned by: | |
---|---|---|---|
Priority: | major | Milestone: | sage-6.4 |
Component: | calculus | Keywords: | incomplete gamma function |
Cc: | Ralf Stephan | 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 8 years ago by
Cc: | Ralf Stephan added |
---|
comment:2 Changed 8 years ago by
comment:3 Changed 8 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 8 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 Changed 8 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 8 years ago by
Branch: | → u/pbruin/17328-incgam |
---|---|
Commit: | → 3cae7c8469d6c2f5dd3d15a872965f528f637c5c |
Priority: | critical → major |
Status: | new → needs_review |
comment:7 Changed 8 years ago by
Authors: | → Peter Bruin |
---|---|
Keywords: | incomplete gamma function added |
Summary: | incomplete gamma in integral is wrong → 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 8 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 8 years ago by
Passes relevant tests, so other than that question we're good to go.
comment:10 Changed 8 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 8 years ago by
Reviewers: | → Karl-Dieter Crisman |
---|---|
Status: | needs_review → positive_review |
comment:12 Changed 8 years ago by
Branch: | u/pbruin/17328-incgam → 3cae7c8469d6c2f5dd3d15a872965f528f637c5c |
---|---|
Resolution: | → fixed |
Status: | positive_review → closed |
comment:13 Changed 8 years ago by
Commit: | 3cae7c8469d6c2f5dd3d15a872965f528f637c5c |
---|
Possibly - even likely - related is #11164.