Opened 8 years ago

Closed 8 years ago

# incomplete gamma function bugs for certain arguments

Reported by: Owned by: Karl-Dieter Crisman major sage-6.4 calculus incomplete gamma function Ralf Stephan Peter Bruin Karl-Dieter Crisman N/A 3cae7c8

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

### comment:2 Changed 8 years ago by Karl-Dieter Crisman

Possibly - even likely - related is #11164.

### comment:3 Changed 8 years ago by Peter Bruin

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 Karl-Dieter Crisman

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 8 years ago by Peter Bruin

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 Peter Bruin

Branch: → u/pbruin/17328-incgam → 3cae7c8469d6c2f5dd3d15a872965f528f637c5c critical → major new → needs_review

### comment:7 in reply to:  4 Changed 8 years ago by Peter Bruin

Authors: → Peter Bruin incomplete gamma function added incomplete gamma in integral is wrong → incomplete gamma function bugs for certain arguments

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 Karl-Dieter Crisman

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.

```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 Karl-Dieter Crisman

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

### comment:10 in reply to:  8 Changed 8 years ago by Peter Bruin

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.

```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 Karl-Dieter Crisman

Reviewers: → Karl-Dieter Crisman needs_review → positive_review

### comment:12 Changed 8 years ago by Volker Braun

Branch: u/pbruin/17328-incgam → 3cae7c8469d6c2f5dd3d15a872965f528f637c5c → fixed positive_review → closed

### comment:13 Changed 8 years ago by Jeroen Demeyer

Commit: 3cae7c8469d6c2f5dd3d15a872965f528f637c5c

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.