Integration result appears to use a different branch of polylog

[normalhuman]

The following is close to being correct:

>>> x = symbols('x', positive=True)
>>> integrate(log(x)/(1+x), (x, 0, 1))
pi**2/6 - polylog(2, 2) + I*pi*log(2)    # -0.822467033424113 + 4.3551721806072*I

The real part is fine. The imaginary part is (obviously) wrong. Indeed, polylog(2, 2) is pi**2/4 - I*pi*log(2) and so we end up with extra 2*I*pi*log(2).

My first guess was that polylog terms comes from somewhere where it's understood to have a different branch cut, equating to pi**2/4 + I*pi*log(2) at 2.

Or perhaps it's the branching of log that's treated wrong here. The indefinite integral is returned as

>>> integrate(log(x)/(1+x), x)
I*pi*log(x + 1) - polylog(2, x + 1)

This is incorrect; the derivative of that function is (log(-x) + i*pi) / (x+1), which isn't equivalent to log(x)/(x+1), even for positive x.

Why SymPy makes this mistake seems related to #13853. From its point of view, the derivative is -polylog(1, x + 1)/(x + 1) + I*pi/(x + 1) (correct) which it turns (with current form of expand_func) into

log((x + 1)*exp_polar(-I*pi) + 1)/(x + 1) + I*pi/(x + 1)

Using that log(exp_polar(-I*pi)) = - I *pi, we get log(x + exp_polar(I*pi) + 1)/(x + 1) which is log(x)/(x+1). So the antiderivative seems correct, even though it isn't. Faulty expand_func makes a wrong antiderivative appear to be correct.

[jksuom]

I suspect that the integration result is oversimplified at some stage. meijerint_indefinite gives a more complete answer (though that is not necessarily correct either).

>>> meijerint_indefinite(log(x)/(1 + x), x)
Piecewise((I*pi*log(x + 1) - polylog(2, x + 1), Abs(x + 1) < 1), (-I*pi*log(1/(x + 1)) - polylog(2, x + 1), Abs(1/(x + 1)) < 1), (-I*pi*meijerg(((), (1, 1)), ((0, 0), ()), x + 1) + I*pi*meijerg(((1, 1), ()), ((), (0, 0)), x + 1) - polylog(2, x + 1), True))

[normalhuman]

It's still the same thing, unfrortunately: the relevant piece is -I*pi*log(1/(x + 1)) - polylog(2, x + 1) which for positive x is equivalent to I*pi*log(x + 1) - polylog(2, x + 1). Would be correct if the logarithm term had opposite sign.

[normalhuman]

I think I found the source of this wrong result:

>>> var('z')
>>> hyperexpand(meijerg(((1, 1, 1), ()), ((1,), (0, 0)), z))
-polylog(2, z*exp_polar(I*pi))

By Wolfram Functions this is -polylog(2, -z) which seems the same but I think the specific exp_polar used here is not always picking the right branch. The choice between exp_polar(I*pi) and exp_polar(-I*pi) may have to depend on |z| being greater or less than 1. Or perhaps -polylog(2, - (z*exp_polar(I*pi))) could be given as the correct answer in all cases, meaning the evaluation of polylog(2, -w) at the point z*exp_polar(I*pi).

---

Why I think so: I traced the computation step by step. First, meijerint_indefinite shifts the function to log(x-1)/x. The power of x is put aside and log(x-1) is expanded as

>>> _rewrite_single(log(x-1), x)
([(I*pi, 0, meijerg(((1,), ()), ((), (0,)), x)), 
(I*pi, 0, meijerg(((), (1,)), ((0,), ()), x)), 
(1, 0, meijerg(((1, 1), ()), ((1,), (0,)), x*exp_polar(-I*pi)))], True)

This reads as I*pi*H(|x| - 1) + I*pi*H(1 - |x|) + log(1 + x*exp_polar(-I*pi)) where H is the Heaviside function. Seems correct so far: in order to use Meijer G representation of logarithm, one needs log(1+z), so log(-1+x) is turned into I*pi + log(1 + x*exp_polar(-I*pi)).

This is going to be tricky: since x > 1 here, the expression 1 + x*exp(-I*pi) is at the negative half-axis, but we don't use the principal branch of log there: we need one that has imaginary part -I*pi. The use of polar number exp_polar(-I*pi) is crucial: it can't be just -1. Unfortunately we are now also in the business of adding 1 to a polar number.

Next, these are integrated.

1. meijerg(((1,), ()), ((), (0,)), x) times 1/x integrates to meijerg(((1, 1), ()), ((), (0, 0)), x)
2. meijerg(((), (1,)), ((0,), ()), x) times 1/x integrates to -meijerg(((), (1, 1)), ((0, 0), ()), x)
3. meijerg(((1, 1), ()), ((1,), (0,)), x*exp_polar(-I*pi)) times 1/x integrates to meijerg(((1, 1, 1), ()), ((1,), (0, 0)), x*exp_polar(-I*pi))

Here, meijerg(((1, 1), ()), ((), (0, 0)), x) is -log(1/x) * H(|x| - 1) and -meijerg(((), (1, 1)), ((0, 0), ()), x) is log(x) * H(1 - |x|). I agree with both of those, considering that x is positive. After shifting the variable back, these produce the term I*pi*log(x+1) that we see in the answer. I see no errors in this part.

So the error has to be in hyperexpand reducing meijerg(((1, 1, 1), ()), ((1,), (0, 0)), x*exp_polar(-I*pi)) to -polylog(2, x). I think this should have been -polylog(2, x*exp_polar(-2*I*pi)), which is another branch because x > 1 (we go around the singularity at 1). It seems hyperexpand is oversimplifying the picture: its answer would be valid for 0 < x < 1.

[normalhuman]

For a concrete example, take x = 3. The value of meijerg(((1, 1, 1), ()), ((1,), (0, 0)), 3*exp_polar(-I*pi)) can be found numerically by taking the limit as

>>> meijerg(((1, 1, 1), ()), ((1,), (0, 0)), 3*exp(-0.99999*I*pi)).evalf()
-2.32008172800931 - 3.4513705193623*I

On the other hand, hyperexpand turns this G function into -polylog(2, x) which for x=3 is -2.3201804233131 + 3.4513922952232*I. Wrong side of branch cut.

I suspect the issue is in do_slater method of hyperexpand which rewrites the G function in terms of hyper((1, 1, 1), (2, 2)). In the rewrite, something seems to be lost from the polarity of the argument.

[jksuom]

_rewrite_single(log(x-1), x) seems to work by writing

log(x - 1) = log((-1)*(1 - x)) = log(-1) + log(1 - x)

Now, the question is, which value should be taken for log(-1). I think that the value should be -I*pi to cancel the extra I*pi in the principal branch of log(1 - x) for x > 1.

On the whole, this computation by means of Meijer G-functions is on uncertain basis since it depends on the analytic continuation of the functions to a boundary cut. For example, the expression of log(1 + x) as meijerg(((1, 1), ()), ((1,), (0,)), x) can only be applied in the domain |x| < 1 where the integral is convergent.

[normalhuman]

Upon further reflection, it seems the implementation of Slater's theorem is correct but perhaps its result is being over-simplified. Since we have a G-function with p=q=3, one has to deal with |z|<1 and |z|>1 separately, which is correctly implemented in hyperexpand which calls do_slater twice.

When reducing meijerg(((1, 1, 1), ()), ((1,), (0, 0)), x*exp_polar(-I*pi)), hyperexpand gets two formulas from do_slater:

1. -polylog(2, x) provided that Abs(x) < 1
2. log(exp_polar(I*pi)/x)**2/2 + polylog(2, exp_polar(2*I*pi)/x) + pi**2/6 provided that Abs(x)>1

But then it considers whether the region of convergence of G-function is connected and if so, decides that the first formula is universally valid by the uniqueness of analytic continuation.

This may be an oversimplification because the second formula had information regarding branch choice that the first one did not. Indeed, evaluating the first formula at x = 3:

-polylog(2, 3).evalf() is -2.3201804233131 + 3.4513922952232*I which leads to incorrect answer; wrong side of branch cut.

Evaluating the second formula at x = 3:

>>> (log(exp_polar(I*pi)/x)**2/2 + polylog(2, exp_polar(2*I*pi)/x) + pi**2/6).subs(x, 3).evalf()
polylog(2, exp_polar(2*I*pi)/3) - 2.68639365329016 - 3.4513922952232*I

Polylog is unbranched in the unit disk, so the first term is same as polylog(2, 1/3). Putting this in we get -2.3201804233131 - 3.4513922952232*I - now that is the correct side of branch cut.

As long as polylog and similar functions are concerned, one should probably prefer piecewise formulas that keep the argument within the unit disk to the simpler formulas that end up evaluating on the branch cut.

[jksuom]

After spending quite a lot of time on various G-functions I noticed that the integral can be computed quite easily by partial integration using log(1+x) = -polylog(1, -x) and diff(polylog(2, x), x) = polylog(1, x)/x to get

Integral(log(t)/(1+t), (t, 0, x)) = log(x)*log(1+x) + polylog(2, -x)

This is valid in the domain |arg(x)| < pi where the principal branches of log(x), log(1+x) and polylog(2, -x) are defined. SymPy can also verify this

>>> f = log(x)*log(1+x) + polylog(2, -x)
>>> expand_func(diff(f, x))
log(x)/(x + 1)

This result can be connected to that given by SymPy by means of the reflection formula of dilogarithm:

polylog(2, x) + polylog(2, 1 - x) = pi**2/6 - log(x)*log(1 - x)

which is valid in the complement of the real axis (and, in addition, for 0 < x < 1 where neither x nor 1 - x is real and negative.) It follows that the value of Integral(log(t)/(1+t), (t, 0, x)) for non-real x can be written

pi**2/6  + s*pi*I*log(1 + x) - polylog(2, 1 + x),

where s = sign(im(x)). So the formula given by SymPy is correct in the upper half plane but wrong in the lower half plane. As the value of polylog(2, 1 + x) on the branch cut 1 + x > 1 is the limit from the lower half plane by convention, the sign of the pi*I*log(1 + x) term is wrong for x real and positive.