incomplete gamma function bugs for certain arguments

[kcrisman]

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.

CC: @rwst

Component: calculus

Keywords: incomplete gamma function

Author: Peter Bruin

Branch: 3cae7c8

Reviewer: Karl-Dieter Crisman

Issue created by migration from https://trac.sagemath.org/ticket/17328

[kcrisman]

comment:2

Possibly - even likely - related is #11164.

[pjbruin]

comment:3

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

[kcrisman]

comment:4

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.

[pjbruin]

comment:5

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.

[pjbruin]

Commit: 3cae7c8

[pjbruin]

Branch: u/pbruin/17328-incgam

[pjbruin]

comment:7

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.

[pjbruin]

Author: Peter Bruin

[pjbruin]

Changed keywords from none to incomplete gamma function

[kcrisman]

comment:8

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?

[kcrisman]

comment:9

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

[pjbruin]

comment:10

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".

[kcrisman]

Reviewer: Karl-Dieter Crisman

[vbraun]

Changed branch from u/pbruin/17328-incgam to 3cae7c8

[jdemeyer]

Changed commit from 3cae7c8 to none

[jdemeyer]

comment:13

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.