Error in integral result using hyper

[DouisLavid]

It seems that when computing the result of integrals that derive from the beta function, some unexpected factors appear :

Python 3.11.7 (main, Dec 15 2023, 18:12:31) [GCC 11.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import sympy as sp
>>> sp.__version__
'1.12'
>>> x = sp.Symbol("x", positive=True)
>>> y = sp.Symbol("y", positive=True)
>>> z = sp.Symbol("z", real=True)
>>> i = sp.Integral(z**(x-1)*(1-z)**(y-1), (z, 0, 1)) # = beta(x, y) by definition
>>> i.doit()
gamma(x)*hyper((x, 1 - y), (x + 1,), 1)/gamma(x + 1)
>>> _.simplify()
gamma(x)*gamma(y)/gamma(x + y) # correct
>>> i = sp.Integral(z**(x+0.5-1)*(1-z)**(y-1), (z, 0, 1)) # = beta(x+0.5, y)
>>> i.doit().simplify()
(-1)**y*exp(I*pi*y)*gamma(y)*gamma(x + 0.5)/gamma(x + y + 0.5) # ???
>>> 

Notice that using a particular value for y can yield results that are not consistent with the previous calculation :

>>> i = sp.Integral(z**(x-0.5)*sp.sqrt(1-z), (z, 0, 1)) # = beta(x+0.5, 1.5)
>>> i.doit().simplify()
-sqrt(pi)*gamma(x + 0.5)/(2*gamma(x + 2.0)) # only a sign error, recall gamma(3/2)=sqrt(pi)/2

[haru-44]

I will share what I have found out.

I don't know if this is the cause of all the issues, but it seems that meijerint_indefinite is not working as it is supposed to. When calculating sp.Integral(z**(x+sp.Rational(1,2)-1)*(1-z)**(sp.Rational(3,2)-1), (z, 0, 1)), it computes the following

i = meijerint_indefinite(z**(x - Rational(1, 2))*sqrt(1 - z), z)
i.subs(z,1) - i.subs(z,0)

Now, if we look at the internals of meijerint_indefinite, it stores the results in results and eventually returns one of them.

https://github.com/sympy/sympy/blob/82e59a73ae89887022c06f6d3b9361b5d7d02c6f/sympy/integrals/meijerint.py#L1687

In this example, we get two results, one of which gives us the correct answer, but actually returns the wrong one.

from sympy import *
from sympy.integrals.meijerint import _meijerint_indefinite_1, _find_splitting_points, _has

def meijerint_indefinite_part(f, x):
    """ part of meijerint_indefinite
    """
    f = sympify(f)
    results = []
    for a in sorted(_find_splitting_points(f, x) | {S.Zero}, key=default_sort_key):
        res = _meijerint_indefinite_1(f.subs(x, x + a), x)
        if not res:
            continue
        res = res.subs(x, x - a)
        if _has(res, hyper, meijerg):
            results.append(res)
        else:
            return res
    return results
    # return next(ordered(results))

x = Symbol("x", positive=True)
y = Symbol("y", positive=True)
z = Symbol("z", real=True)

i = meijerint_indefinite_part(z**(x - Rational(1, 2))*sqrt(1 - z), z)
(i[0].subs(z,1) - i[0].subs(z,0)).simplify() # sqrt(pi)*gamma(x + 1/2)/(2*gamma(x + 2))  # correct
(i[1].subs(z,1) - i[1].subs(z,0)).simplify() # -sqrt(pi)*gamma(x + 1/2)/(2*gamma(x + 2)) # bad