Opened 7 years ago

Closed 7 years ago

# incomplete gamma function bugs for certain arguments

Reported by: Owned by: kcrisman major sage-6.4 calculus incomplete gamma function rws 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),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 7 years ago by kcrisman

Possibly - even likely - related is #11164.

### comment:3 Changed 7 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: ↓ 5 ↓ 7 Changed 7 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 7 years ago by pbruin

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

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

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

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

Seems to give correct answers in various precisions.

```plot(lambda t: gamma(-1, -log(t)).real() - numerical_integral(1/ln(x)^2,2,t),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 7 years ago by kcrisman

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

### comment:12 Changed 7 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 7 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.