[sympy] Fix issue 3503 (#1662)

[raoulb]

Fix Issue 3503: Special values for Si are not implemented.

This fix adds values for Si, Ci, Shi, Chi at oo and -oo.

---

You can merge this Pull Request by running:

  git pull https://github.com/raoulb/sympy fix_issue_3503

Or view, comment on, or merge it at:

  https://github.com/sympy/sympy/pull/1662

Commit Summary

• Implement known values of Si, Ci, Shi, Chi at oo and -oo
• Tests for values of Si, Ci, Shi, Chi at oo and -oo

File Changes

• M sympy/functions/special/error_functions.py (13)
• M sympy/functions/special/tests/test_error_functions.py (10)

Patch Links

• https://github.com/sympy/sympy/pull/1662.patch
• https://github.com/sympy/sympy/pull/1662.diff

—
Reply to this email directly or view it on GitHub.

[Christopher Smith]

In sympy/functions/special/error_functions.py:

> @@ -669,6 +676,8 @@ class Ci(TrigonometricIntegral):
>  
>      _trigfunc = C.cos
>      _atzero = S.ComplexInfinity
> +    _atinf = 0

0->S.Zero

[Christopher Smith]

Other than the small change indicated, this looks fine (assuming the values are correct).

[raoulb]

assuming the values are correct

I checked them with MMA (Wolfram Alpha and Tables) and Fricas (Axiom).

The only issue is limit(Chi(x), x, -oo). The tables at:

http://functions.wolfram.com/GammaBetaErf/CoshIntegral/03/02/

say it's oo but if we evaluate numerically (with both, MMA and mpmath)
we get an imaginary part if pi. (In Fricas this value is undefined.)

[Aaron Meurer]

This is obviously unrelated to this PR, but at some point, we should refactor these class-level definitions into @property methods. The reason is that they are run at import time, and the creation of objects like 2*pi takes more time then compiling a function (and it also leads to certain bugs making SymPy not import, like issue 3235).

But good work here. If you're confident that these values are correct, this can be merged.

[raoulb]

we should refactor these class-level definitions into @property methods.

I never used the proptery. Should it be something like:

    @property
    def _atzero(self):
        return S(0)

But this does not seem to work.

If you're confident that these values are correct, this can be merged.

The only one I'd like to get a second opinion is Chi(-oo). Here is the comparison to mpmath:

In [17]: mpmath.chi(-10000)
Out[17]: mpc(real='4.4038495418373575e+4338', imag='3.1415926535897931')

In [18]: Chi(-oo)
Out[18]: oo

But good work here.

Thanks.

[Aaron Meurer]

Is Chi defined by an integral for negative x? If so, what does meijerint give?

[raoulb]

> Is Chi defined by an integral for negative x?

The general definition is:

http://upload.wikimedia.org/math/9/0/3/9036287a03d2334c8db2af61b1b83300.png
http://functions.wolfram.com/GammaBetaErf/CoshIntegral/02/

The issues may arise from the log(x) in there.

Sympy, Maxima and Fricas as well as MMA/WA all give:

In [12]: limit(log(x), x, -oo)
Out[12]: oo

http://www.wolframalpha.com/input/?i=limit%28log%28x%29%2C%20x%2C%20-oo%29&dataset=

In contrast numerical evaluation shows:

In [3]: mpmath.log(-100000)
Out[3]: mpc(real='11.512925464970229', imag='3.1415926535897931')

and in Maxima

log(-1000), bfloat;
3.141592653589793b0*%i+6.907755278982137b0

Therefore we have the very same problem with the common log(x)
function. Any solution we agree upon there should be applied
consistently to Chi(x).

To me it seems that numerics add a +I*pi and symbolics
keep the limit real.

> If so, what does> meijerint give?

As strange as it may sound: the only useful solution
is what we get without explicitly invoking the MeijerG code!

In [13]: integrate((cosh(t)-1)/t, (t, 0, x), meijerg=True)
Out[13]: Integral((cosh(t) - 1)/t, (t, 0, x))

In [14]: integrate((cosh(t)-1)/t, (t, 0, oo), meijerg=True)
Out[14]: Integral((cosh(t) - 1)/t, (t, 0, oo))

In [15]: integrate((cosh(t)-1)/t, (t, 0, oo))
Out[15]: Integral((cosh(t) - 1)/t, (t, 0, oo))

In [16]: integrate((cosh(t)-1)/t, (t, 0, x))
Out[16]: -2*log(x) + log(x**2)/2 + Chi(x) - EulerGamma

[Aaron Meurer]

According to the docstring, that's the definition for positive z. For negative z, we have to use analytic continuation, or the second formula at http://docs.sympy.org/0.7.2/modules/functions/special.html?highlight=chi#sympy.functions.special.error_functions.Chi. I guess there's a branch cut on the negative axis, so Chi(-oo) is not well-defined.

@ness01 any thoughts?

[raoulb]

According to the docstring, that's the definition for positive z. For negative z, we have to use analytic continuation

OTOH Wolfram functions (who seem to be very concerned about branch cuts and such things) do
not give any restrictions on z.

And there is also this one:

http://functions.wolfram.com/GammaBetaErf/CoshIntegral/introductions/ExpIntegrals/04/04/imagetext/0087/text87.gif

yielding

In [9]: mpmath.ci(-1000*1.0j) + mpmath.log(-1000) - mpmath.log(-1000*1.0j)
Out[9]: mpc(real='9.8602256857061915e+430', imag='3.1415926535897931')

I guess there's a branch cut on the negative axis, so Chi(-oo) is not well-defined.

Yes. This is probably the source of all confusion.

[Aaron Meurer]

Not using MeijerG is not helpful, because we aren't as assured of correctness in these corner cases then.

By the way, the graph of http://www.wolframalpha.com/input/?i=Chi%28z%29 seems to agree with -oo + pi.

[raoulb]

> Not using MeijerG is not helpful, because we aren't as assured of
> correctness in these corner cases then.

All others are unambiguous manual table lookups.

> By the way, the graph of
> http://www.wolframalpha.com/input/?i=Chi%28z%29 seems to agree with
> -oo + pi.

But the point is that we have the same behavior for log:

http://www.wolframalpha.com/input/?i=log%28z%29

and return oo as limit for log(-oo).

[Aaron Meurer]

I don't see how log(-oo) can be oo, given:

In [15]: log(-1e10000 - 0.00000001*I).evalf()
Out[15]: -+inf - 3.14159265358979⋅ⅈ

In [16]: log(-1e10000 + 0.00000001*I).evalf()
Out[16]: -+inf + 3.14159265358979⋅ⅈ

So below the branch, it is oo + pi*I and above it it is oo - pi*I. Actually, it appears to actually be computing this as -oo, not oo (the -+ is a printing bug I guess).

In [17]: log(-1e10000 + 0.00000001*I).evalf().args
Out[17]: (-inf, 3.14159265358979⋅ⅈ)

In [19]: log(-1e10000 + 0.00000001*I).evalf().args[0] == oo
Out[19]: False

In [20]: log(-1e10000 + 0.00000001*I).evalf().args[0] == -oo
Out[20]: True

[raoulb]

> I don't see how log(-oo) can be oo, given:
>

> ```py

> In [15]: log(-1e10000 - 0.00000001*I).evalf()
> Out[15]: -+inf - 3.14159265358979⋅ⅈ
>
> In [16]: log(-1e10000 + 0.00000001*I).evalf()
> Out[16]: -+inf + 3.14159265358979⋅ⅈ
>
> ```
>
> So below the branch, it is oo + pi*I and above it it is oo - pi*I.
> Actually, it appears to actually be computing this as -oo, not oo
> (the -+ is a printing bug I guess).
>

> ```py

> In [17]: log(-1e10000 + 0.00000001*I).evalf().args
> Out[17]: (-inf, 3.14159265358979⋅ⅈ)
>
> In [19]: log(-1e10000 + 0.00000001*I).evalf().args[0] == oo
> Out[19]: False
>
> In [20]: log(-1e10000 + 0.00000001*I).evalf().args[0] == -oo
> Out[20]: True

> ```

I was referring to the symbolic values:

```
In [5]: gruntz(log(x), x, -oo)
Out[5]: oo

In [6]: log(-oo)
Out[6]: oo
```

This is what sympy returns. Independently of the question if this is
correct, I think Chi(x) should have the same behavior, either both
add the I*pi or none.

[Julien Rioux]

SymPy Bot Summary: Passed after merging raoulb/fix_issue_3503 (a961b1c) into master (211f0f7).
Python 2.7.2-final-0: pass
Python 3.2.1-final-0: pass
Sphinx 1.1.3: pass

[Aaron Meurer]

The symbolic values and numerical values should agree.

Maybe the real answer is zoo, as at 0 (the other branch point of log).

[raoulb]

> The symbolic values and numerical values should agree.

They should. But compare to:

http://www.wolframalpha.com/input/?i=log%28-oo%29&dataset=

http://www.wolframalpha.com/input/?i=N[log%28-100000000%29]


> Maybe the *real* answer is zoo, as at 0 (the other branch point of
> log).

At least these tables say otherwise:

http://functions.wolfram.com/ElementaryFunctions/Log/03/ShowAll.html

There is no complex infinity in the values.
And I'm not sure if ``oo+i*pi`` is really equivalent
to a complex infinity. Because ``zoo`` is a complex "number"
with infinite magnitude and *arbitrary* angle.

[Aaron Meurer]

I rather see zoo as the identification of all complex numbers of infinite magnitude. Like zero, the magnitude of zoo is meaningless.

And by the way, wouldn't arg(oo + pi*I) == arg(oo)?

[Aaron Meurer]

Anyway, lets go ahead and agree with Wolfram on the symbolic values. But I would still like to understand better why log(-oo) is evaluated this way.

[Aaron Meurer]

In other words, I'm +1 to this branch, but I'm still curious. Maybe someone has written a paper on this?

[Julien Rioux]

I would remove the symbolic values Ci(-oo) and Chi(-oo) because they lie on the (usual choice of) branch cut for ln. ln(-oo) seems wrong too, by the way.

And by the way, wouldn't arg(oo + pi*I) == arg(oo)?

It depends where you put your branch cut.

[raoulb]

What is the state of this PR now? To what can we agree in order to merge this?

[Aaron Meurer]

I stand by what I said above. I'm running the sympy-bot tests now. When they come in, let's merge this, unless they show errors, or if someone feels that we should go a different way and discuss this further.

[raoulb]

Then let's merge, no?

[Aaron Meurer]

SymPy Bot Summary: Failed after merging raoulb/fix_issue_3503 (a961b1c) into master (2abef4f).
@raoulb: Please fix the test failures.
Python 2.5.0-final-0: pass
Python 2.6.6-final-0: pass

Python 2.7.2-final-0: pass

Python 2.6.8-final-0: pass
Python 2.7.3-final-0: fail
PyPy 2.0.0-beta-1; 2.7.3-final-42: pass
Python 3.2.2-final-0: pass
Python 3.3.0-final-0: pass
Python 3.2.3-final-0: pass
Python 3.3.0-final-0: pass
Python 3.3.0-final-0: pass
**Sphinx 1.1.3:** pass

[raoulb]

That are plotting failures. I think this is not related to my changes, is it?

[Aaron Meurer]

Yeah let's merge.

[Aaron Meurer]

Merged #1662.