Opened 9 years ago

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

Reported by: Owned by: benreynwar major calculus maxima, integration Reported upstream. Developers acknowledge bug. todo

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

### comment:1 Changed 9 years ago by kcrisman

• Component changed from PLEASE CHANGE to calculus
• Owner changed from tbd to burcin
• Priority changed from minor to major
• Report Upstream changed from N/A to Fixed upstream, in a later stable release.

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.

### comment:2 Changed 8 years ago by dsm

• Status changed from new to needs_review

Might as well doctest and close.

### comment:3 Changed 8 years ago by kcrisman

• Status changed from needs_review to needs_work

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!

### comment:4 Changed 8 years ago by dsm

• Status changed from needs_work to needs_info

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?

### comment:5 follow-up: ↓ 6 Changed 8 years ago by kcrisman

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.

### comment:6 in reply to: ↑ 5 Changed 8 years ago by kcrisman

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

### comment:7 Changed 8 years ago by kcrisman

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

### comment:8 follow-up: ↓ 9 Changed 7 years ago by kcrisman

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

### comment:9 in reply to: ↑ 8 Changed 6 years ago by kcrisman

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.

### comment:10 Changed 5 years ago by tscrim

• Description modified (diff)

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

### comment:11 follow-up: ↓ 12 Changed 5 years ago by pbruin

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.

### comment:12 in reply to: ↑ 11 ; follow-up: ↓ 13 Changed 5 years ago by kcrisman

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?

### comment:13 in reply to: ↑ 12 ; follow-up: ↓ 14 Changed 5 years ago by kcrisman

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?

### comment:14 in reply to: ↑ 13 ; follow-up: ↓ 15 Changed 5 years ago by kcrisman

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

Last edited 5 years ago by kcrisman (previous) (diff)

### comment:15 in reply to: ↑ 14 ; follow-up: ↓ 17 Changed 5 years ago by kcrisman

• Report Upstream changed from Fixed upstream, in a later stable release. to Reported upstream. No feedback yet.

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/

Last edited 5 years ago by kcrisman (previous) (diff)

### comment:16 Changed 5 years ago by kcrisman

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

```

### comment:17 in reply to: ↑ 15 Changed 5 years ago by kcrisman

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.

### comment:18 follow-up: ↓ 19 Changed 4 years ago by jakobkroeker

• Stopgaps set to todo

### comment:19 in reply to: ↑ 18 Changed 5 months ago by gh-JohnScott623

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.

Last edited 5 months ago by gh-JohnScott623 (previous) (diff)