hyperexpand bad output called on I*meijerg(((1, x + 2), ()), ((x + S(1)/2,), (0,)), z*exp_polar(I*pi))*exp(-I*pi*x)/(2*sqrt(pi))

[ehren]

Here's an example where hyperexpand changes the numerical evaluation of its input. The origin of the expression in the title is obscure (see below). I'm also not sure if the bad output of hyperexpand is just a bug or if it represents an undocumented could-not-expand error state (I'd guess it's a bug).

from sympy import *
x = Symbol("x", real=True)
z = Symbol("z", real=True)
res = I*meijerg(((1, x + 2), ()), ((x + S(1)/2,), (0,)), z*exp_polar(I*pi))*exp(-I*pi*x)/(2*sqrt(pi))
hyp = hyperexpand(res)  # I*(-2*sqrt(pi)*z**(x + 1)*exp(-I*pi*(-x - 1))*gamma(x + 1)*hyper((-1/2, -x - 1), (-x,), 1/z)/gamma(x + 2) + gamma(-x - 1)*gamma(x + 1/2))*exp(-I*pi*x)/(2*sqrt(pi))
print(res.subs(x,0).subs(z,0).evalf())  # 0
print(hyp.subs(x,0).subs(z,0).evalf())  # zoo
print(res.subs(x,0).subs(z,1))  # I*meijerg(((1, 2), ()), ((1/2,), (0,)), exp_polar(I*pi))/(2*sqrt(pi))
print(res.subs(x,0).subs(z,1).evalf())  # 1.57079632679490
print(hyp.subs(x,0).subs(z,1))  # I*(2*sqrt(pi)*hyper((-1, -1/2), (0,), 1) + zoo)/(2*sqrt(pi))
print(hyp.subs(x,0).subs(z,1).evalf())  # nan  
print(res.subs(x,1).subs(z,1).evalf())  # 0.392699081698724
print(hyp.subs(x,1).subs(z,1).evalf())  # nan

This is related to #26477 but is a separate issue. Fixing the above bug won't affect #26477 for various reasons.

Aside: Steps to obtain the expression in issue title:

1. use an older copy of sympy. say about 2 years old.
2. use the modified meijerint_indefinite (that returns 2 results) provided by @haru-44 in Error in integral result using hyper #26477
3. disable the call to hyperexpand here:

   sympy/sympy/integrals/meijerint.py

   Line 1736 in a3659be

   r = hyperexpand(r.subs(t, a*x**b), place=place)

   
   that is, replace

r = hyperexpand(r.subs(t, a*x**b), place=place)

with

r = r.subs(t, a*x**b)

4. repeat the call to meijerint_indefinite from issue Error in integral result using hyper #26477. that is call meijerint_indefinite(z**(x - Rational(1, 2))*sqrt(1 - z), z) . One of the 2 results is in the issue title (the other result isn't mangled by hyperexpand).

edit: The "could-not-expand" comments above are based on meijerint checking for and discarding symbolic nan and zoo in various places - but that seems to be a red herring wrt this issue.

[ehren]

Here's a simpler example that avoids exp_polar. hyperexpand introduces a sign change:

>>> m = I*meijerg(((1, x + 2), ()), ((x + S(1)/2,), (0,)), -1)/(2*sqrt(pi))
>>> m
I*meijerg(((1, x + 2), ()), ((x + 1/2,), (0,)), -1)/(2*sqrt(pi))
>>> m.subs(x,0)
I*meijerg(((1, 2), ()), ((1/2,), (0,)), -1)/(2*sqrt(pi))
>>> _.evalf()
1.57079632679490
>>> hyperexpand(m)
I*sqrt(pi)*exp(I*pi*(1/2 - x))*gamma(x + 1/2)/(2*gamma(x + 2))
>>> _.subs(x,0)
-pi/2
>>> m.subs(x,0)
I*meijerg(((1, 2), ()), ((1/2,), (0,)), -1)/(2*sqrt(pi))
>>> hyperexpand(_)
-pi/2

[ehren]

Wolfram agrees with pi/2

[oscarbenjamin]

Distilling this down I got to

In [52]: m = meijerg(((1,), ()), ((1,), ()), 0)

In [53]: m
Out[53]:
╭─╮1, 1 ⎛1    │  ⎞
│╶┐     ⎜     │ 0⎟
╰─╯1, 1 ⎝1    │  ⎠

In [54]: hyperexpand(m)
Out[54]: 1

In [55]: m.n()
Out[55]: 0

[ehren]

two more with the wrong sign after hyperexpand:

meijerg(((0,), ()), ((1/2,), (1,)), I)

and

meijerg(((0,), ()), ((1/2,), (1,)), -1)

[ehren]

It seems the issue is the z param of meijerg. There's already quite a bit of "polarifying" in meijerint and hyperexpand but in the case of e.g. do_slater it seems we need to work with the polar lifted version throughout (too early "unpolarify" in do_slater brings back the old behaviour). I can illustrate with a pull (with polariify applied to the last arg of meijerg in hyperexpand) but this script should suffice for now (and show that this approach is not quite correct ending up with the wrong sign vs wolfram in one case at least).

from sympy import *
x, z = symbols("x z")

test1_problematic = meijerg(((1, x + 2), ()), ((x + S(1)/2,), (0,)), -I*z)
test1a_problematic = meijerg(((1, 2), ()), ((S(1)/2,), (0,)), -I)
test1b = meijerg(((1, 2), ()), ((S(1)/2,), (0,)), -1)

test2 = meijerg(((0,), ()), ((S(1)/2,), (1,)), I)
test2a = meijerg(((0,), ()), ((S(1)/2,), (1,)), -1)

test3 = meijerg(((1,), ()), ((1,), ()), 0)


def polarified_meijerg_z(g: meijerg):
    assert isinstance(g, meijerg)

    z = g.args[2]
    factors = z.as_ordered_factors()
    zp = 1
    zz = 1
    for f in factors:
        # lowercase is_number (rather than is_Number) so that I.is_number is True (just a bad heuristic)
        if f.is_number:
            zp *= polarify(f, lift=True)
        else:
            zz *= f

    return meijerg(g.args[0], g.args[1], zp*zz)


def test_hyperexpand(gs: str, polarify_z=False):
    print(gs)
    g = eval(gs)

    if polarify_z:
        g = polarified_meijerg_z(g)
    print("g:", g)

    h = hyperexpand(g)
    print("h:", h)

    if g.free_symbols:
        print("g.subs(x, 0).subs(z, 1).evalf():", g.subs(x, 0).subs(z, 1).evalf())
        print("h.subs(x, 0).subs(z, 1).evalf():", h.subs(x, 0).subs(z, 1).evalf())
    else:
        print("g.evalf():", g.evalf())
        print("h.evalf():", h.evalf())
    print()


tests = ["test1_problematic", "test1a_problematic", "test1b", "test2", "test2a", "test3"]

print("WITHOUT POLARIFY")

for test in tests:
    test_hyperexpand(test)

print("WITH POLARIFY")

for test in tests:
    test_hyperexpand(test, polarify_z=True)

output:

WITHOUT POLARIFY
test1_problematic
g: meijerg(((1, x + 2), ()), ((x + 1/2,), (0,)), -I*z)
h: -2*sqrt(pi)*z**(x + 1/2)*exp(3*I*pi*(x + 1/2)/2)*gamma(x + 1/2)*hyper((-1/2, x + 1/2), (x + 3/2,), z*exp_polar(5*I*pi/2))/gamma(x + 3/2)
g.subs(x, 0).subs(z, 1).evalf(): -4.32257296433943 + 5.92192068901469*I
h.subs(x, 0).subs(z, 1).evalf(): 4.32257296433943 - 5.92192068901469*I

test1a_problematic
g: meijerg(((1, 2), ()), ((1/2,), (0,)), -I)
h: -4*sqrt(pi)*((1/2 - I/2)/sqrt(1 - exp_polar(I*pi/2)) + exp_polar(-I*pi/4)*asin(exp_polar(I*pi/4))/2)*exp(-I*pi/4)
g.evalf(): -4.32257296433943 + 5.92192068901469*I
h.evalf(): -4.32257296433943 + 5.92192068901469*I

test1b
g: meijerg(((1, 2), ()), ((1/2,), (0,)), -1)
h: I*pi**(3/2)
g.evalf(): -5.56832799683171*I
h.evalf(): 5.56832799683171*I

test2
g: meijerg(((0,), ()), ((1/2,), (1,)), I)
h: (1/2 - I)*exp(-I)*exp(-3*I*pi/4)
g.evalf(): 0.275572216790442 - 1.08354047147912*I
h.evalf(): -0.275572216790442 + 1.08354047147912*I

test2a
g: meijerg(((0,), ()), ((1/2,), (1,)), -1)
h: -3*E*I/2
g.evalf(): 4.07742274268857*I
h.evalf(): -4.07742274268857*I

test3
g: meijerg(((1,), ()), ((1,), ()), 0)
h: 1
g.evalf(): 0
h.evalf(): 1.00000000000000

WITH POLARIFY
test1_problematic
g: meijerg(((1, x + 2), ()), ((x + 1/2,), (0,)), z*exp_polar(3*I*pi/2))
h: -2*sqrt(pi)*z**(x + 1/2)*exp(3*I*pi*(x + 1/2)/2)*gamma(x + 1/2)*hyper((-1/2, x + 1/2), (x + 3/2,), z*exp_polar(5*I*pi/2))/gamma(x + 3/2)
g.subs(x, 0).subs(z, 1).evalf(): 4.32257296433943 - 5.92192068901469*I
h.subs(x, 0).subs(z, 1).evalf(): 4.32257296433943 - 5.92192068901469*I

test1a_problematic
g: meijerg(((1, 2), ()), ((1/2,), (0,)), exp_polar(3*I*pi/2))
h: -4*sqrt(pi)*((1/2 - I/2)/sqrt(1 - exp_polar(I*pi/2)) + exp_polar(-I*pi/4)*asin(exp_polar(I*pi/4))/2)*exp(3*I*pi/4)
g.evalf(): 4.32257296433943 - 5.92192068901469*I
h.evalf(): 4.32257296433943 - 5.92192068901469*I

test1b
g: meijerg(((1, 2), ()), ((1/2,), (0,)), exp_polar(I*pi))
h: -I*pi**(3/2)
g.evalf(): -5.56832799683171*I
h.evalf(): -5.56832799683171*I

test2
g: meijerg(((0,), ()), ((1/2,), (1,)), exp_polar(I*pi/2))
h: (1/2 - I)*exp(-I)*exp(I*pi/4)
g.evalf(): 0.275572216790442 - 1.08354047147912*I
h.evalf(): 0.275572216790442 - 1.08354047147912*I

test2a
g: meijerg(((0,), ()), ((1/2,), (1,)), exp_polar(I*pi))
h: 3*E*I/2
g.evalf(): 4.07742274268857*I
h.evalf(): 4.07742274268857*I

test3
g: meijerg(((1,), ()), ((1,), ()), polar_lift(0))
h: 0
g.evalf(): 0
h.evalf(): 0

The results now are a little more consistent. But in the case of test1a_problematic (special case of test1_problematic) the results were consistent before (same with or without hyperexpand) and, even though they're still 'consistent', we've now introduced a sign change vs Wolfram. I wouldn't be suprised if there are cases where for c*z (c = 1, -1, I, -I) at least one c will always be problematic (for minor perturbations of this naive fix?)

However, in the case of @oscarbenjamin's z=0 example (test3), I do like the way the extra polar_lift handles things. Wolfram calculates that as 0 even though https://en.wikipedia.org/wiki/Meijer_G-function#Definition_of_the_Meijer_G-function has a z != 0 disclaimer with no discussion (?) of the z=0 case.

(it works because)

>>> unpolarify(polar_lift(0))
1

Unless there's an easy fix to polarified_meijerg_z (or something like trying multiple variations could be devised?) it doesn't make sense to add these changes in hyperexpand now. The z=0 case could be done but sympy already handles Integral(0) without meijerint!

There's also the problem of properly unpolarifying z after hyperexpand has completed. It's not so much a numerical correctness issue but e.g. if we just left z polar lifted like the above we'll get some uglified output e.g. test_holonomic:

assert hyperexpand(expr_to_holonomic(log(x)).to_meijerg()).simplify() == log(x)

will end up (depending on the approach) as

log(polar_lift(x - 1) + 1)

or even

Piecewise((log(-x) - I*pi, (x > 2) | (x < 0)), (log(x), True))

There is a halfway working solution to the uglification. If one unpolarifies z too early (e.g. in do_slater) we'll end up with the same results as in the tests "WITHOUT POLARIFY" above. But something along these lines (applied when hyperexpand is fully completed) seems to fix the introduced uglification in the test suite:

def deep_unpolarify_args(e):
    return e.replace(lambda f: isinstance(f, Function), 
                     lambda f: f.func(*[unpolarify(a) for a in f.args]))

^^ this may be totally invalid but is a bit better than just applying unpolarify to all args of any Expr (I think the warning about Add in help(polarify) and https://mathoverflow.net/questions/291622/defining-addition-on-the-riemann-surface-of-logz are relevant).

[ehren]

The above also doesn't work for meijerg(((0,), ()), ((1/2,), (1,)), -I) (even though it works in current sympy) but interestingly test1a_problematic is also z=-I. I seem to have a fix for the above using z.extract_multiplicatively(-I). It appears to work integrated in hyperexpand. Also the comments above re "uglification": I had the extra polarification also in the "hyper" case, now that it's only in the meijerg case the existing test suite is unaltered. Also "deep_unpolarify_args" (though no longer needed) is either wrong or unnecessary - it breaks (or maybe just uglifies) test_bessel in test_meijerint.

This doesn't fix the formula in the issue title. @skirpichev's commit seems to hint that's a different issue than the sign changes.

[ehren]

unpolarify(polar_lift(0))
1

I'm not sure how I arrived at this (maybe confused with polar_lift(1) == exp_polar(0)) but it is not what sympy actually produces:

>>> unpolarify(polar_lift(0))
0

(perhaps the z=0 case is a little more confusing than the others)