Merge branch 'u/rws/jacobi_p_polynomials_without_pexpect_maxima' of g…

…it://trac.sagemath.org/sage into t/22325/replace_all_pexpect_maxima_calls_in_function_orthogonal_polys

src/sage/functions/orthogonal_polys.py

@@ -1868,9 +1868,10 @@ def __init__(self):
 
 hermite = Func_hermite()
 
-def jacobi_P(n, a, b, x):
+
+class Func_jacobi_P(OrthogonalFunction):
     r"""
-    Returns the Jacobi polynomial `P_n^{(a,b)}(x)` for
+    Return the Jacobi polynomial `P_n^{(a,b)}(x)` for
     integers `n > -1` and a and b symbolic or `a > -1`
     and `b > -1`. The Jacobi polynomials are actually defined
     for all a and b. However, the Jacobi polynomial weight
@@ -1886,31 +1887,94 @@ def jacobi_P(n, a, b, x):
         sage: x = PolynomialRing(QQ, 'x').gen()
         sage: jacobi_P(2,0,0,x)
         3/2*x^2 - 1/2
-        sage: jacobi_P(2,1,2,1.2)        # random output of low order bits
-        5.009999999999998
+        sage: jacobi_P(2,1,2,1.2)
+        5.01000000000000
+    """
+    def __init__(self):
+        r"""
+        Init method for the Jacobi polynomials.
 
-    Check that :trac:`17192` is fixed::
+        EXAMPLES::
 
-        sage: x = PolynomialRing(QQ, 'x').gen()
-        sage: jacobi_P(0,0,0,x)
-        1
+            sage: _ = var('n a b x')
+            sage: loads(dumps(jacobi_P))
+            jacobi_P
+            sage: jacobi_P(n, a, b, x, hold=True)._sympy_()
+            jacobi(n, a, b, x)
+        """
+        OrthogonalFunction.__init__(self, "jacobi_P", nargs=4, latex_name=r"P",
+                conversions={'maxima':'jacobi_p', 'mathematica':'JacobiP',
+                    'maple':'JacobiP', 'sympy':'jacobi'})
 
-        sage: jacobi_P(-1,0,0,x)
-        Traceback (most recent call last):
-        ...
-        ValueError: n must be greater than -1, got n = -1
+    def _eval_(self, n, a, b, x):
+        """
+        EXAMPLES::
 
-        sage: jacobi_P(-7,0,0,x)
-        Traceback (most recent call last):
-        ...
-        ValueError: n must be greater than -1, got n = -7
-    """
-    if not (n > -1):
-        raise  ValueError("n must be greater than -1, got n = {0}".format(n))
+            sage: _ = var('n a b x')
+            sage: jacobi_P(1,n,n,n)
+            (n + 1)*n
+            sage: jacobi_P(2,n,n,n)
+            1/4*(2*n - 1)*(n + 2)*(n + 1)^2
+            sage: jacobi_P(1,n,n,x)
+            (n + 1)*x
+            sage: jacobi_P(3,2,1,x)
+            21/2*x^3 + 7/2*x^2 - 7/2*x - 1/2
+            sage: jacobi_P(1,a,b,x)
+            1/2*a*x + 1/2*b*x + 1/2*a - 1/2*b + x
+
+        TESTS:
 
-    _init()
-    return sage_eval(maxima.eval('jacobi_p(%s,%s,%s,x)'%(ZZ(n),a,b)), locals={'x':x})
+        Check that :trac:`17192` is fixed::
+
+            sage: x = PolynomialRing(QQ, 'x').gen()
+            sage: jacobi_P(0,0,0,x)
+            1
+            sage: jacobi_P(-1,0,0,x)
+            1
+            sage: jacobi_P(-1,1,1,x)
+            Traceback (most recent call last):
+            ...
+            ValueError: n must be greater than -1, got n = -1
+
+            sage: jacobi_P(-7,0,0,x)
+            231/16*x^6 - 315/16*x^4 + 105/16*x^2 - 5/16
+            sage: jacobi_P(-7,0,2,x)
+            Traceback (most recent call last):
+            ...
+            ValueError: n must be greater than -1, got n = -7
+        """
+        if SR(a).is_trivial_zero() and SR(b).is_trivial_zero():
+            return legendre_P(n, x)
+        if SR(n).is_numeric() and not (n > -1):
+            raise ValueError("n must be greater than -1, got n = {0}".format(n))
+        if not n in ZZ:
+            return
+        from sage.functions.other import gamma
+        s = sum(binomial(n,m) * gamma(a+b+n+m+1) / gamma(a+m+1) * ((x-1)/2)**m for m in range(n+1))
+        r = gamma(a+n+1) / factorial(n) / gamma(n+a+b+1) * s
+        return r.to_gamma().gamma_normalize()
+
+    def _evalf_(self, n, a, b, x, **kwds):
+        """
+        EXAMPLES::
+
+            sage: jacobi_P(2, 1, 2, 1.2)
+            5.01000000000000
+            sage: jacobi_P(2, 1, 2, 1.2, hold=True).n(20)
+            5.0100
+            sage: jacobi_P(2, 1, 2, pi+I, hold=True).n(100)
+            41.103034125334442891187112674 + 31.486722862692829003857755524*I
+        """
+        from sage.rings.complex_arb import ComplexBallField as CBF
+        the_parent = kwds.get('parent', None)
+        if the_parent is None:
+            the_parent = parent(x)
+        prec = the_parent.precision()
+        BF = CBF(prec+5)
+        ret = BF(x).jacobi_P(BF(n), BF(a), BF(b))
+        return the_parent(ret)
 
+jacobi_P = Func_jacobi_P()
 
 
 class Func_ultraspherical(GinacFunction):

src/sage/libs/pynac/pynac.pxd

@@ -144,7 +144,9 @@ cdef extern from "sage/libs/pynac/wrap.h":
     object py_object_from_numeric(GEx e)     except +
 
     # Algorithms
-    GEx g_gcd "gcd"(GEx a, GEx b) except +
+    GEx g_gcd "gcd"(GEx a, GEx b)   except +
+    GEx to_gamma(GEx expr)          except +
+    GEx gamma_normalize(GEx expr)   except +
 
     # Pattern matching wildcards
     GEx g_wild "wild"(unsigned int label) except +

src/sage/symbolic/expression.pyx

@@ -9614,6 +9614,53 @@ cdef class Expression(CommutativeRingElement):
 
     factorial_simplify = simplify_factorial
 
+    def to_gamma(self):
+        """
+        Convert factorial, binomial, and Pochhammer symbol
+        expressions to their gamma function equivalents.
+
+        EXAMPLES::
+
+            sage: m,n = var('m n', domain='integer')
+            sage: factorial(n).to_gamma()
+            gamma(n + 1)
+            sage: binomial(m,n).to_gamma()
+            gamma(m + 1)/(gamma(m - n + 1)*gamma(n + 1))
+        """
+        cdef GEx x
+        sig_on()
+        try:
+            x = to_gamma(self._gobj)
+        finally:
+            sig_off()
+        return new_Expression_from_GEx(self._parent, x)
+
+    def gamma_normalize(self):
+        """
+        Return the expression with any gamma functions that have
+        a common base converted to that base.
+
+        Addtionally the expression is normalized so any fractions
+        can be simplified through cancellation.
+
+        EXAMPLES::
+
+            sage: m,n = var('m n', domain='integer')
+            sage: (gamma(n+2)/gamma(n)).gamma_normalize()
+            (n + 1)*n
+            sage: (gamma(n+2)*gamma(n)).gamma_normalize()
+            (n + 1)*n*gamma(n)^2
+            sage: (gamma(n+2)*gamma(m-1)/gamma(n)/gamma(m+1)).gamma_normalize()
+            (n + 1)*n/((m - 1)*m)
+        """
+        cdef GEx x
+        sig_on()
+        try:
+            x = gamma_normalize(self._gobj)
+        finally:
+            sig_off()
+        return new_Expression_from_GEx(self._parent, x)
+
     def expand_sum(self):
         r"""
         For every symbolic sum in the given expression, try to expand it,