Sage: Ticket #11164: Integral of sin(x)/x gives false result.
https://trac.sagemath.org/ticket/11164
<pre class="wiki">> eq = sin(x)/x
> integrate(eq, x, -1e-6, 1e-6)
-3.14159065359
</pre><p>
I expected an answer of about 2e-6.
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/11164
Trac 1.2Karl-Dieter CrismanThu, 14 Apr 2011 01:51:49 GMTpriority, owner, component, upstream changed
https://trac.sagemath.org/ticket/11164#comment:1
https://trac.sagemath.org/ticket/11164#comment:1
<ul>
<li><strong>priority</strong>
changed from <em>minor</em> to <em>major</em>
</li>
<li><strong>owner</strong>
changed from <em>tbd</em> to <em>Burcin Erocal</em>
</li>
<li><strong>component</strong>
changed from <em>PLEASE CHANGE</em> to <em>calculus</em>
</li>
<li><strong>upstream</strong>
changed from <em>N/A</em> to <em>Fixed upstream, in a later stable release.</em>
</li>
</ul>
<p>
In the future, be sure to pick a component (for instance, calculus or symbolics); that will help people find it more easily.
</p>
<p>
There are two problems here.
</p>
<ul><li>The actual problem - the report unfortunately has a typo.
<pre class="wiki">sage: integrate(eq,x,-1e-6,1e-6)
-3.14159065359
</pre>This is definitely a problem. However, it's fixed in the latest Maxima. Here is the old one, in Sage:
<pre class="wiki">Maxima 5.22.1 http://maxima.sourceforge.net
using Lisp ECL 10.4.1
<snip>
(%i2) keepfloat:true;
(%o2) true
(%i3) integrate(sin(x)/x,x,-1.00000000000000e-6,1.00000000000000e-6)
;
(%o3) - 3.141590653589793
</pre>And here is 5.24.0. I assume this isn't a problem of which Lisp we are using?
<pre class="wiki">Maxima 5.24.0 http://maxima.sourceforge.net
using Lisp SBCL 1.0.24
(%i3) keepfloat:true;
(%o3) true
(%i8) integrate(sin(x)/x,x,-1.00000000000000e-6,1.00000000000000e-6);
rat: replaced 9.999999999999999e-7 by 1/1000000 = 9.999999999999999e-7
rat: replaced -9.999999999999999e-7 by -1/1000000 = -9.999999999999999e-7
rat: replaced -9.999999999999999e-7 by -1/1000000 = -9.999999999999999e-7
rat: replaced 9.999999999999999e-7 by 1/1000000 = 9.999999999999999e-7
%i %i
(%o8) %i gamma_incomplete(0, - -------) - %i gamma_incomplete(0, -------)
1000000 1000000
</pre></li><li>A problem revealed, but also fixed in the latest. Old:
<pre class="wiki">(%i4) integrate(sin(x)/x,x,-1e-6,1e6);
Maxima encountered a Lisp error:
#<a FLOATING-POINT-OVERFLOW>
Automatically continuing.
To enable the Lisp debugger set *debugger-hook* to nil.
</pre>New:
<pre class="wiki">(%i10) integrate(sin(x)/x,x,-1e-6,1e6);
rat: replaced 9.999999999999999e-7 by 1/1000000 = 9.999999999999999e-7
rat: replaced -9.999999999999999e-7 by -1/1000000 = -9.999999999999999e-7
rat: replaced -1000000.0 by -1000000/1 = -1000000.0
rat: replaced 1000000.0 by 1000000/1 = 1000000.0
%i
%i gamma_incomplete(0, -------)
%i gamma_incomplete(0, 1000000 %i) 1000000
(%o10) - ---------------------------------- - -------------------------------
2 2
%i
%i gamma_incomplete(0, - -------)
1000000 %i gamma_incomplete(0, - 1000000 %i)
+ --------------------------------- + ------------------------------------
2 2
</pre></li></ul><p>
We'd need doctesting to test both of these once we upgrade Maxima.
</p>
TicketD.S. McNeilThu, 15 Dec 2011 05:56:42 GMTattachment set
https://trac.sagemath.org/ticket/11164
https://trac.sagemath.org/ticket/11164
<ul>
<li><strong>attachment</strong>
set to <em>trac_11164_verify_integration_sinx_x.patch</em>
</li>
</ul>
<p>
add doctest to confirm fix
</p>
TicketD.S. McNeilThu, 15 Dec 2011 05:57:52 GMTstatus changed
https://trac.sagemath.org/ticket/11164#comment:2
https://trac.sagemath.org/ticket/11164#comment:2
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
</ul>
<p>
Might as well doctest and close.
</p>
TicketKarl-Dieter CrismanThu, 15 Dec 2011 14:17:41 GMTstatus changed
https://trac.sagemath.org/ticket/11164#comment:3
https://trac.sagemath.org/ticket/11164#comment:3
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
<p>
I may have misstated whether this is "fixed".
</p>
<pre class="wiki">sage: integrate(sin(x)/x,x,-1e-6,1e-6)
I*Ei(-1/1000000*I) - I*Ei(1/1000000*I)
sage: integrate(sin(x)/x,x,-1e-6,1e-6).n()
3.14159465358979
sage: numerical_integral(sin(x)/x,-1e-6,1e-6)
(1.9999999999998889e-06, nan)
</pre><p>
The graph makes it very obvious the second (numerical) answer is correct. Unsurprisingly, the other answer is off by exactly <code>pi</code>...
</p>
<p>
Presumably the first answer is "correct" in the sense that it is an alternate representation of the <code>Si</code> sine integral function. But I think there are some branch issues - <code>Ei(I*x)</code> only has an alternate formula for <code>x>0</code>.
</p>
<p>
So we still need to do something with this. Or at the very least warn about branches. But really, <code>sin(x)/x</code> shouldn't have to worry about this!
</p>
TicketD.S. McNeilThu, 15 Dec 2011 14:30:41 GMTstatus changed
https://trac.sagemath.org/ticket/11164#comment:4
https://trac.sagemath.org/ticket/11164#comment:4
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_info</em>
</li>
</ul>
<p>
You're right, of course. It's silly to start with sin(x)/x and then generalize to something which doesn't reduce to the expected value, but far, far sillier to add a doctest enshrining the goofy behaviour as the expected. :<sup>)
</sup></p>
<p>
Does this work on the current maxima trunk now?
</p>
TicketKarl-Dieter CrismanThu, 15 Dec 2011 14:53:18 GMT
https://trac.sagemath.org/ticket/11164#comment:5
https://trac.sagemath.org/ticket/11164#comment:5
<p>
This is the previous version (5.25.1):
</p>
<pre class="wiki">
(%i5) integrate(sin(x)/x,x,-1/1000000,1/1000000);
(%o5) %i*gamma_incomplete(0,-%i/1000000)-%i*gamma_incomplete(0,%i/1000000)
</pre><p>
W|A agrees that this is about <code>-pi</code>. (To be precise, <code>2*Si(1/10^6)-pi</code>.)
</p>
<pre class="wiki">
(%i4) display2d:false;
(%o4) false
(%i5) integrate(sin(x)/x,x);
(%o5) -(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2
</pre><p>
Here is the issue; this function may or may not be the same as the sine integral. The problem seems to be that although Maxima has the sine integral (<code>expintegral_si</code>), it isn't really connected with the rest of Maxima (?).
</p>
<p>
What do you think a good solution to this is? I suppose we could open a ticket for it on the Maxima site, but I'm not sure whether the <code>Si</code> there is robust enough to handle all this yet.
</p>
TicketKarl-Dieter CrismanThu, 15 Dec 2011 14:54:44 GMT
https://trac.sagemath.org/ticket/11164#comment:6
https://trac.sagemath.org/ticket/11164#comment:6
<blockquote class="citation">
<p>
Here is the issue; this function may or may not be the same as the sine integral.
</p>
</blockquote>
<p>
By which I mean it isn't the same, of course, but also that (as defined) it isn't even always the same constant away. Antiderivatives can differ by a constant, but when the constant is different depending on what <code>x</code> is, that's annoying. I think the Maxima folks would say this is a feature (for symbolic reasons).
</p>
TicketKarl-Dieter CrismanSat, 24 Dec 2011 04:45:32 GMT
https://trac.sagemath.org/ticket/11164#comment:7
https://trac.sagemath.org/ticket/11164#comment:7
<p>
By the way, achrzesz has a great one-liner for <a class="ext-link" href="http://ask.sagemath.org/question/6/how-do-i-understand-the-result-of-symbolic"><span class="icon"></span>a related question</a>:
</p>
<pre class="wiki"> maxima('integrate(diff((exp(x)-1)/x,x),x),gamma_expand:true,factor').sage()
</pre><p>
Though here it doesn't solve the problem, but maybe it will trigger someone's memory as to another Maxima flag to set...
</p>
<pre class="wiki">sage: maxima('integrate(sin(x)/x,x,-1/1000000,1/1000000),gamma_expand:true').sage()
I*expintegral_ei(-1/1000000*I) - I*expintegral_ei(1/1000000*I) # currently not translated into Sage, but see #11143
sage: a = I*Ei(-1/1000000*I)-I*Ei(1/1000000*I)
sage: a.n()
3.14159465358979
</pre>
TicketKarl-Dieter CrismanSun, 20 Jan 2013 01:38:52 GMT
https://trac.sagemath.org/ticket/11164#comment:8
https://trac.sagemath.org/ticket/11164#comment:8
<p>
I wonder whether this is potentially fixed in <a class="closed ticket" href="https://trac.sagemath.org/ticket/13364" title="#13364: enhancement: Upgrade Maxima to 5.29.1 (closed: fixed)">#13364</a> - I don't know if we ever reported upstream, or whether we should have...
</p>
TicketKarl-Dieter CrismanWed, 26 Jun 2013 14:48:17 GMT
https://trac.sagemath.org/ticket/11164#comment:9
https://trac.sagemath.org/ticket/11164#comment:9
<blockquote class="citation">
<p>
I wonder whether this is potentially fixed in <a class="closed ticket" href="https://trac.sagemath.org/ticket/13364" title="#13364: enhancement: Upgrade Maxima to 5.29.1 (closed: fixed)">#13364</a> - I don't know if we ever reported upstream, or whether we should have...
</p>
</blockquote>
<p>
Apparently not yet.
</p>
TicketTravis ScrimshawMon, 23 Jun 2014 18:31:03 GMTdescription changed
https://trac.sagemath.org/ticket/11164#comment:10
https://trac.sagemath.org/ticket/11164#comment:10
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/11164?action=diff&version=10">diff</a>)
</li>
</ul>
<p>
In 6.3.beta4 (so with <a class="closed ticket" href="https://trac.sagemath.org/ticket/13973" title="#13973: defect: Upgrade Maxima to 5.33.0 (closed: fixed)">#13973</a>), I get:
</p>
<pre class="wiki">sage: eq = sin(x) / x
sage: integrate(eq, x, -1e-6, 1e-6)
-I*Ei(1/1000000*I) + I*Ei(-1/1000000*I)
sage: _.n()
3.14159465358979
sage: integrate(eq, x, -1e-12, 1e-12)
-I*Ei(1/1000000000000*I) + I*Ei(-1/1000000000000*I)
sage: _.n()
3.14159265359179
</pre><p>
Something bad happens with the integration computation when crossing over 0 as we have:
</p>
<pre class="wiki">sage: integrate(eq, x, 1e-12, 1e-6)
-1/2*I*Ei(1/1000000*I) + 1/2*I*Ei(1/1000000000000*I) - 1/2*I*Ei(-1/1000000000000*I) + 1/2*I*Ei(-1/1000000*I)
sage: _.n()
9.99999000050877e-7
</pre>
TicketPeter BruinTue, 24 Jun 2014 00:37:52 GMT
https://trac.sagemath.org/ticket/11164#comment:11
https://trac.sagemath.org/ticket/11164#comment:11
<p>
The problem is caused by the integral being computed via symbolic integration; see also <a class="ticket" href="https://trac.sagemath.org/ticket/11164#comment:5" title="Comment 5">comment:5</a>. Maxima integrates <code>exp(x)/x</code> to <code>-gamma_incomplete(0, -x)</code>, which is in principle correct, but only gives a certain branch of the integral. The branch cut chosen for the incomplete Gamma function is along the negative real axis. When using this to compute <code>integrate(sin(x)/x, x, -1, 1)</code>, the two terms coming from <code>exp(%i*x)</code> and <code>exp(-%i*x)</code> have branch cuts along the negative and positive imaginary axis, respectively, so the sum has a branch cut along the whole imaginary axis and therefore can't be the antiderivative of <code>sin(x)/x</code>, which is an entire function.
</p>
TicketKarl-Dieter CrismanFri, 14 Nov 2014 21:08:25 GMT
https://trac.sagemath.org/ticket/11164#comment:12
https://trac.sagemath.org/ticket/11164#comment:12
<blockquote class="citation">
<p>
The problem is caused by the integral being computed via symbolic integration; see also <a class="ticket" href="https://trac.sagemath.org/ticket/11164#comment:5" title="Comment 5">comment:5</a>. Maxima integrates <code>exp(x)/x</code> to <code>-gamma_incomplete(0, -x)</code>, which is in principle correct, but only gives a certain branch of the integral. The branch cut chosen for the incomplete Gamma function is along the negative real axis. When using this to compute <code>integrate(sin(x)/x, x, -1, 1)</code>, the two terms coming from <code>exp(%i*x)</code> and <code>exp(-%i*x)</code> have branch cuts along the negative and positive imaginary axis, respectively, so the sum has a branch cut along the whole imaginary axis and therefore can't be the antiderivative of <code>sin(x)/x</code>, which is an entire function.
</p>
</blockquote>
<p>
How much you wanna bet this is the same underlying problem as in <a class="closed ticket" href="https://trac.sagemath.org/ticket/17328" title="#17328: defect: incomplete gamma function bugs for certain arguments (closed: fixed)">#17328</a>?
</p>
TicketKarl-Dieter CrismanSat, 15 Nov 2014 02:58:00 GMT
https://trac.sagemath.org/ticket/11164#comment:13
https://trac.sagemath.org/ticket/11164#comment:13
<blockquote class="citation">
<p>
How much you wanna bet this is the same underlying problem as in <a class="closed ticket" href="https://trac.sagemath.org/ticket/17328" title="#17328: defect: incomplete gamma function bugs for certain arguments (closed: fixed)">#17328</a>?
</p>
</blockquote>
<p>
Nah, I was thrown off by the gammas, I think.
</p>
<p>
What if we don't allow incomplete gamma to become <code>Ei</code> when the first arg is zero - would that solve the branch problem?
</p>
TicketKarl-Dieter CrismanMon, 17 Nov 2014 16:19:24 GMT
https://trac.sagemath.org/ticket/11164#comment:14
https://trac.sagemath.org/ticket/11164#comment:14
<blockquote class="citation">
<p>
What if we don't allow incomplete gamma to become <code>Ei</code> when the first arg is zero - would that solve the branch problem?
</p>
</blockquote>
<p>
Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).
</p>
TicketKarl-Dieter CrismanMon, 17 Nov 2014 16:39:01 GMTupstream changed
https://trac.sagemath.org/ticket/11164#comment:15
https://trac.sagemath.org/ticket/11164#comment:15
<ul>
<li><strong>upstream</strong>
changed from <em>Fixed upstream, in a later stable release.</em> to <em>Reported upstream. No feedback yet.</em>
</li>
</ul>
<blockquote class="citation">
<blockquote class="citation">
<p>
What if we don't allow incomplete gamma to become <code>Ei</code> when the first arg is zero - would that solve the branch problem?
</p>
</blockquote>
<p>
Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).
</p>
</blockquote>
<p>
Or maybe not;
</p>
<pre class="wiki">sage: maxima('integrate(sin(x)/x,x)')
-(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2
sage: maxima('integrate(sin(x)/x,x), gamma_expand:true')
-(%i*expintegral_ei(%i*x)-%i*expintegral_ei(-%i*x)+2*%pi)/2
</pre><p>
so maybe we are already setting <code>gamma_expand:true</code> somewhere nowadays. Reported upstream at <a class="ext-link" href="https://sourceforge.net/p/maxima/bugs/2842/"><span class="icon"></span>https://sourceforge.net/p/maxima/bugs/2842/</a>
</p>
TicketKarl-Dieter CrismanMon, 17 Nov 2014 16:47:36 GMT
https://trac.sagemath.org/ticket/11164#comment:16
https://trac.sagemath.org/ticket/11164#comment:16
<p>
See also the documentation for <a class="ext-link" href="http://maxima.sourceforge.net/docs/manual/maxima_15.html#Item_003a-gamma_005fexpand"><span class="icon"></span>gamma_expand</a>. The <a class="ext-link" href="http://www.walkingrandomly.com/images/maxima/ChangeLog-5.17-special-functions.txt"><span class="icon"></span>changelog says</a>
</p>
<pre class="wiki">
Implementation of the Incomplete Gamma function
New Maxima User function: gamma_incomplete(a,z)
The following features are implemented:
- Evaluation for real and complex numbers in double float and
bigfloat precision
- Special values for gamma_incomplete(a,0) and gamma_incomplete(a,inf)
- When $gamma_expand T expand the following expressions:
gamma_incomplete(0,z)
gamma_incomplete(n+1/2)
gamma_incomplete(1/2-n)
gamma_incomplete(n,z)
gamma_incomplete(-n,z)
gamma_incomplete(a+n,z)
gamma_incomplete(a-n,z)
- Mirror symmetry
- Derivative wrt the arguments a and z
--------------------------------------------------------------------------------
Implementation of the Generalized Incomplete Gamma function
New Maxima User function: gamma_incomplete_generalized(a,z1,z2)
The following features are implemented:
- Evaluation for real and complex numbers in double float and
bigfloat precision
- Special values for:
gamma_incomplete_generalized(a,z1,0)
gamma_incomplete_generalized(a,0,z2),
gamma_incomplete_generalized(a,z1,inf)
gamma_incomplete_generalized(a,inf,z2)
gamma_incomplete_generalized(a,0,inf)
gamma_incomplete_generalized(a,x,x)
- When $gamma_expand T and n an integer expand
gamma_incomplete_generalized(a+n,z1,z2)
- Implementation of Mirror symmetry
- Derivative wrt the arguments a, z1 and z2
--------------------------------------------------------------------------------
Implementation of the Regularized Incomplete Gamma function
New Maxima User function: gamma_incomplete_regularized(a,z)
The following features are implemented:
- Evaluation for real and complex numbers in double float and
bigfloat precision
- Special values for:
gamma_incomplete_regularized(a,0)
gamma_incomplete_regularized(0,z)
gamma_incomplete_regularized(a,inf)
- When $gamma_expand T and n a positive integer expansions for
gamma_incomplete_regularized(n+1/2,z)
gamma_incomplete_regularized(1/2-n,z)
gamma_incomplete_regularized(n,z)
gamma_incomplete_regularized(a+n,z)
gamma_incomplete_regularized(a-n,z)
- Derivative wrt the arguments a and z
- Implementation of Mirror symmetry
</pre>
TicketKarl-Dieter CrismanMon, 17 Nov 2014 17:01:45 GMT
https://trac.sagemath.org/ticket/11164#comment:17
https://trac.sagemath.org/ticket/11164#comment:17
<blockquote class="citation">
<blockquote class="citation">
<blockquote class="citation">
<p>
What if we don't allow incomplete gamma to become <code>Ei</code> when the first arg is zero - would that solve the branch problem?
</p>
</blockquote>
<p>
Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).
</p>
</blockquote>
<p>
Or maybe not;
</p>
</blockquote>
<p>
Drivel. We could try this, though the branch problem is likely still there with incomplete gamma.
</p>
<pre class="wiki">sage: maxima_calculus('integrate(sin(x)/x,x)')
-(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2
sage: maxima_calculus('integrate(sin(x)/x,x)').sage()
-1/2*I*Ei(I*x) + 1/2*I*Ei(-I*x)
</pre><p>
because we do
</p>
<pre class="wiki"> if x == 0:
return -Ei(-y)
</pre><p>
in sage/functions/other.py
</p>
<p>
But unfortunately
</p>
<pre class="wiki">sage: i*CC(0).gamma_inc(-i/100000)-i*CC(0).gamma_inc(i/100000)
-3.14157265358979
</pre><p>
So same problem there. Yuck. We really need Maxima to provide a proper <code>Si</code> function, or somehow magically translate such combos ourselves with pattern-matching, which sounds horrible too.
</p>
TicketJakob KroekerTue, 25 Aug 2015 15:48:24 GMTstopgaps set
https://trac.sagemath.org/ticket/11164#comment:18
https://trac.sagemath.org/ticket/11164#comment:18
<ul>
<li><strong>stopgaps</strong>
set to <em>todo</em>
</li>
</ul>
TicketJohn ScottSat, 20 Jul 2019 16:17:32 GMT
https://trac.sagemath.org/ticket/11164#comment:19
https://trac.sagemath.org/ticket/11164#comment:19
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/11164#comment:17" title="Comment 17">jakobkroeker</a>:
</p>
<blockquote class="citation">
<p>
We really need Maxima to provide a proper <code>Si</code> function
</p>
</blockquote>
<p>
Maxima has provided a proper Si function for a while called <code>expintegral_si(x)</code> which seems to fit the bill. Maxima doesn't seem to use it though, and the behavior you described still occurs:
</p>
<p>
<code>sage: maxima_calculus('integrate(sin(x)/x,x)')</code>
<code>-(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2</code>
</p>
<p>
Edit: upstream said this would be fixed <a class="ext-link" href="https://sourceforge.net/p/maxima/bugs/3079/#94bd"><span class="icon"></span>here</a> in 2016.
</p>
TicketMarkus WageringelSun, 21 Jul 2019 06:30:10 GMTupstream changed; keywords set
https://trac.sagemath.org/ticket/11164#comment:20
https://trac.sagemath.org/ticket/11164#comment:20
<ul>
<li><strong>keywords</strong>
<em>maxima</em> added
</li>
<li><strong>upstream</strong>
changed from <em>Reported upstream. No feedback yet.</em> to <em>Reported upstream. Developers acknowledge bug.</em>
</li>
</ul>
TicketJohn ScottThu, 01 Aug 2019 13:56:09 GMTowner deleted
https://trac.sagemath.org/ticket/11164#comment:21
https://trac.sagemath.org/ticket/11164#comment:21
<ul>
<li><strong>owner</strong>
<em>Burcin Erocal</em> deleted
</li>
</ul>
TicketFrédéric ChapotonSat, 24 Aug 2019 19:13:02 GMTkeywords changed
https://trac.sagemath.org/ticket/11164#comment:22
https://trac.sagemath.org/ticket/11164#comment:22
<ul>
<li><strong>keywords</strong>
<em>integration</em> added
</li>
</ul>
TicketPeter BruinThu, 20 Aug 2020 06:33:13 GMT
https://trac.sagemath.org/ticket/11164#comment:23
https://trac.sagemath.org/ticket/11164#comment:23
<p>
Similar bug (regression after upgrade to Maxima 5.44.0): <a class="new ticket" href="https://trac.sagemath.org/ticket/30389" title="#30389: defect: Wrong integrate(log(cot(x)-1),x,0,pi/4) (regression) (new)">#30389</a>
</p>
Ticket