Ticket #2516: hyper3.patch

doc/en/reference/functions/index.rst

@@ -11 +11 @@
    sage/functions/orthogonal_polys
    sage/functions/other
    sage/functions/special
+   sage/functions/hypergeometric
    sage/functions/exp_integral
    sage/functions/wigner
    sage/functions/generalized

sage/calculus/calculus.py

@@ -1685 +1685 @@
 
 maxima_qp = re.compile("\?\%[a-z|A-Z|0-9|_]*")  # e.g., ?%jacobi_cd
 
-maxima_var = re.compile("\%[a-z|A-Z|0-9|_]*")  # e.g., ?%jacobi_cd
+maxima_var = re.compile("\%[a-z|A-Z|0-9|_]*")  # e.g., %jacobi_cd
 
 sci_not = re.compile("(-?(?:0|[1-9]\d*))(\.\d+)?([eE][-+]\d+)")
 
@@ -1693 +1693 @@
 
 maxima_polygamma = re.compile("psi\[(\d*)\]\(")  # matches psi[n]( where n is a number
 
+maxima_hyper = re.compile("\%f\[\d+,\d+\]")  # matches %f[m,n]
+
 def symbolic_expression_from_maxima_string(x, equals_sub=False, maxima=maxima):
     """
     Given a string representation of a Maxima expression, parse it and
@@ -1778 +1780 @@
                 pass
             else:
                 syms[X[2:]] = function_factory(X[2:])
-        s = s.replace("?%","")
+        s = s.replace("?%", "")
 
-    s = polylog_ex.sub('polylog(\\1,',s)
+    s = maxima_hyper.sub('hypergeometric', s)
+
+    s = polylog_ex.sub('polylog(\\1,', s)
     s = multiple_replace(symtable, s)
     s = s.replace("%","")
 
     s = s.replace("#","!=") # a lot of this code should be refactored somewhere...
 
-    s = maxima_polygamma.sub('psi(\g<1>,',s) # this replaces psi[n](foo) with psi(n,foo), ensuring that derivatives of the digamma function are parsed properly below
+    s = maxima_polygamma.sub('psi(\g<1>,', s) # this replaces psi[n](foo) with psi(n,foo), ensuring that derivatives of the digamma function are parsed properly below
 
     if equals_sub:
         s = s.replace('=','==')

sage/functions/all.py

@@ -67 +67 @@
                           sin_integral, cos_integral, Si, Ci,
                           sinh_integral, cosh_integral, Shi, Chi,
                           exponential_integral_1, Ei)
+
+from hypergeometric import hypergeometric

new file sage/functions/hypergeometric.py

@@ +1 @@
+r"""
+Hypergeometric Functions
+
+This module implements manipulation of infinite hypergeometric series
+represented in standard parametric form (as $\,_pFq$ functions).
+
+AUTHORS:
+
+- Fredrik Johansson (2010): initial version
+
+- Eviatar Bach (2013): major changes
+
+EXAMPLES:
+
+Examples from :trac:`9908`::
+
+    sage: maxima('integrate(bessel_j(2, x), x)').sage()
+    1/24*x^3*hypergeometric((3/2,), (5/2, 3), -1/4*x^2)
+    sage: sum(((2 * I)^x / (x^3 + 1) * (1/4)^x), x, 0, oo)
+    hypergeometric((1, 1, -1/2*I*sqrt(3) - 1/2, 1/2*I*sqrt(3) - 1/2), (2, -1/2*I*sqrt(3) + 1/2, 1/2*I*sqrt(3) + 1/2), 1/2*I)
+    sage: sum((-1)^x / ((2 * x + 1) * factorial(2 * x + 1)), x, 0, oo)
+    hypergeometric((1/2,), (3/2, 3/2), -1/4)
+
+Simplification::
+
+    sage: hypergeometric((hypergeometric((), (), x),), (),
+    ....:                hypergeometric((), (), x)).simplify_hypergeometric()
+    1/((-e^x + 1)^e^x)
+    sage: hypergeometric([-2], [], x).simplify_hypergeometric()
+    x^2 - 2*x + 1
+    sage: hypergeometric([], [], x).simplify_hypergeometric()
+    e^x
+    sage: hypergeometric((hypergeometric((), (), x),), (),
+    ....:                hypergeometric((), (), x)).simplify_hypergeometric(algorithm='sage')
+    (-e^x + 1)^(-e^x)
+
+Equality testing::
+
+    sage: bool(hypergeometric([], [], x).derivative(x) ==
+    ....:      hypergeometric([], [], x))  # diff(e^x, x) == e^x
+    True
+    sage: bool(hypergeometric([], [], x) == hypergeometric([], [1], x))
+    False
+
+Computing terms and series::
+
+    sage: z = var('z')
+    sage: hypergeometric([], [], z).series(z, 0)
+    Order(1)
+    sage: hypergeometric([], [], z).series(z, 1)
+    1 + Order(z)
+    sage: hypergeometric([], [], z).series(z, 2)
+    1 + 1*z + Order(z^2)
+    sage: hypergeometric([], [], z).series(z, 3)
+    1 + 1*z + 1/2*z^2 + Order(z^3)
+
+    sage: hypergeometric([-2], [], z).series(z, 3)
+    1 + (-2)*z + 1*z^2
+    sage: hypergeometric([-2], [], z).series(z, 6)
+    1 + (-2)*z + 1*z^2
+    sage: hypergeometric([-2], [], z).series(z, 6).is_terminating_series()
+    True
+    sage: hypergeometric([-2], [], z).series(z, 2)                        
+    1 + (-2)*z + Order(z^2)
+    sage: hypergeometric([-2], [], z).series(z, 2).is_terminating_series()
+    False
+
+    sage: hypergeometric([1], [], z).series(z, 6) 
+    1 + 1*z + 1*z^2 + 1*z^3 + 1*z^4 + 1*z^5 + Order(z^6)
+    sage: hypergeometric([], [1/2], -z^2/4).series(z, 11)
+    1 + (-1/2)*z^2 + 1/24*z^4 + (-1/720)*z^6 + 1/40320*z^8 + (-1/3628800)*z^10 + Order(z^11)
+
+    sage: hypergeometric([1],[5],x).series(x,5)
+    1 + 1/5*x + 1/30*x^2 + 1/210*x^3 + 1/1680*x^4 + Order(x^5)
+
+    sage: sum(hypergeometric([1, 2], [3], 1/3).terms(6)).n()
+    1.29788359788360
+    sage: hypergeometric([1, 2], [3], 1/3).n()
+    1.29837194594696
+    sage: hypergeometric([], [], x).series(x, 20)(x=1).n() == e.n()
+    True
+
+Plotting::
+
+    sage: plot(hypergeometric([1, 1], [3, 3, 3], x), x, -30, 30)
+
+Numeric evaluation::
+
+    sage: hypergeometric([1], [], 1/10).n()  # geometric series
+    1.11111111111111
+    sage: hypergeometric([], [], 1).n()  # e
+    2.71828182845905
+    sage: hypergeometric([], [], 3., hold=True)
+    hypergeometric((), (), 3.00000000000000)
+    sage: hypergeometric([1,2,3],[4,5,6],1/2).n()
+    1.02573619590134
+    sage: hypergeometric([1,2,3],[4,5,6],1/2).n(digits=30)
+    1.02573619590133865036584139535
+    sage: hypergeometric([5-3*I],[3/2,2+I,sqrt(2)],4+I).n() #TOLERANCE?
+    5.52605111678805 - 7.86331357527544*I
+
+SymPy conversion::
+
+    sage: hypergeometric((5,4), (4,4), 3)._sympy_()
+    hyper((5, 4), (4, 4), 3)
+
+"""
+#*****************************************************************************
+#       Copyright (C) 2010 Fredrik Johansson <fredrik.johansson@gmail.com>
+#       Copyright (C) 2013 Eviatar Bach <eviatarbach@gmail.com>
+#
+#  Distributed under the terms of the GNU General Public License (GPL)
+#  as published by the Free Software Foundation; either version 2 of
+#  the License, or (at your option) any later version.
+#                  http://www.gnu.org/licenses/
+#*****************************************************************************
+
+from sage.rings.all import (Integer, ZZ, QQ, Infinity, binomial,
+                            rising_factorial, factorial)
+from sage.functions.other import sqrt, gamma, real_part
+from sage.functions.log import exp, log
+from sage.functions.trig import cos, sin
+from sage.functions.hyperbolic import cosh, sinh
+from sage.functions.other import erf
+from sage.symbolic.constants import pi
+from sage.symbolic.all import I
+from sage.symbolic.function import BuiltinFunction, is_inexact
+from sage.symbolic.ring import SR
+from sage.structure.element import get_coercion_model 
+from sage.misc.latex import latex
+from sage.misc.flatten import flatten
+from sage.misc.misc_c import prod
+from mpmath import hyper
+from sage.libs.mpmath import utils as mpmath_utils
+from sage.symbolic.expression import Expression
+from sage.interfaces.maxima import maxima
+from sage.calculus.functional import derivative
+
+def rational_param_as_tuple(x):
+    """
+    Utility function for converting rational PFQ parameters to
+    tuples (which mpmath handles more efficiently).
+
+        sage: from sage.functions.hypergeometric import rational_param_as_tuple
+        sage: rational_param_as_tuple(1/2)
+        (1, 2)
+        sage: rational_param_as_tuple(3)
+        3
+        sage: rational_param_as_tuple(pi)
+        pi
+
+    """
+    try:
+        x = x.pyobject()
+    except AttributeError:
+        pass
+    try:
+        if x.parent() is QQ:
+            p = int(x.numer())
+            q = int(x.denom())
+            return p, q
+    except AttributeError:
+        pass
+    return x
+
+
+class Hypergeometric(BuiltinFunction):
+    r"""
+    Represents a (formal) generalized infinite hypergeometric series.
+    """
+
+    def __init__(self):
+        """
+        EXAMPLES::
+
+            sage: hypergeometric([], [], 1)
+            hypergeometric((), (), 1)
+            sage: hypergeometric([], [1], 1)
+            hypergeometric((), (1,), 1)
+            sage: hypergeometric([2, 3], [1], 1)
+            hypergeometric((2, 3), (1,), 1)
+
+        """
+        BuiltinFunction.__init__(self, 'hypergeometric', nargs=3,
+                                 conversions={'mathematica':'HypergeometricPFQ',
+                                              'maxima':'hypergeometric',
+                                              'sympy':'hyper'})
+
+    def __call__(self, a, b, z, **kwargs):
+        return BuiltinFunction.__call__(self,
+                                        SR._force_pyobject(a),
+                                        SR._force_pyobject(b),
+                                        z, **kwargs)
+
+    def _print_latex_(self, a, b, z):
+        r"""
+        TESTS::
+
+            sage: latex(hypergeometric([1, 1], [2], -1))
+            \,_2F_1\left(\begin{matrix} 1,1 \\ 2 \end{matrix} ; -1 \right)
+
+        """
+        aa = ",".join(latex(c) for c in a)
+        bb = ",".join(latex(c) for c in b)
+        z = latex(z)
+        return (r"\,_{}F_{}\left(\begin{{matrix}} {} \\ {} \end{{matrix}} ; "
+                r"{} \right)").format(len(a), len(b), aa, bb, z)
+
+    def _eval_(self, a, b, z, **kwargs):
+        coercion_model = get_coercion_model()
+        co = reduce(lambda x, y: coercion_model.canonical_coercion(x, y)[0],
+                    a + b + (z,))
+        if is_inexact(co) and not isinstance(co, Expression):
+            from sage.structure.coerce import parent
+            return self._evalf_(a, b, z, parent=parent(co))
+        if z == 0:
+            return Integer(1)
+        return
+
+    def _evalf_(self, a, b, z, parent):
+        """
+        TESTS::
+
+            sage: hypergeometric([1, 1], [2], -1).n()
+            0.693147180559945
+            sage: hypergeometric([], [], RealField(100)(1))
+            2.7182818284590452353602874714
+
+        """
+        aa = [rational_param_as_tuple(c) for c in a]
+        bb = [rational_param_as_tuple(c) for c in b]
+        return mpmath_utils.call(hyper, aa, bb, z, parent=parent)
+
+    def _tderivative_(self, a, b, z, *args, **kwargs):
+        """
+        EXAMPLES::
+
+            sage: hypergeometric([1/3, 2/3], [5], x^2).diff(x)
+            4/45*x*hypergeometric((4/3, 5/3), (6,), x^2)
+            sage: hypergeometric([1, 2], [x], 2).diff(x)
+            Traceback (most recent call last):
+                ...
+            NotImplementedError: derivative of hypergeometric function with respect to parameters. Try calling `.simplify_hypergeometric()` first.
+        """
+        diff_param = kwargs['diff_param']
+        if diff_param in hypergeometric(a, b, 1).variables():
+            raise NotImplementedError("derivative of hypergeometric function "
+                                      "with respect to parameters. Try calling"                                      " `.simplify_hypergeometric()` first.")
+        t = (reduce(lambda x, y: x * y, a, 1) *
+             reduce(lambda x, y: x / y, b, Integer(1)))
+        return (t * derivative(z, diff_param) *
+                hypergeometric([c + 1 for c in a], [c + 1 for c in b], z))
+
+    '''def _maxima_init_evaled_(self, a, b, z):
+        return "hypergeometric({},{},{})".format(repr(maxima(a)),
+                                                 repr(maxima(b)),
+                                                 repr(maxima(z)))
+
+    def _maxima_lib_(self, a, b, z):
+        return "hypergeometric({},{},{})".format(repr(maxima(a.operands())),
+                                                 repr(maxima(b.operands())),
+                                                 repr(maxima(z)))'''
+
+    class EvaluationMethods:
+        '''def operands(self, ex, a, b, z):
+            return [(a), list(b), z]'''
+
+        def sorted_parameters(self, ex, a, b, z):
+            """
+            Return with parameters sorted in a canonical order
+
+            EXAMPLES::
+
+                sage: hypergeometric([2,1,3], [5,4], 1/2).sorted_parameters()
+                hypergeometric((1, 2, 3), (4, 5), 1/2)
+
+            """
+            return hypergeometric(sorted(a), sorted(b), z)
+
+        def eliminate_parameters(self, ex, a, b, z):
+            """
+            Eliminate repeated parameters by pairwise cancellation of identical
+            terms in ``a`` and ``b``
+
+            EXAMPLES::
+
+                sage: hypergeometric([1, 1, 2, 5], [5, 1, 4],
+                ....:                1/2).eliminate_parameters()
+                hypergeometric((1, 2), (4,), 1/2)
+                sage: hypergeometric([x], [x], x).eliminate_parameters()
+                hypergeometric((), (), x)
+                sage: hypergeometric((5, 4), (4, 4), 3).eliminate_parameters()
+                hypergeometric((5,), (4,), 3)
+
+            """
+            aa = list(a)
+            bb = list(b)
+            p = pp = len(aa)
+            q = qq = len(bb)
+            i = 0
+            while i < qq and aa:
+                bbb = bb[i]
+                if bbb in aa:
+                    aa.remove(bbb)
+                    bb.remove(bbb)
+                    pp -= 1
+                    qq -= 1
+                else:
+                    i += 1
+            if (pp, qq) != (p, q):
+                return hypergeometric(aa, bb, z)
+            return ex
+
+        def is_termwise_finite(self, ex, a, b, z):
+            """
+            Determine whether all terms of self are finite. Any infinite
+            terms or ambiguous terms beyond the first zero, if one exists,
+            are ignored.
+
+            Ambiguous cases (where a term is the product of both zero
+            and an infinity) are not considered finite.
+
+
+                sage: hypergeometric([2], [3, 4], 5).is_termwise_finite()
+                True
+                sage: hypergeometric([2], [-3, 4], 5).is_termwise_finite()
+                False
+                sage: hypergeometric([-2], [-3, 4], 5).is_termwise_finite()
+                True
+                sage: hypergeometric([-3], [-3, 4], 5).is_termwise_finite()  # ambiguous
+                False
+
+                sage: hypergeometric([0], [-1], 5).is_termwise_finite()
+                True
+                sage: hypergeometric([0], [0], 5).is_termwise_finite()  # ambiguous
+                False
+                sage: hypergeometric([1], [2], Infinity).is_termwise_finite()
+                False
+                sage: hypergeometric([0], [0], Infinity).is_termwise_finite()  # ambiguous
+                False
+                sage: hypergeometric([0], [], Infinity).is_termwise_finite()  # ambiguous
+                False
+
+            """
+            if z == 0:
+                return 0 not in b
+            if abs(z) == Infinity:
+                return False
+            if abs(z) == Infinity:
+                return False
+            for bb in b:
+                if bb in ZZ and bb <= 0:
+                    if any((aa in ZZ) and (bb < aa <= 0) for aa in a):
+                        continue
+                    return False
+            return True
+
+        def is_terminating(self, ex, a, b, z):
+            """
+            Determine whether the series represented by self terminates
+            after a finite number of terms, i.e. whether any of the
+            numerator parameters are nonnegative integers (with no
+            preceding nonnegative denominator parameters), or z = 0.
+
+            If terminating, the series represents a polynomial of z
+
+            EXAMPLES::
+
+                sage: hypergeometric([1, 2], [3, 4], x).is_terminating()
+                False
+                sage: hypergeometric([1, -2], [3, 4], x).is_terminating()
+                True
+                sage: hypergeometric([1, -2], [], x).is_terminating()
+                True
+
+            """
+            if z == 0:
+                return True
+            for aa in a:
+                if (aa in ZZ) and (aa <= 0):
+                    return ex.is_termwise_finite()
+            return False
+
+        def is_absolutely_convergent(self, ex, a, b, z):
+            """
+            Determine whether self converges absolutely as an infinite series.
+            False is returned if not all terms are finite.
+
+            EXAMPLES:
+
+            Degree giving infinite radius of convergence::
+
+                sage: hypergeometric([2,3],[4,5],6).is_absolutely_convergent()
+                True
+                sage: hypergeometric([2,3],[-4,5],6).is_absolutely_convergent()  # undefined
+                False
+                sage: hypergeometric([2,3],[-4,5],Infinity).is_absolutely_convergent()  # undefined
+                False
+
+            Ordinary geometric series (unit radius of convergence)::
+
+                sage: hypergeometric([1],[],1/2).is_absolutely_convergent()
+                True
+                sage: hypergeometric([1],[],2).is_absolutely_convergent()
+                False
+                sage: hypergeometric([1],[],1).is_absolutely_convergent()
+                False
+                sage: hypergeometric([1],[],-1).is_absolutely_convergent()
+                False
+                sage: hypergeometric([1],[],-1).n()    # Sum still exists
+                0.500000000000000
+
+            Degree $p = q+1$ (unit radius of convergence)::
+
+                sage: hypergeometric([2,3],[4],6).is_absolutely_convergent()
+                False
+                sage: hypergeometric([2,3],[4],1).is_absolutely_convergent()
+                False
+                sage: hypergeometric([2,3],[5],1).is_absolutely_convergent()
+                False
+                sage: hypergeometric([2,3],[6],1).is_absolutely_convergent()
+                True
+                sage: hypergeometric([-2,3],[4],5).is_absolutely_convergent()
+                True
+                sage: hypergeometric([2,-3],[4],5).is_absolutely_convergent()
+                True
+                sage: hypergeometric([2,-3],[-4],5).is_absolutely_convergent()
+                True
+                sage: hypergeometric([2,-3],[-1],5).is_absolutely_convergent()
+                False
+
+            Degree giving zero radius of convergence::
+
+                sage: hypergeometric([1,2,3],[4],2).is_absolutely_convergent()
+                False
+                sage: hypergeometric([1,2,3],[4],1/2).is_absolutely_convergent()
+                False
+                sage: hypergeometric([1,2,-3],[4],1/2).is_absolutely_convergent()  # polynomial
+                True
+
+            """
+            p, q = len(a), len(b)
+            if not ex.is_termwise_finite():
+                return False
+            if p <= q:
+                return True
+            if ex.is_terminating():
+                return True
+            if p == q + 1:
+                if abs(z) < 1:
+                    return True
+                if abs(z) == 1:
+                    if real_part(sum(b) - sum(a)) > 0:
+                        return True
+            return False
+
+        def terms(self, ex, a, b, z, n=None):
+            """
+            Generate the terms of self (optionally only n terms).
+
+            EXAMPLES::
+
+                sage: list(hypergeometric([-2, 1], [3, 4], x).terms())
+                [1, -1/6*x, 1/120*x^2]
+                sage: list(hypergeometric([-2, 1], [3, 4], x).terms(2))
+                [1, -1/6*x]
+                sage: list(hypergeometric([-2, 1], [3, 4], x).terms(0))
+                []
+
+            """
+            if n is None:
+                n = Infinity
+            t = Integer(1)
+            k = 1
+            while k <= n:
+                yield t
+                for aa in a:
+                    t *= (aa + k - 1)
+                for bb in b:
+                    t /= (bb + k - 1)
+                t *= z
+                if t == 0:
+                    break
+                t /= k
+                k += 1
+
+        def deflated(self, ex, a, b, z):
+            """
+            Rewrite as a linear combination of functions of strictly lower
+            degree by eliminating all parameters a[i] and b[j] such that
+            a[i] = b[i] + m for nonnegative integer m.
+
+            EXAMPLES::
+
+                sage: x = hypergeometric([5], [4], 3)
+                sage: y = x.deflated()
+                sage: y
+                7/4*hypergeometric((), (), 3)
+                sage: x.n(); y.n()
+                35.1496896155784
+                35.1496896155784
+
+                sage: x = hypergeometric([6, 1], [3, 4, 5], 10)
+                sage: y = x.deflated()
+                sage: y
+                hypergeometric((1,), (4, 5), 10) + 1/2*hypergeometric((2,), (5, 6), 10) + 1/12*hypergeometric((3,), (6, 7), 10) + 1/252*hypergeometric((4,), (7, 8), 10)
+                sage: x.n(); y.n()
+                2.87893612686782
+                2.87893612686782
+
+                sage: x = hypergeometric([6, 7], [3, 4, 5], 10)
+                sage: y = x.deflated()
+                sage: y
+                hypergeometric((), (5,), 10) + 5*hypergeometric((), (6,), 10) + 19/3*hypergeometric((), (7,), 10) + 181/63*hypergeometric((), (8,), 10) + 265/504*hypergeometric((), (9,), 10) + 25/648*hypergeometric((), (10,), 10) + 25/27216*hypergeometric((), (11,), 10)
+                sage: x.n(); y.n()
+                63.0734110716969
+                63.0734110716969
+
+            """
+            return sum(map(prod, ex._deflated()))
+
+        def _deflated(self, ex, a, b, z):
+            new = ex.eliminate_parameters()
+            aa = new.operands()[0].operands()
+            bb = new.operands()[1].operands()
+            for i, aaa in enumerate(aa):
+                for j, bbb in enumerate(bb):
+                    m = aaa - bbb
+                    if m in ZZ and m > 0:
+                        aaaa = aa[:i] + aa[i + 1:]
+                        bbbb = bb[:j] + bb[j + 1:]
+                        terms = []
+                        for k in xrange(m + 1):
+                            # TODO: could rewrite prefactors as recurrence
+                            term = binomial(m, k)
+                            for c in aaaa:
+                                term *= rising_factorial(c, k)
+                            for c in bbbb:
+                                term /= rising_factorial(c, k)
+                            term *= z ** k
+                            term /= rising_factorial(aaa - m, k)
+                            F = hypergeometric([c + k for c in aaaa],
+                                               [c + k for c in bbbb], z)
+                            unique = []
+                            counts = []
+                            for c, f in F._deflated():
+                                if f in unique:
+                                    counts[unique.index(f)] += c
+                                else:
+                                    unique.append(f)
+                                    counts.append(c)
+                            Fterms = zip(counts, unique)
+                            terms += [(term*termG, G) for (termG, G) in
+                                      Fterms]
+                        return terms
+            return ((1, new),)
+
+hypergeometric = Hypergeometric()
+
+def closed_form(ex):
+    """
+    Try to evaluate self in closed form using elementary
+    (and other simple) functions.
+
+    It may be necessary to call :meth:`deflated` first to
+    find some closed forms.
+
+    EXAMPLES::
+
+        sage: from sage.functions.hypergeometric import closed_form
+        sage: var('a b c z')
+        (a, b, c, z)
+        sage: closed_form(hypergeometric([1], [], 1 + z))
+        -1/z
+        sage: closed_form(hypergeometric([], [], 1 + z))
+        e^(z + 1)
+        sage: closed_form(hypergeometric([], [1/2], 4))
+        cosh(4)
+        sage: closed_form(hypergeometric([], [3/2], 4))
+        1/4*sinh(4)
+        sage: closed_form(hypergeometric([], [5/2], 4))
+        -3/64*sinh(4) + 3/16*cosh(4)
+        sage: closed_form(hypergeometric([], [-3/2], 4))
+        -4*sinh(4) + 19/3*cosh(4)
+        sage: closed_form(hypergeometric([-3, 1], [var('a')], z))
+        -6*z^3/((a + 1)*(a + 2)*a) + 6*z^2/((a + 1)*a) - 3*z/a + 1
+        sage: closed_form(hypergeometric([-3,1/3],[-4],z))
+        7/162*z^3 + 1/9*z^2 + 1/4*z + 1
+        sage: closed_form(hypergeometric([],[],z))
+        e^z
+        sage: closed_form(hypergeometric([a],[],z))
+        (-z + 1)^(-a)
+        sage: closed_form(hypergeometric([1,1,2],[1,1],z))
+        (z - 1)^(-2)
+        sage: closed_form(hypergeometric([2,3],[1],x))
+        -1/(x - 1)^3 + 3*x/(x - 1)^4
+        sage: closed_form(hypergeometric([1/2],[3/2],-5))
+        1/10*sqrt(pi)*sqrt(5)*erf(sqrt(5))
+        sage: closed_form(hypergeometric([2],[5],3))
+        4
+        sage: closed_form(hypergeometric([2],[5],5))
+        48/625*e^5 + 612/625
+        sage: closed_form(hypergeometric([1/2,7/2],[3/2],z))
+        1/sqrt(-z + 1) + 2/3*z/(-z + 1)^(3/2) + 1/5*z^2/(-z + 1)^(5/2)
+        sage: closed_form(hypergeometric([1/2,1],[2],z))
+        -2*(sqrt(-z + 1) - 1)/z
+        sage: closed_form(hypergeometric([1,1],[2],z))
+        -log(-z + 1)/z
+        sage: closed_form(hypergeometric([1,1],[3],z))
+        -2*((z - 1)*log(-z + 1)/z - 1)/z
+
+    TODO: Try Python int
+    """
+    if ex.is_terminating():
+        return sum(ex.terms())
+
+    # necessary?
+    new = ex.eliminate_parameters()
+    #new = ex.deflated()
+    def _closed_form(ex):
+        a, b, z = ex.operands()
+        a, b = a.operands(), b.operands()
+        p, q = len(a), len(b)
+
+        if z == 0:
+            return Integer(1)
+        if p == q == 0:
+            return exp(z)
+        if p == 1 and q == 0:
+            return (1 - z) ** (-a[0])
+
+        if p == 0 and q == 1:
+            # TODO: make this require only linear time
+            def _0f1(b, z):
+                F12 = cosh(2 * sqrt(z))
+                F32 = sinh(2 * sqrt(z)) / (2 * sqrt(z))
+                if 2 * b == 1:
+                    return F12
+                if 2 * b == 3:
+                    return F32
+                if 2 * b > 3:
+                    return ((b - 2) * (b - 1) / z * (_0f1(b - 2, z) -
+                            _0f1(b - 1, z)))
+                if 2 * b < 1:
+                    return (_0f1(b + 1, z) + z / (b * (b + 1)) *
+                            _0f1(b + 2, z))
+                raise ValueError
+            # Can evaluate 0F1 in terms of elementary functions when
+            # the parameter is a half-integer
+            if 2*b[0] in ZZ and b[0] not in ZZ:
+                return _0f1(b[0], z)
+
+        # Confluent hypergeometric function
+        if p == 1 and q == 1:
+            aa, bb = a[0], b[0]
+            if aa * 2 == 1 and bb * 2 == 3:
+                t = sqrt(-z)
+                return sqrt(pi) / 2 * erf(t) / t
+            if a == 1 and b == 2:
+                return (exp(z) - 1) / z
+            n, m = aa, bb
+            if n in ZZ and m in ZZ and m > 0 and n > 0:
+                rf = rising_factorial
+                if m <= n:
+                    return (exp(z) * sum(rf(m - n, k) * (-z) ** k /
+                            factorial(k) / rf(m, k) for k in
+                            xrange(n - m + 1)))
+                else:
+                    T = sum(rf(n - m + 1, k) * z ** k /
+                            (factorial(k) * rf(2 - m, k)) for k in
+                            xrange(m - n))
+                    U = sum(rf(1 - n, k) * (-z) ** k /
+                            (factorial(k) * rf(2 - m, k)) for k in
+                            xrange(n))
+                    return (factorial(m - 2) * rf(1 - m, n) *
+                            z ** (1 - m) / factorial(n - 1) *
+                            (T - exp(z) * U))
+
+        if p == 2 and q == 1:
+            R12 = QQ('1/2')
+            R32 = QQ('3/2')
+            def _2f1(a, b, c, z):
+                """
+                Evaluation of 2F1(a, b, c, z), assuming a, b, c positive
+                integers or half-integers
+                """
+                if b == c:
+                    return (1 - z) ** (-a)
+                if a == c:
+                    return (1 - z) ** (-b)
+                if a == 0 or b == 0:
+                    return Integer(1)
+                if a > b:
+                    a, b = b, a
+                if b >= 2:
+                    F1 = _2f1(a, b - 1, c, z)
+                    F2 = _2f1(a, b - 2, c, z)
+                    q = (b - 1) * (z - 1)
+                    return (((c - 2 * b + 2 + (b - a - 1) * z) * F1 +
+                            (b - c - 1) * F2) / q)
+                if c > 2:
+                    # how to handle this case?
+                    if a - c + 1 == 0 or b - c + 1 == 0:
+                        raise NotImplementedError
+                    F1 = _2f1(a, b, c - 1, z)
+                    F2 = _2f1(a, b, c - 2, z)
+                    r1 = (c - 1) * (2 - c - (a + b - 2 * c + 3) * z)
+                    r2 = (c - 1) * (c - 2) * (1 - z)
+                    q = (a - c + 1) * (b - c + 1) * z
+                    return (r1 * F1 + r2 * F2) / q
+
+                if (a, b, c) == (R12, 1, 2):
+                    return (2 - 2 * sqrt(1 - z)) / z
+                if (a, b, c) == (1, 1, 2):
+                    return -log(1 - z) / z
+                if (a, b, c) == (1, R32, R12):
+                    return (1 + z) / (1 - z) ** 2
+                if (a, b, c) == (1, R32, 2):
+                    return 2 * (1 / sqrt(1 - z) - 1) / z
+                if (a, b, c) == (R32, 2, R12):
+                    return (1 + 3 * z) / (1 - z) ** 3
+                if (a, b, c) == (R32, 2, 1):
+                    return (2 + z) / (2 * (sqrt(1 - z) * (1 - z) ** 2))
+                if (a, b, c) == (2, 2, 1):
+                    return (1 + z) / (1 - z) ** 3
+                raise NotImplementedError
+            aa, bb = a
+            cc, = b
+            if z == 1:
+                return (gamma(cc) * gamma(cc - aa - bb) / gamma(cc - aa) /
+                        gamma(cc - bb))
+            if ((aa * 2) in ZZ and (bb * 2) in ZZ and (cc * 2) in ZZ and
+                aa > 0 and bb > 0 and cc > 0):
+                try:
+                    return _2f1(aa, bb, cc, z)
+                except NotImplementedError:
+                    pass
+    return sum([coeff * _closed_form(pfq) for coeff, pfq in new._deflated()])
+#    return new

sage/interfaces/maxima_lib.py

@@ -1158 +1158 @@
 
 
 ## Here we build dictionaries for operators needing special conversions.
-ratdisrep=EclObject("ratdisrep")
-mrat=EclObject("MRAT")
-mqapply=EclObject("MQAPPLY")
-max_li=EclObject("$LI")
-max_psi=EclObject("$PSI")
-max_array=EclObject("ARRAY")
-mdiff=EclObject("%DERIVATIVE")
-max_lambert_w=sage_op_dict[sage.functions.log.lambert_w]
+ratdisrep = EclObject("ratdisrep")
+mrat = EclObject("MRAT")
+mqapply = EclObject("MQAPPLY")
+max_li = EclObject("$LI")
+max_psi = EclObject("$PSI")
+max_hyper = EclObject("$%F")
+max_array = EclObject("ARRAY")
+mdiff = EclObject("%DERIVATIVE")
+max_lambert_w = sage_op_dict[sage.functions.log.lambert_w]
 
 def mrat_to_sage(expr):
     r"""
@@ -1221 +1222 @@
     """
     if caaadr(expr) == max_li:
         return sage.functions.log.polylog(max_to_sr(cadadr(expr)),
-                                           max_to_sr(caddr(expr)))
+                                          max_to_sr(caddr(expr)))
     if caaadr(expr) == max_psi:
         return sage.functions.other.psi(max_to_sr(cadadr(expr)),
-                                         max_to_sr(caddr(expr)))
+                                        max_to_sr(caddr(expr)))
+    if caaadr(expr) == max_hyper:
+        return sage.functions.hypergeometric.hypergeometric(mlist_to_sage(car(cdr(cdr(expr)))),
+                                                            mlist_to_sage(car(cdr(cdr(cdr(expr))))),
+                                                            max_to_sr(car(cdr(cdr(cdr(cdr(expr)))))))
     else:
         op=max_to_sr(cadr(expr))
         max_args=cddr(expr)
@@ -1260 +1265 @@
     
     - ``expr`` - ECL object; a Maxima MLIST expression (i.e., a list)
 
-    OUTPUT: a python list of converted expressions.
+    OUTPUT: a Python list of converted expressions.
     
     EXAMPLES::
     
@@ -1479 +1484 @@
             return EclObject(l)
         elif (op in special_sage_to_max):
             return EclObject(special_sage_to_max[op](*[sr_to_max(o) for o in expr.operands()]))
+        elif op == tuple:
+            return maxima(expr.operands()).ecl()
         elif not (op in sage_op_dict):
             # Maxima does some simplifications automatically by default
             # so calling maxima(expr) can change the structure of expr

sage/symbolic/expression.pyx

@@ -7469 +7469 @@
             factorial(n)
         """
         x = self
+        x = x.simplify_hypergeometric()
         x = x.simplify_factorial()
         x = x.simplify_trig()
         x = x.simplify_rational()
@@ -7479 +7480 @@
 
     full_simplify = simplify_full
 
+    def simplify_hypergeometric(self, algorithm='maxima'):
+        """
+        EXAMPLES::
+
+            sage: hypergeometric((5,4), (4,1,2,3),x).simplify_hypergeometric()          
+            1/144*x^2*hypergeometric((), (3, 4), x) + 1/3*x*hypergeometric((), (2, 3), x) + hypergeometric((), (1, 2), x)
+            sage: (2*hypergeometric((), (), x)).simplify_hypergeometric()
+            2*e^x
+        """
+        from sage.functions.hypergeometric import hypergeometric, closed_form
+        from sage.calculus.calculus import maxima
+        try:
+            op = self.operator()
+        except RuntimeError:
+            return self
+        ops = self.operands()
+        if op == hypergeometric:
+            if algorithm == 'maxima':
+                return self.parent()(maxima.hgfred(map(lambda o: o.simplify_hypergeometric(algorithm), ops[0].operands()), map(lambda o: o.simplify_hypergeometric(algorithm), ops[1].operands()), ops[2].simplify_hypergeometric(algorithm)))
+            elif algorithm == 'sage':
+                return closed_form(hypergeometric(map(lambda o: o.simplify_hypergeometric(algorithm), ops[0].operands()), map(lambda o: o.simplify_hypergeometric(algorithm), ops[1].operands()), ops[2].simplify_hypergeometric(algorithm)))
+            else:
+                return NotImplementedError('unknown algorithm')
+        if not op:
+            return self
+        return op(*map(lambda o: o.simplify_hypergeometric(algorithm), ops))
+
+    hypergeometric_simplify = simplify_hypergeometric
+
     def simplify_trig(self,expand=True):
         r""" 
         Optionally expands and then employs identities such as

sage/symbolic/expression_conversions.py

@@ -216 +216 @@
             return self.relation(ex, operator)
         elif isinstance(operator, FDerivativeOperator):
             return self.derivative(ex, operator)
+        elif operator == tuple:
+            return self.tuple(ex, operator)
         else:
             return self.composition(ex, operator)
 
@@ -459 +461 @@
         return "%s %s %s"%(self(ex.lhs()), self.relation_symbols[operator],
                            self(ex.rhs()))
 
+    def tuple(self, ex, operator):
+        x = map(self, ex.operands())
+        X = ','.join(x)
+        return '%s%s%s'%(self.interface._left_list_delim(), X, self.interface._right_list_delim())
+
     def derivative(self, ex, operator):
         """
         EXAMPLES::
@@ -1238 +1245 @@
             1.41421356237309...
         """
         f = operator
-        g = map(self, ex.operands())
+        try:
+            g = map(self, ex.operands())
+        except AttributeError, error:
+            from sage.functions.hypergeometric import hypergeometric
+            if f is hypergeometric:
+                ops = ex.operands()
+                return hypergeometric(*[map(lambda x: self.pyobject(None, x), ops[0]),
+                                        map(lambda x: self.pyobject(None, x), ops[1]),
+                                        self(ops[2])])
+            else:
+                raise AttributeError(error)
         try:
             return f(*g)
         except TypeError:
@@ -1247 +1264 @@
                 return abs(*g) # special case
             else:
                 return self.ff.fast_float_func(f, *g)
-            
+           
 def fast_float(ex, *vars):
     """
     Returns an object which provides fast floating point evaluation of