Opened 5 years ago

Closed 5 years ago

Last modified 5 years ago

#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) 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...

See possibly also here and here and possibly even #16697.

Change History (13)

comment:1 Changed 5 years ago by rws

  • Cc rws added

comment:2 Changed 5 years ago by kcrisman

Possibly - even likely - related is #11164.

comment:3 Changed 5 years ago by pbruin

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: Changed 5 years ago by 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.

comment:5 in reply to: ↑ 4 Changed 5 years ago by pbruin

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 5 years ago by pbruin

  • 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 5 years ago by pbruin

  • Authors set to Peter Bruin
  • 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: Changed 5 years ago by 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.

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 5 years ago by kcrisman

Passes relevant tests, so other than that question we're good to go.

comment:10 in reply to: ↑ 8 Changed 5 years ago by pbruin

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 5 years ago by kcrisman

  • Reviewers set to Karl-Dieter Crisman
  • Status changed from needs_review to positive_review

comment:12 Changed 5 years ago by vbraun

  • Branch changed from u/pbruin/17328-incgam to 3cae7c8469d6c2f5dd3d15a872965f528f637c5c
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:13 Changed 5 years ago by jdemeyer

  • Commit 3cae7c8469d6c2f5dd3d15a872965f528f637c5c deleted

Just to let you know: I already had a patch for this at #17130, so the code of this patch will be overwritten again in #17130. I noticed this ticket only now because of a merge conflict.

Note: See TracTickets for help on using tickets.