Fix issue 3503

[raoulb]

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

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

[smichr]

0->S.Zero

[smichr]

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

[asmeurer]

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

[asmeurer]

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, 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 evaluations (Maxima, Fricas, mpmath, MMA/WA) show:

In [6]: mpmath.log(-1000)
Out[6]: mpc(real='6.9077552789821368', imag='3.1415926535897931')

and in Maxima

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

and Fricas

(6) -> log(-1000::Complex(Float))
(6) 6.9077552789_821370521 + 3.1415926535_897932385 %i

Therefore we have the very same problem with the more 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

[asmeurer]

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.

[asmeurer]

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

[asmeurer]

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:

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

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.

[jrioux]

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

[asmeurer]

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.

[asmeurer]

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

[asmeurer]

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.

[asmeurer]

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

[jrioux]

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?

[asmeurer]

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?

[asmeurer]

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?

[asmeurer]

Yeah let's merge.