Integral of sin(x)/x gives false result.

[sagetrac-benreynwar]

> eq = sin(x)/x
> integrate(eq, x, -1e-6, 1e-6)
-3.14159065359

I expected an answer of about 2e-6.

Upstream: Reported upstream. Developers acknowledge bug.

Component: calculus

Keywords: maxima, integration

Stopgaps: todo

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

[kcrisman]

Upstream: Fixed upstream, in a later stable release.

[kcrisman]

comment:1

In the future, be sure to pick a component (for instance, calculus or symbolics); that will help people find it more easily.

There are two problems here.

• The actual problem - the report unfortunately has a typo.

sage: integrate(eq,x,-1e-6,1e-6)
-3.14159065359

This is definitely a problem. However, it's fixed in the latest Maxima. Here is the old one, in Sage:

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

And here is 5.24.0. I assume this isn't a problem of which Lisp we are using?

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

• A problem revealed, but also fixed in the latest. Old:

(%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.

New:

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

We'd need doctesting to test both of these once we upgrade Maxima.

[sagetrac-dsm]

add doctest to confirm fix

[sagetrac-dsm]

comment:2

Attachment: trac_11164_verify_integration_sinx_x.patch.gz

Might as well doctest and close.

[kcrisman]

comment:3

I may have misstated whether this is "fixed".

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)

The graph makes it very obvious the second (numerical) answer is correct. Unsurprisingly, the other answer is off by exactly pi...

Presumably the first answer is "correct" in the sense that it is an alternate representation of the Si sine integral function. But I think there are some branch issues - Ei(I*x) only has an alternate formula for x>0.

So we still need to do something with this. Or at the very least warn about branches. But really, sin(x)/x shouldn't have to worry about this!

[sagetrac-dsm]

comment:4

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. :^)

Does this work on the current maxima trunk now?

[kcrisman]

comment:5

This is the previous version (5.25.1):

(%i5) integrate(sin(x)/x,x,-1/1000000,1/1000000);

(%o5) %i*gamma_incomplete(0,-%i/1000000)-%i*gamma_incomplete(0,%i/1000000)

W|A agrees that this is about -pi. (To be precise, 2*Si(1/10^6)-pi.)

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

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 (expintegral_si), it isn't really connected with the rest of Maxima (?).

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 Si there is robust enough to handle all this yet.

[kcrisman]

comment:6

Here is the issue; this function may or may not be the same as the sine integral.

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 x is, that's annoying. I think the Maxima folks would say this is a feature (for symbolic reasons).

[kcrisman]

comment:7

By the way, achrzesz has a great one-liner for a related question:

 maxima('integrate(diff((exp(x)-1)/x,x),x),gamma_expand:true,factor').sage()

Though here it doesn't solve the problem, but maybe it will trigger someone's memory as to another Maxima flag to set...

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

[kcrisman]

comment:8

I wonder whether this is potentially fixed in #13364 - I don't know if we ever reported upstream, or whether we should have...

[kcrisman]

comment:9

I wonder whether this is potentially fixed in #13364 - I don't know if we ever reported upstream, or whether we should have...

Apparently not yet.

[tscrim]

comment:10

In 6.3.beta4 (so with #13973), I get:

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

Something bad happens with the integration computation when crossing over 0 as we have:

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

[pjbruin]

comment:11

The problem is caused by the integral being computed via symbolic integration; see also comment:5. Maxima integrates exp(x)/x to -gamma_incomplete(0, -x), 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 integrate(sin(x)/x, x, -1, 1), the two terms coming from exp(%i*x) and exp(-%i*x) 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 sin(x)/x, which is an entire function.

[kcrisman]

comment:12

The problem is caused by the integral being computed via symbolic integration; see also comment:5. Maxima integrates exp(x)/x to -gamma_incomplete(0, -x), 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 integrate(sin(x)/x, x, -1, 1), the two terms coming from exp(%i*x) and exp(-%i*x) 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 sin(x)/x, which is an entire function.

How much you wanna bet this is the same underlying problem as in #17328?

[kcrisman]

comment:13

How much you wanna bet this is the same underlying problem as in #17328?

Nah, I was thrown off by the gammas, I think.

What if we don't allow incomplete gamma to become Ei when the first arg is zero - would that solve the branch problem?

[kcrisman]

comment:14

What if we don't allow incomplete gamma to become Ei when the first arg is zero - would that solve the branch problem?

Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).

[kcrisman]

Changed upstream from Fixed upstream, in a later stable release. to Reported upstream. No feedback yet.

[kcrisman]

comment:15

What if we don't allow incomplete gamma to become Ei when the first arg is zero - would that solve the branch problem?

Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).

Or maybe not;

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

so maybe we are already setting gamma_expand:true somewhere nowadays. Reported upstream at https://sourceforge.net/p/maxima/bugs/2842/

[kcrisman]

comment:16

See also the documentation for gamma_expand. The changelog says

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

[kcrisman]

comment:17

What if we don't allow incomplete gamma to become Ei when the first arg is zero - would that solve the branch problem?

Looks like if that was once possible, it isn't any more (already done during integration in Maxima, I mean).

Or maybe not;

Drivel. We could try this, though the branch problem is likely still there with incomplete gamma.

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)

because we do

        if x == 0:
            return -Ei(-y)

in sage/functions/other.py

But unfortunately

sage: i*CC(0).gamma_inc(-i/100000)-i*CC(0).gamma_inc(i/100000)
-3.14157265358979

So same problem there. Yuck. We really need Maxima to provide a proper Si function, or somehow magically translate such combos ourselves with pattern-matching, which sounds horrible too.

[sagetrac-jakobkroeker]

Stopgaps: todo

[JohnScott623]

comment:19

Replying to @sagetrac-jakobkroeker:

We really need Maxima to provide a proper Si function

Maxima has provided a proper Si function for a while called expintegral_si(x) which seems to fit the bill. Maxima doesn't seem to use it though, and the behavior you described still occurs:

sage: maxima_calculus('integrate(sin(x)/x,x)')
-(%i*gamma_incomplete(0,%i*x)-%i*gamma_incomplete(0,-%i*x))/2

Edit: upstream said this would be fixed here in 2016.

[mwageringel]

Changed keywords from none to maxima

[mwageringel]

Changed upstream from Reported upstream. No feedback yet. to Reported upstream. Developers acknowledge bug.

[fchapoton]

Changed keywords from maxima to maxima, integration

[pjbruin]

comment:23

Similar bug (regression after upgrade to Maxima 5.44.0): #30389