Sage: Ticket #17328: incomplete gamma function bugs for certain arguments
https://trac.sagemath.org/ticket/17328
<p>
See <a class="ext-link" href="https://groups.google.com/forum/#!topic/sage-support/lGGO_q_NVzg"><span class="icon"></span>this sage-support thread</a> for details.
</p>
<pre class="wiki">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))
</pre><p>
We get the antiderivative and answer from Maxima, which evaluates this correctly numerically.
</p>
<pre class="wiki"> (%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
</pre><p>
See in particular this ugly plot.
</p>
<pre class="wiki">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')
</pre><p>
And even that we have to call <code>real()</code> to remove the numerical noise is not good...
</p>
<p>
See possibly also <a class="ext-link" href="https://groups.google.com/forum/#!topic/sage-support/RrveCpPgoKU"><span class="icon"></span>here</a> and <a class="ext-link" href="https://groups.google.com/forum/#!topic/sage-support/RD3YH3jb3qo"><span class="icon"></span>here</a> and possibly even <a class="closed ticket" href="https://trac.sagemath.org/ticket/16697" title="defect: implement symbolic lower incomplete gamma function (closed: fixed)">#16697</a>.
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/17328
Trac 1.1.6rwsFri, 14 Nov 2014 08:47:48 GMTcc set
https://trac.sagemath.org/ticket/17328#comment:1
https://trac.sagemath.org/ticket/17328#comment:1
<ul>
<li><strong>cc</strong>
<em>rws</em> added
</li>
</ul>
TicketkcrismanFri, 14 Nov 2014 21:08:37 GMT
https://trac.sagemath.org/ticket/17328#comment:2
https://trac.sagemath.org/ticket/17328#comment:2
<p>
Possibly - even likely - related is <a class="needs_info ticket" href="https://trac.sagemath.org/ticket/11164" title="defect: Integral of sin(x)/x gives false result. (needs_info)">#11164</a>.
</p>
TicketpbruinFri, 14 Nov 2014 22:08:36 GMT
https://trac.sagemath.org/ticket/17328#comment:3
https://trac.sagemath.org/ticket/17328#comment:3
<p>
I can't seem to reproduce this in Sage 6.4:
</p>
<pre class="wiki">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
</pre><p>
In fact, the incomplete gamma function is evaluated using PARI, and the bug appears to have been fixed by the latest PARI upgrade (see <a class="closed ticket" href="https://trac.sagemath.org/ticket/15767" title="enhancement: Upgrade PARI to 2.7.1 (closed: fixed)">#15767</a>):
</p>
<pre class="wiki">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
</pre>
TicketkcrismanSat, 15 Nov 2014 02:58:05 GMT
https://trac.sagemath.org/ticket/17328#comment:4
https://trac.sagemath.org/ticket/17328#comment:4
<p>
Interesting and confirmed. Of course, possibly due to something completely different, we now have
</p>
<pre class="wiki">sage: plot(lambda t: gamma(-1, -log(t)),2,3,color='red')
AttributeError: type object 'float' has no attribute 'precision'
</pre><p>
as opposed to a wacky plot, no plot.
</p>
TicketpbruinSat, 15 Nov 2014 11:24:23 GMT
https://trac.sagemath.org/ticket/17328#comment:5
https://trac.sagemath.org/ticket/17328#comment:5
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/17328#comment:4" title="Comment 4">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
Interesting and confirmed. Of course, possibly due to something completely different, we now have
</p>
<pre class="wiki">sage: plot(lambda t: gamma(-1, -log(t)),2,3,color='red')
AttributeError: type object 'float' has no attribute 'precision'
</pre><p>
as opposed to a wacky plot, no plot.
</p>
</blockquote>
<p>
That is because of a patch I made at <a class="closed ticket" href="https://trac.sagemath.org/ticket/7099" title="defect: serious incomplete gamma function precision bugs (closed: fixed)">#7099</a>; it assumed that the parent had a <code>precision()</code> method (previously the precision was ignored when converting the arguments into a <code>ComplexField</code>). I will fix this as soon as I have time.
</p>
<p>
We should also add a doctest for the original bug.
</p>
TicketpbruinMon, 17 Nov 2014 12:00:58 GMTstatus, priority changed; commit, branch set
https://trac.sagemath.org/ticket/17328#comment:6
https://trac.sagemath.org/ticket/17328#comment:6
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
<li><strong>commit</strong>
set to <em>3cae7c8469d6c2f5dd3d15a872965f528f637c5c</em>
</li>
<li><strong>branch</strong>
set to <em>u/pbruin/17328-incgam</em>
</li>
<li><strong>priority</strong>
changed from <em>critical</em> to <em>major</em>
</li>
</ul>
TicketpbruinMon, 17 Nov 2014 12:03:45 GMTsummary changed; keywords, author set
https://trac.sagemath.org/ticket/17328#comment:7
https://trac.sagemath.org/ticket/17328#comment:7
<ul>
<li><strong>keywords</strong>
<em>incomplete</em> <em>gamma</em> <em>function</em> added
</li>
<li><strong>summary</strong>
changed from <em>incomplete gamma in integral is wrong</em> to <em>incomplete gamma function bugs for certain arguments</em>
</li>
<li><strong>author</strong>
set to <em>Peter Bruin</em>
</li>
</ul>
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/17328#comment:4" title="Comment 4">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
Interesting and confirmed. Of course, possibly due to something completely different, we now have
</p>
<pre class="wiki">sage: plot(lambda t: gamma(-1, -log(t)),2,3,color='red')
AttributeError: type object 'float' has no attribute 'precision'
</pre><p>
as opposed to a wacky plot, no plot.
</p>
</blockquote>
<p>
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.
</p>
TicketkcrismanMon, 17 Nov 2014 14:54:41 GMT
https://trac.sagemath.org/ticket/17328#comment:8
https://trac.sagemath.org/ticket/17328#comment:8
<p>
Seems to give correct answers in various precisions.
</p>
<pre class="wiki">plot(lambda t: gamma(-1, -log(t)).real() - numerical_integral(1/ln(x)^2,2,t)[0],2,3,color='red')
</pre><p>
is a very nice straight line <em>just</em> below zero now. And the imaginary part seems to be consistently pi for t>1. And the fix is a good one.
</p>
<blockquote class="citation">
<p>
We should also add a doctest for the original bug.
</p>
</blockquote>
<p>
Did you want to add
</p>
<pre class="wiki">N(integral(1/log(x)^2,(x,2,3)))
</pre><p>
as well then, or is the example you have with the actual antideriv and floats sufficient?
</p>
TicketkcrismanMon, 17 Nov 2014 16:07:06 GMT
https://trac.sagemath.org/ticket/17328#comment:9
https://trac.sagemath.org/ticket/17328#comment:9
<p>
Passes relevant tests, so other than that question we're good to go.
</p>
TicketpbruinMon, 17 Nov 2014 19:27:44 GMT
https://trac.sagemath.org/ticket/17328#comment:10
https://trac.sagemath.org/ticket/17328#comment:10
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/17328#comment:8" title="Comment 8">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
Seems to give correct answers in various precisions.
</p>
<pre class="wiki">plot(lambda t: gamma(-1, -log(t)).real() - numerical_integral(1/ln(x)^2,2,t)[0],2,3,color='red')
</pre><p>
is a very nice straight line <em>just</em> below zero now. And the imaginary part seems to be consistently pi for t>1. And the fix is a good one.
</p>
</blockquote>
<p>
OK, thanks!
</p>
<blockquote class="citation">
<blockquote class="citation">
<p>
We should also add a doctest for the original bug.
</p>
</blockquote>
<p>
Did you want to add
</p>
<pre class="wiki">N(integral(1/log(x)^2,(x,2,3)))
</pre><p>
as well then, or is the example you have with the actual antideriv and floats sufficient?
</p>
</blockquote>
<p>
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".
</p>
TicketkcrismanMon, 17 Nov 2014 21:54:40 GMTstatus changed; reviewer set
https://trac.sagemath.org/ticket/17328#comment:11
https://trac.sagemath.org/ticket/17328#comment:11
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
<li><strong>reviewer</strong>
set to <em>Karl-Dieter Crisman</em>
</li>
</ul>
TicketvbraunWed, 19 Nov 2014 08:32:10 GMTstatus, branch changed; resolution set
https://trac.sagemath.org/ticket/17328#comment:12
https://trac.sagemath.org/ticket/17328#comment:12
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>closed</em>
</li>
<li><strong>resolution</strong>
set to <em>fixed</em>
</li>
<li><strong>branch</strong>
changed from <em>u/pbruin/17328-incgam</em> to <em>3cae7c8469d6c2f5dd3d15a872965f528f637c5c</em>
</li>
</ul>
TicketjdemeyerWed, 26 Nov 2014 08:18:51 GMTcommit deleted
https://trac.sagemath.org/ticket/17328#comment:13
https://trac.sagemath.org/ticket/17328#comment:13
<ul>
<li><strong>commit</strong>
<em>3cae7c8469d6c2f5dd3d15a872965f528f637c5c</em> deleted
</li>
</ul>
<p>
Just to let you know: I already had a patch for this at <a class="closed ticket" href="https://trac.sagemath.org/ticket/17130" title="defect: Fix coercion bugs in symbolic functions (closed: fixed)">#17130</a>, so the code of this patch will be overwritten again in <a class="closed ticket" href="https://trac.sagemath.org/ticket/17130" title="defect: Fix coercion bugs in symbolic functions (closed: fixed)">#17130</a>. I noticed this ticket only now because of a merge conflict.
</p>
Ticket