Integral with float exponent leads has incorrect sign

[oscarbenjamin]

From SO:
https://stackoverflow.com/questions/72223270/sympy-evaluates-integral-sign-wrong

With a float exponent 0.5 the integral comes out with the wrong sign:

In [11]: i = Integral((1-x**2)**y/(4*x), (x, 1/sqrt(2), 1))

In [12]: i
Out[12]: 
1              
⌠              
⎮          y   
⎮  ⎛     2⎞    
⎮  ⎝1 - x ⎠    
⎮  ───────── dx
⎮     4⋅x      
⌡              
√2             
──             
2              

In [13]: i.subs(y, .5).doit()
Out[13]: 
        ┌─  ⎛-0.5, -0.5 │  ⎞                        ┌─  ⎛-0.5, -0.5 │  ⎞
0.25⋅ⅈ⋅ ├─  ⎜           │ 1⎟ - 0.176776695296637⋅ⅈ⋅ ├─  ⎜           │ 2⎟
       2╵ 1 ⎝   0.5     │  ⎠                       2╵ 1 ⎝   0.5     │  ⎠

In [14]: N(_)    #  wrong sign:
Out[14]: -0.0435667014582489 - 2.68e-17⋅ⅈ

In [15]: i.subs(y, S.Half).doit()
Out[15]: 
  √2   acosh(√2)
- ── + ─────────
  8        4    

In [16]: N(_)
Out[16]: 0.0435667014582489

In [17]: i.subs(y, S.Half).evalf()
Out[17]: 0.0435667014582489

The integral seems to be evaluated by meijerg so I guess that the bug is in there somewhere.

[oscarbenjamin]

The bug seems to be in hyperexpand:

In [32]: G = -meijerg(((1.5,), (1,)), ((0, 0), ()), x**2*exp_polar(I*pi))

In [33]: F = hyperexpand(G)

In [34]: G
Out[34]: 
 ╭─╮2, 1 ⎛1.5   1 │  2  ⅈ⋅π⎞
-│╶┐     ⎜        │ x ⋅ℯ   ⎟
 ╰─╯2, 2 ⎝0, 0    │        ⎠

In [35]: F
Out[35]: 
                              ⎛           │ 1 ⎞
                     1.0  ┌─  ⎜-0.5, -0.5 │ ──⎟
-7.08981540362206⋅ⅈ⋅x   ⋅ ├─  ⎜           │  2⎟
                         2╵ 1 ⎝   0.5     │ x ⎠

In [36]: G.subs(x, 0.5).n()
Out[36]: 3.19702813586035 - 11.1366559936634⋅ⅈ

In [37]: F.subs(x, 0.5).n()
Out[37]: -3.19702813586035 - 11.1366559936634⋅ⅈ

When the input to the G function is rational a different expansion is found:

In [41]: Gr = -meijerg(((Rational(3, 2),), (1,)), ((0, 0), ()), x**2*exp_polar(I*pi))

In [42]: Gr
Out[42]: 
 ╭─╮2, 1 ⎛3/2   1 │  2  ⅈ⋅π⎞
-│╶┐     ⎜        │ x ⋅ℯ   ⎟
 ╰─╯2, 2 ⎝0, 0    │        ⎠

In [43]: print(hyperexpand(Gr))
-4*I*sqrt(pi)*x*(((1/4 - 1/(4*x**2))*(1 - 1/x**2) + (1/2 - 1/(2*x**2))*(1 + 1/(2*x**2)) + 1/4 - 1/(4*x**2))*Piecewise((I/sqrt(-1 + x**(-2)), 1/Abs(x**2) > 1), (1/sqrt(1 - 1/x**2), True)) + Piecewise((I*acosh(sqrt(x**(-2)))/sqrt(x**(-2)) + pi/(2*sqrt(x**(-2))), 1/Abs(x**2) > 1), (asin(sqrt(x**(-2)))/sqrt(x**(-2)), True))/x**2)

This expansion is undefined at x=1 but taking limits gives the right results:

In [44]: N(Gr.subs(x, 1) - Gr.subs(x, 1/sqrt(2)))
Out[44]: -1.23551948433479 - 0.e-20⋅ⅈ

In [45]: Fr = hyperexpand(Gr)

In [46]: N(Fr.subs(x, 1) - Fr.subs(x, 1/sqrt(2)))
Out[46]: nan

In [47]: N(Fr.limit(x, 1) - Fr.limit(x, 1/sqrt(2)))
Out[47]: -1.23551948433479 + 0.e-21⋅ⅈ

Presumably that's how the integral succeeds when the exponent is rational.

I think the problem is maybe somewhere around here:

sympy/sympy/simplify/hyperexpand.py

Lines 2309 to 2316 in 68c37df

	k = polar_lift(S.NegativeOne**(len(ap) - len(bm)))
	harg = k*zfinal
	# NOTE even though k "is" +-1, this has to be t/k instead of
	# t*k ... we are using polar numbers for consistency!
	premult = (t/k)**bh
	hyp = _hyperexpand(Hyper_Function(nap, nbq), harg, ops,
	t, premult, bh, rewrite=None)
	res += fac * hyp

At least it is that call to _hyperexpand that seems to return different things depending on whether the input is float or rational.

[oscarbenjamin]

The lookup for inv in self.concreteformulae[sizes] fails here:

sympy/sympy/simplify/hyperexpand.py

Lines 853 to 855 in 68c37df

	if sizes in self.concrete_formulae and \
	inv in self.concrete_formulae[sizes]:
	return self.concrete_formulae[sizes][inv]

In the debugger:

(Pdb) p inv
(0, ((0.500000000000000, 2),), ((0.500000000000000, 1),))
(Pdb) p sizes
(2, 1)
(Pdb) p list(self.concrete_formulae[sizes])
[(0, ((0, 2),), ((0, 1),)), (0, ((0, 1), (1/2, 1)), ((1/2, 1),)), (0, ((1/2, 2),), ((1/2, 1),)), (0, ((0, 2),), ((1/2, 1),)), (0, ((1/2, 2),), ((0, 1),))]
(Pdb) p inv in self.concrete_formulae[sizes]
False
(Pdb) p inv in list(self.concrete_formulae[sizes])
True

The dict list lookup succeeds and the float inv compares equal to the entry (0, ((1/2, 2),), ((1/2, 1),)). However the dict lookup fails. That's because Float.__hash__ doesn't work (#11707). This is the cause of the difference in behaviour between rational and float.

However the fallback case when a formula isn't found should still be able to return a result with the correct sign. The fallback path goes through build_hypergeometric_formula which tries to rewrite a G-function in terms of an F-function.

[oscarbenjamin]

This is a way to reproduce the bug without using floats:

In [74]: G = meijerg(((Rational(1, 2),), (1,)), ((0, 0), ()), x*exp_polar(I*pi))

In [75]: F = hyperexpand(G)

In [76]: G
Out[76]: 
╭─╮2, 1 ⎛1/2   1 │    ⅈ⋅π⎞
│╶┐     ⎜        │ x⋅ℯ   ⎟
╰─╯2, 2 ⎝0, 0    │       ⎠

In [77]: F
Out[77]: 
         ┌─  ⎛1/2, 1/2 │ 1⎞ 
-2⋅ⅈ⋅√π⋅ ├─  ⎜         │ ─⎟ 
        2╵ 1 ⎝  3/2    │ x⎠ 
────────────────────────────
             √x             

In [78]: G.subs(x, -1).n()
Out[78]: 3.12438801679839

In [79]: F.subs(x, -1).n()
Out[79]: -3.12438801679839

[oscarbenjamin]

This is a way to reproduce the bug without using floats:

Actually I realise I used this diff when running that:

diff --git a/sympy/simplify/hyperexpand.py b/sympy/simplify/hyperexpand.py
index aa0a534..5f93211 100644
--- a/sympy/simplify/hyperexpand.py
+++ b/sympy/simplify/hyperexpand.py
@@ -850,9 +850,9 @@ def lookup_origin(self, func):
         """
         inv = func.build_invariants()
         sizes = func.sizes
-        if sizes in self.concrete_formulae and \
-                inv in self.concrete_formulae[sizes]:
-            return self.concrete_formulae[sizes][inv]
+        #if sizes in self.concrete_formulae and \
+        #        inv in self.concrete_formulae[sizes]:
+        #    return self.concrete_formulae[sizes][inv]
 
         # We don't have a concrete formula. Try to instantiate.
         if sizes not in self.symbolic_formulae:

Without the diff hyperexpand will remove the hypergeometric functions but the result still has the wrong sign (at least for x=-1:

In [1]: G = meijerg(((Rational(1, 2),), (1,)), ((0, 0), ()), x*exp_polar(I*pi))

In [2]: F = hyperexpand(G)

In [3]: G
Out[3]: 
╭─╮2, 1 ⎛1/2   1 │    ⅈ⋅π⎞
│╶┐     ⎜        │ x⋅ℯ   ⎟
╰─╯2, 2 ⎝0, 0    │       ⎠

In [4]: F
Out[4]: 
        ⎛⎧       ⎛    ___⎞                         ⎞ 
        ⎜⎪       ⎜   ╱ 1 ⎟                         ⎟ 
        ⎜⎪ⅈ⋅acosh⎜  ╱  ─ ⎟                         ⎟ 
        ⎜⎪       ⎝╲╱   x ⎠       π           1     ⎟ 
        ⎜⎪──────────────── + ─────────  for ─── > 1⎟ 
        ⎜⎪        ___              ___      │x│    ⎟ 
        ⎜⎪       ╱ 1              ╱ 1              ⎟ 
        ⎜⎪      ╱  ─         2⋅  ╱  ─              ⎟ 
        ⎜⎪    ╲╱   x           ╲╱   x              ⎟ 
-2⋅ⅈ⋅√π⋅⎜⎨                                         ⎟ 
        ⎜⎪           ⎛    ___⎞                     ⎟ 
        ⎜⎪           ⎜   ╱ 1 ⎟                     ⎟ 
        ⎜⎪       asin⎜  ╱  ─ ⎟                     ⎟ 
        ⎜⎪           ⎝╲╱   x ⎠                     ⎟ 
        ⎜⎪       ─────────────           otherwise ⎟ 
        ⎜⎪              ___                        ⎟ 
        ⎜⎪             ╱ 1                         ⎟ 
        ⎜⎪            ╱  ─                         ⎟ 
        ⎝⎩          ╲╱   x                         ⎠ 
─────────────────────────────────────────────────────
                          √x                         

In [5]: G.subs(x, -1).n()
Out[5]: 3.12438801679839

In [6]: F.subs(x, -1).n()
Out[6]: -3.12438801679839

[jksuom]

These issues are caused by the multiple values of G-functions on branch cuts. Most G-functions are branched at 0 and infinity, and also at either 1 or -1 if the numbers of a-parameters and b-parameters are the same.

I have not studied the code in detail but it seems that mpmath implements G-functions with branch cuts along the negative real axis and also along (1, oo) if 1 is a branch point. mpmath will return one of the boundary values but that may be the wrong one.

SymPy has a (partial) solution to the problem: HyperRep functions whose arguments are polar numbers. Using exp_polar(I*pi) or exp_polar(-I*pi) makes it possible to choose the desired boundary value on the branch cut at -1. More generally, an argument like x*exp_polar(I*pi) will be resolved by rewrite('nonrep') to a branch with branch cut along the positive real axis (i.e., for negative real x).

Representations by HyperRep functions are currently implemented only for exact parameter values. There is no polar number support for plain hyper functions. It seems that the issue of OP would be solved if integer and half-integer Float parameters would be converted to exact ones in _hyperexpand. Otherwise, it seems that plain hyper functions created by do_slater should not be analytically continued as in this block:

sympy/sympy/simplify/hyperexpand.py

Lines 2397 to 2408 in 918eb81

	if m.delta > 0 or \
	(m.delta == 0 and len(m.ap) == len(m.bq) and
	(re(m.nu) < -1) is not False and polar_lift(z0) == polar_lift(1)):
	# The condition delta > 0 means that the convergence region is
	# connected. Any expression we find can be continued analytically
	# to the entire convergence region.
	# The conditions delta==0, p==q, re(nu) < -1 imply that G is continuous
	# on the positive reals, so the values at z=1 agree.
	if cond1 is not False:
	cond1 = True
	if cond2 is not False:
	cond2 = True

One way to implement that seems to be this:

 --- a/sympy/simplify/hyperexpand.py
+++ b/sympy/simplify/hyperexpand.py
@@ -2402,9 +2402,9 @@ def tr(l):
         # to the entire convergence region.
         # The conditions delta==0, p==q, re(nu) < -1 imply that G is continuous
         # on the positive reals, so the values at z=1 agree.
-        if cond1 is not False:
+        if cond1 is not False and not slater1.has(hyper):
             cond1 = True
-        if cond2 is not False:
+        if cond2 is not False and not slater2.has(hyper):
             cond2 = True
 
     if cond1 is True:

[oscarbenjamin]

One way to implement that seems to be this:

Does that just prevent hyperexpand from converting a G function into anything involving hyper?

[jksuom]

No, it just prevents the hyper components from being analytically continued outside the unit disk in case the continuation is not unique, The hyper components may be included in the output by setting allow_hyper=True.

[oscarbenjamin]

That seems reasonable. To be clear is the rationale for that something like:

With MeijerG it is possible to control the desired branch using exp_polar but that does not work with hyper since it has no support for exp_polar. That means that under certain conditions conversion from MeijerG to hyper would lose information about which branch was intended and potentially result in selecting a different branch so that the hyper function would not be equivalent to the original MeijerG function. To be safe we prevent the conversion from MeijerG to hyper unless the intended branch is unambiguous.

mpmath will return one of the boundary values but that may be the wrong one.

The problem as described above seems to be purely to do with sympy's hyperexpand. Is there also a problem with mpmath's definition/implementation of MeijerG or hyper?

[jksuom]

Is there also a problem with mpmath's definition/implementation of MeijerG or hyper?

G-functions are defined by contour integrals over s of z**s = exp(s*log(z)) multiplied by a meromorphic function of s (a 'gamma product'). It is analytic in z apart from some singularities and multivalued with {-oo, 0) as a branch cut (assuming that the principal value of log(z) is used).

mpmath does not implement G-functions by means of these integrals; instead, it computes their values by means of (generalized) hypergeometric functions (hyper) as explained here. Actually, the description is slightly inaccurate: a G-function is represented by a sum of hypergeometric series multiplied by a power of the argument z. Such a product is branched at 0 and at infinity with branch cut {-oo, 0) if the exponent of the power is not an integer.

If the number p of a-parameters of a G-function is different from the number q of its b-parameters, then (the principal branch of) the function is well defined with branch cut (-oo, 0).
If p < q, then the contour defining the integral representation is a loop starting and ending at +oo. The hypergeometric expansion consists of power series in z or -z with infinite radius of convergence (series=1 in mpmath).
If p > q, then the contour is a loop starting and ending at -oo. The expansion consists of power series in 1/z or -1/z (series=2 in mpmath).

If p = q, there is no single contour defining (the principal branch of) a G-function, in general. The integral over the loop passing through +oo converges if |z| < 1, and can be represented by hyper functions of argument z or -z. The integral over the loop passing through -oo converges if |z |> 1, and can be represented by hyper functions of argument 1/z or -1/z. The series representing these hyper functions converge in the unit disc and have a branch point at 1. They can be analytically continued to the complex plane with branch cut along (1, oo). Hence there are three possibilities for defining the principal branch of a G-function with p = q:

A. Use the series converging in the unit disk and extend them analytically to the whole plane. There will be a branch cut along (1, oo), if the argument of the hyper functions is z. If the argument is -z, the cut from -oo to -1 will be covered by (-oo, 0) if 0 is a branch point. This is the series=1 case of mpmath.

B. Use the series converging in the exterior of the unit disk and extend then analytically to the disk. This will introduce a possible branch cut (0,1) (extending (-oo, 0)). This is the series=2 case of mpmath.

C. Use the series converging in the unit disk if |z| < 1 and the series converging in its exterior if |z| > 1. There will be branch cuts on the unit circle if 1 is a branch point.

From the code it is seen that mpmath favors A in case p = q. B is not used and C is only used if m+n == p. This condition says that delta == 0, which means that the integral along the contour from -oo*I to oo*I is convergent at most on the positive real axis (|arg(z)| <= delta*pi). I do not understand the role of that condition here. (Maybe that is coming from Mathematica.) In any case, it seems that these choices might not be optimal for SymPy.

EDIT: Link to mpmath code corrected.