Ticket #2516: hyper3.patch

File hyper3.patch, 37.1 KB (added by eviatarbach, 9 years ago)
  • doc/en/reference/functions/index.rst

    # HG changeset patch
    # User fredrik@scv
    # Date 1280114065 -7200
    # Branch hyper
    # Node ID c5d9b1672bc133593098c60fad47903ff2bb46e1
    # Parent  42c2de7ba6a5461b04018606bdca93da60991117
    working on hypergeometric
    
    diff --git a/doc/en/reference/functions/index.rst b/doc/en/reference/functions/index.rst
    a b  
    1111   sage/functions/orthogonal_polys
    1212   sage/functions/other
    1313   sage/functions/special
     14   sage/functions/hypergeometric
    1415   sage/functions/exp_integral
    1516   sage/functions/wigner
    1617   sage/functions/generalized
  • sage/calculus/calculus.py

    diff --git a/sage/calculus/calculus.py b/sage/calculus/calculus.py
    a b  
    16851685
    16861686maxima_qp = re.compile("\?\%[a-z|A-Z|0-9|_]*")  # e.g., ?%jacobi_cd
    16871687
    1688 maxima_var = re.compile("\%[a-z|A-Z|0-9|_]*")  # e.g., ?%jacobi_cd
     1688maxima_var = re.compile("\%[a-z|A-Z|0-9|_]*")  # e.g., %jacobi_cd
    16891689
    16901690sci_not = re.compile("(-?(?:0|[1-9]\d*))(\.\d+)?([eE][-+]\d+)")
    16911691
     
    16931693
    16941694maxima_polygamma = re.compile("psi\[(\d*)\]\(")  # matches psi[n]( where n is a number
    16951695
     1696maxima_hyper = re.compile("\%f\[\d+,\d+\]")  # matches %f[m,n]
     1697
    16961698def symbolic_expression_from_maxima_string(x, equals_sub=False, maxima=maxima):
    16971699    """
    16981700    Given a string representation of a Maxima expression, parse it and
     
    17781780                pass
    17791781            else:
    17801782                syms[X[2:]] = function_factory(X[2:])
    1781         s = s.replace("?%","")
     1783        s = s.replace("?%", "")
    17821784
    1783     s = polylog_ex.sub('polylog(\\1,',s)
     1785    s = maxima_hyper.sub('hypergeometric', s)
     1786
     1787    s = polylog_ex.sub('polylog(\\1,', s)
    17841788    s = multiple_replace(symtable, s)
    17851789    s = s.replace("%","")
    17861790
    17871791    s = s.replace("#","!=") # a lot of this code should be refactored somewhere...
    17881792
    1789     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
     1793    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
    17901794
    17911795    if equals_sub:
    17921796        s = s.replace('=','==')
  • sage/functions/all.py

    diff --git a/sage/functions/all.py b/sage/functions/all.py
    a b  
    6767                          sin_integral, cos_integral, Si, Ci,
    6868                          sinh_integral, cosh_integral, Shi, Chi,
    6969                          exponential_integral_1, Ei)
     70
     71from hypergeometric import hypergeometric
  • new file sage/functions/hypergeometric.py

    diff --git a/sage/functions/hypergeometric.py b/sage/functions/hypergeometric.py
    new file mode 100644
    - +  
     1r"""
     2Hypergeometric Functions
     3
     4This module implements manipulation of infinite hypergeometric series
     5represented in standard parametric form (as $\,_pFq$ functions).
     6
     7AUTHORS:
     8
     9- Fredrik Johansson (2010): initial version
     10
     11- Eviatar Bach (2013): major changes
     12
     13EXAMPLES:
     14
     15Examples from :trac:`9908`::
     16
     17    sage: maxima('integrate(bessel_j(2, x), x)').sage()
     18    1/24*x^3*hypergeometric((3/2,), (5/2, 3), -1/4*x^2)
     19    sage: sum(((2 * I)^x / (x^3 + 1) * (1/4)^x), x, 0, oo)
     20    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)
     21    sage: sum((-1)^x / ((2 * x + 1) * factorial(2 * x + 1)), x, 0, oo)
     22    hypergeometric((1/2,), (3/2, 3/2), -1/4)
     23
     24Simplification::
     25
     26    sage: hypergeometric((hypergeometric((), (), x),), (),
     27    ....:                hypergeometric((), (), x)).simplify_hypergeometric()
     28    1/((-e^x + 1)^e^x)
     29    sage: hypergeometric([-2], [], x).simplify_hypergeometric()
     30    x^2 - 2*x + 1
     31    sage: hypergeometric([], [], x).simplify_hypergeometric()
     32    e^x
     33    sage: hypergeometric((hypergeometric((), (), x),), (),
     34    ....:                hypergeometric((), (), x)).simplify_hypergeometric(algorithm='sage')
     35    (-e^x + 1)^(-e^x)
     36
     37Equality testing::
     38
     39    sage: bool(hypergeometric([], [], x).derivative(x) ==
     40    ....:      hypergeometric([], [], x))  # diff(e^x, x) == e^x
     41    True
     42    sage: bool(hypergeometric([], [], x) == hypergeometric([], [1], x))
     43    False
     44
     45Computing terms and series::
     46
     47    sage: z = var('z')
     48    sage: hypergeometric([], [], z).series(z, 0)
     49    Order(1)
     50    sage: hypergeometric([], [], z).series(z, 1)
     51    1 + Order(z)
     52    sage: hypergeometric([], [], z).series(z, 2)
     53    1 + 1*z + Order(z^2)
     54    sage: hypergeometric([], [], z).series(z, 3)
     55    1 + 1*z + 1/2*z^2 + Order(z^3)
     56
     57    sage: hypergeometric([-2], [], z).series(z, 3)
     58    1 + (-2)*z + 1*z^2
     59    sage: hypergeometric([-2], [], z).series(z, 6)
     60    1 + (-2)*z + 1*z^2
     61    sage: hypergeometric([-2], [], z).series(z, 6).is_terminating_series()
     62    True
     63    sage: hypergeometric([-2], [], z).series(z, 2)                       
     64    1 + (-2)*z + Order(z^2)
     65    sage: hypergeometric([-2], [], z).series(z, 2).is_terminating_series()
     66    False
     67
     68    sage: hypergeometric([1], [], z).series(z, 6)
     69    1 + 1*z + 1*z^2 + 1*z^3 + 1*z^4 + 1*z^5 + Order(z^6)
     70    sage: hypergeometric([], [1/2], -z^2/4).series(z, 11)
     71    1 + (-1/2)*z^2 + 1/24*z^4 + (-1/720)*z^6 + 1/40320*z^8 + (-1/3628800)*z^10 + Order(z^11)
     72
     73    sage: hypergeometric([1],[5],x).series(x,5)
     74    1 + 1/5*x + 1/30*x^2 + 1/210*x^3 + 1/1680*x^4 + Order(x^5)
     75
     76    sage: sum(hypergeometric([1, 2], [3], 1/3).terms(6)).n()
     77    1.29788359788360
     78    sage: hypergeometric([1, 2], [3], 1/3).n()
     79    1.29837194594696
     80    sage: hypergeometric([], [], x).series(x, 20)(x=1).n() == e.n()
     81    True
     82
     83Plotting::
     84
     85    sage: plot(hypergeometric([1, 1], [3, 3, 3], x), x, -30, 30)
     86
     87Numeric evaluation::
     88
     89    sage: hypergeometric([1], [], 1/10).n()  # geometric series
     90    1.11111111111111
     91    sage: hypergeometric([], [], 1).n()  # e
     92    2.71828182845905
     93    sage: hypergeometric([], [], 3., hold=True)
     94    hypergeometric((), (), 3.00000000000000)
     95    sage: hypergeometric([1,2,3],[4,5,6],1/2).n()
     96    1.02573619590134
     97    sage: hypergeometric([1,2,3],[4,5,6],1/2).n(digits=30)
     98    1.02573619590133865036584139535
     99    sage: hypergeometric([5-3*I],[3/2,2+I,sqrt(2)],4+I).n() #TOLERANCE?
     100    5.52605111678805 - 7.86331357527544*I
     101
     102SymPy conversion::
     103
     104    sage: hypergeometric((5,4), (4,4), 3)._sympy_()
     105    hyper((5, 4), (4, 4), 3)
     106
     107"""
     108#*****************************************************************************
     109#       Copyright (C) 2010 Fredrik Johansson <fredrik.johansson@gmail.com>
     110#       Copyright (C) 2013 Eviatar Bach <eviatarbach@gmail.com>
     111#
     112#  Distributed under the terms of the GNU General Public License (GPL)
     113#  as published by the Free Software Foundation; either version 2 of
     114#  the License, or (at your option) any later version.
     115#                  http://www.gnu.org/licenses/
     116#*****************************************************************************
     117
     118from sage.rings.all import (Integer, ZZ, QQ, Infinity, binomial,
     119                            rising_factorial, factorial)
     120from sage.functions.other import sqrt, gamma, real_part
     121from sage.functions.log import exp, log
     122from sage.functions.trig import cos, sin
     123from sage.functions.hyperbolic import cosh, sinh
     124from sage.functions.other import erf
     125from sage.symbolic.constants import pi
     126from sage.symbolic.all import I
     127from sage.symbolic.function import BuiltinFunction, is_inexact
     128from sage.symbolic.ring import SR
     129from sage.structure.element import get_coercion_model
     130from sage.misc.latex import latex
     131from sage.misc.flatten import flatten
     132from sage.misc.misc_c import prod
     133from mpmath import hyper
     134from sage.libs.mpmath import utils as mpmath_utils
     135from sage.symbolic.expression import Expression
     136from sage.interfaces.maxima import maxima
     137from sage.calculus.functional import derivative
     138
     139def rational_param_as_tuple(x):
     140    """
     141    Utility function for converting rational PFQ parameters to
     142    tuples (which mpmath handles more efficiently).
     143
     144        sage: from sage.functions.hypergeometric import rational_param_as_tuple
     145        sage: rational_param_as_tuple(1/2)
     146        (1, 2)
     147        sage: rational_param_as_tuple(3)
     148        3
     149        sage: rational_param_as_tuple(pi)
     150        pi
     151
     152    """
     153    try:
     154        x = x.pyobject()
     155    except AttributeError:
     156        pass
     157    try:
     158        if x.parent() is QQ:
     159            p = int(x.numer())
     160            q = int(x.denom())
     161            return p, q
     162    except AttributeError:
     163        pass
     164    return x
     165
     166
     167class Hypergeometric(BuiltinFunction):
     168    r"""
     169    Represents a (formal) generalized infinite hypergeometric series.
     170    """
     171
     172    def __init__(self):
     173        """
     174        EXAMPLES::
     175
     176            sage: hypergeometric([], [], 1)
     177            hypergeometric((), (), 1)
     178            sage: hypergeometric([], [1], 1)
     179            hypergeometric((), (1,), 1)
     180            sage: hypergeometric([2, 3], [1], 1)
     181            hypergeometric((2, 3), (1,), 1)
     182
     183        """
     184        BuiltinFunction.__init__(self, 'hypergeometric', nargs=3,
     185                                 conversions={'mathematica':'HypergeometricPFQ',
     186                                              'maxima':'hypergeometric',
     187                                              'sympy':'hyper'})
     188
     189    def __call__(self, a, b, z, **kwargs):
     190        return BuiltinFunction.__call__(self,
     191                                        SR._force_pyobject(a),
     192                                        SR._force_pyobject(b),
     193                                        z, **kwargs)
     194
     195    def _print_latex_(self, a, b, z):
     196        r"""
     197        TESTS::
     198
     199            sage: latex(hypergeometric([1, 1], [2], -1))
     200            \,_2F_1\left(\begin{matrix} 1,1 \\ 2 \end{matrix} ; -1 \right)
     201
     202        """
     203        aa = ",".join(latex(c) for c in a)
     204        bb = ",".join(latex(c) for c in b)
     205        z = latex(z)
     206        return (r"\,_{}F_{}\left(\begin{{matrix}} {} \\ {} \end{{matrix}} ; "
     207                r"{} \right)").format(len(a), len(b), aa, bb, z)
     208
     209    def _eval_(self, a, b, z, **kwargs):
     210        coercion_model = get_coercion_model()
     211        co = reduce(lambda x, y: coercion_model.canonical_coercion(x, y)[0],
     212                    a + b + (z,))
     213        if is_inexact(co) and not isinstance(co, Expression):
     214            from sage.structure.coerce import parent
     215            return self._evalf_(a, b, z, parent=parent(co))
     216        if z == 0:
     217            return Integer(1)
     218        return
     219
     220    def _evalf_(self, a, b, z, parent):
     221        """
     222        TESTS::
     223
     224            sage: hypergeometric([1, 1], [2], -1).n()
     225            0.693147180559945
     226            sage: hypergeometric([], [], RealField(100)(1))
     227            2.7182818284590452353602874714
     228
     229        """
     230        aa = [rational_param_as_tuple(c) for c in a]
     231        bb = [rational_param_as_tuple(c) for c in b]
     232        return mpmath_utils.call(hyper, aa, bb, z, parent=parent)
     233
     234    def _tderivative_(self, a, b, z, *args, **kwargs):
     235        """
     236        EXAMPLES::
     237
     238            sage: hypergeometric([1/3, 2/3], [5], x^2).diff(x)
     239            4/45*x*hypergeometric((4/3, 5/3), (6,), x^2)
     240            sage: hypergeometric([1, 2], [x], 2).diff(x)
     241            Traceback (most recent call last):
     242                ...
     243            NotImplementedError: derivative of hypergeometric function with respect to parameters. Try calling `.simplify_hypergeometric()` first.
     244        """
     245        diff_param = kwargs['diff_param']
     246        if diff_param in hypergeometric(a, b, 1).variables():
     247            raise NotImplementedError("derivative of hypergeometric function "
     248                                      "with respect to parameters. Try calling"                                      " `.simplify_hypergeometric()` first.")
     249        t = (reduce(lambda x, y: x * y, a, 1) *
     250             reduce(lambda x, y: x / y, b, Integer(1)))
     251        return (t * derivative(z, diff_param) *
     252                hypergeometric([c + 1 for c in a], [c + 1 for c in b], z))
     253
     254    '''def _maxima_init_evaled_(self, a, b, z):
     255        return "hypergeometric({},{},{})".format(repr(maxima(a)),
     256                                                 repr(maxima(b)),
     257                                                 repr(maxima(z)))
     258
     259    def _maxima_lib_(self, a, b, z):
     260        return "hypergeometric({},{},{})".format(repr(maxima(a.operands())),
     261                                                 repr(maxima(b.operands())),
     262                                                 repr(maxima(z)))'''
     263
     264    class EvaluationMethods:
     265        '''def operands(self, ex, a, b, z):
     266            return [(a), list(b), z]'''
     267
     268        def sorted_parameters(self, ex, a, b, z):
     269            """
     270            Return with parameters sorted in a canonical order
     271
     272            EXAMPLES::
     273
     274                sage: hypergeometric([2,1,3], [5,4], 1/2).sorted_parameters()
     275                hypergeometric((1, 2, 3), (4, 5), 1/2)
     276
     277            """
     278            return hypergeometric(sorted(a), sorted(b), z)
     279
     280        def eliminate_parameters(self, ex, a, b, z):
     281            """
     282            Eliminate repeated parameters by pairwise cancellation of identical
     283            terms in ``a`` and ``b``
     284
     285            EXAMPLES::
     286
     287                sage: hypergeometric([1, 1, 2, 5], [5, 1, 4],
     288                ....:                1/2).eliminate_parameters()
     289                hypergeometric((1, 2), (4,), 1/2)
     290                sage: hypergeometric([x], [x], x).eliminate_parameters()
     291                hypergeometric((), (), x)
     292                sage: hypergeometric((5, 4), (4, 4), 3).eliminate_parameters()
     293                hypergeometric((5,), (4,), 3)
     294
     295            """
     296            aa = list(a)
     297            bb = list(b)
     298            p = pp = len(aa)
     299            q = qq = len(bb)
     300            i = 0
     301            while i < qq and aa:
     302                bbb = bb[i]
     303                if bbb in aa:
     304                    aa.remove(bbb)
     305                    bb.remove(bbb)
     306                    pp -= 1
     307                    qq -= 1
     308                else:
     309                    i += 1
     310            if (pp, qq) != (p, q):
     311                return hypergeometric(aa, bb, z)
     312            return ex
     313
     314        def is_termwise_finite(self, ex, a, b, z):
     315            """
     316            Determine whether all terms of self are finite. Any infinite
     317            terms or ambiguous terms beyond the first zero, if one exists,
     318            are ignored.
     319
     320            Ambiguous cases (where a term is the product of both zero
     321            and an infinity) are not considered finite.
     322
     323
     324                sage: hypergeometric([2], [3, 4], 5).is_termwise_finite()
     325                True
     326                sage: hypergeometric([2], [-3, 4], 5).is_termwise_finite()
     327                False
     328                sage: hypergeometric([-2], [-3, 4], 5).is_termwise_finite()
     329                True
     330                sage: hypergeometric([-3], [-3, 4], 5).is_termwise_finite()  # ambiguous
     331                False
     332
     333                sage: hypergeometric([0], [-1], 5).is_termwise_finite()
     334                True
     335                sage: hypergeometric([0], [0], 5).is_termwise_finite()  # ambiguous
     336                False
     337                sage: hypergeometric([1], [2], Infinity).is_termwise_finite()
     338                False
     339                sage: hypergeometric([0], [0], Infinity).is_termwise_finite()  # ambiguous
     340                False
     341                sage: hypergeometric([0], [], Infinity).is_termwise_finite()  # ambiguous
     342                False
     343
     344            """
     345            if z == 0:
     346                return 0 not in b
     347            if abs(z) == Infinity:
     348                return False
     349            if abs(z) == Infinity:
     350                return False
     351            for bb in b:
     352                if bb in ZZ and bb <= 0:
     353                    if any((aa in ZZ) and (bb < aa <= 0) for aa in a):
     354                        continue
     355                    return False
     356            return True
     357
     358        def is_terminating(self, ex, a, b, z):
     359            """
     360            Determine whether the series represented by self terminates
     361            after a finite number of terms, i.e. whether any of the
     362            numerator parameters are nonnegative integers (with no
     363            preceding nonnegative denominator parameters), or z = 0.
     364
     365            If terminating, the series represents a polynomial of z
     366
     367            EXAMPLES::
     368
     369                sage: hypergeometric([1, 2], [3, 4], x).is_terminating()
     370                False
     371                sage: hypergeometric([1, -2], [3, 4], x).is_terminating()
     372                True
     373                sage: hypergeometric([1, -2], [], x).is_terminating()
     374                True
     375
     376            """
     377            if z == 0:
     378                return True
     379            for aa in a:
     380                if (aa in ZZ) and (aa <= 0):
     381                    return ex.is_termwise_finite()
     382            return False
     383
     384        def is_absolutely_convergent(self, ex, a, b, z):
     385            """
     386            Determine whether self converges absolutely as an infinite series.
     387            False is returned if not all terms are finite.
     388
     389            EXAMPLES:
     390
     391            Degree giving infinite radius of convergence::
     392
     393                sage: hypergeometric([2,3],[4,5],6).is_absolutely_convergent()
     394                True
     395                sage: hypergeometric([2,3],[-4,5],6).is_absolutely_convergent()  # undefined
     396                False
     397                sage: hypergeometric([2,3],[-4,5],Infinity).is_absolutely_convergent()  # undefined
     398                False
     399
     400            Ordinary geometric series (unit radius of convergence)::
     401
     402                sage: hypergeometric([1],[],1/2).is_absolutely_convergent()
     403                True
     404                sage: hypergeometric([1],[],2).is_absolutely_convergent()
     405                False
     406                sage: hypergeometric([1],[],1).is_absolutely_convergent()
     407                False
     408                sage: hypergeometric([1],[],-1).is_absolutely_convergent()
     409                False
     410                sage: hypergeometric([1],[],-1).n()    # Sum still exists
     411                0.500000000000000
     412
     413            Degree $p = q+1$ (unit radius of convergence)::
     414
     415                sage: hypergeometric([2,3],[4],6).is_absolutely_convergent()
     416                False
     417                sage: hypergeometric([2,3],[4],1).is_absolutely_convergent()
     418                False
     419                sage: hypergeometric([2,3],[5],1).is_absolutely_convergent()
     420                False
     421                sage: hypergeometric([2,3],[6],1).is_absolutely_convergent()
     422                True
     423                sage: hypergeometric([-2,3],[4],5).is_absolutely_convergent()
     424                True
     425                sage: hypergeometric([2,-3],[4],5).is_absolutely_convergent()
     426                True
     427                sage: hypergeometric([2,-3],[-4],5).is_absolutely_convergent()
     428                True
     429                sage: hypergeometric([2,-3],[-1],5).is_absolutely_convergent()
     430                False
     431
     432            Degree giving zero radius of convergence::
     433
     434                sage: hypergeometric([1,2,3],[4],2).is_absolutely_convergent()
     435                False
     436                sage: hypergeometric([1,2,3],[4],1/2).is_absolutely_convergent()
     437                False
     438                sage: hypergeometric([1,2,-3],[4],1/2).is_absolutely_convergent()  # polynomial
     439                True
     440
     441            """
     442            p, q = len(a), len(b)
     443            if not ex.is_termwise_finite():
     444                return False
     445            if p <= q:
     446                return True
     447            if ex.is_terminating():
     448                return True
     449            if p == q + 1:
     450                if abs(z) < 1:
     451                    return True
     452                if abs(z) == 1:
     453                    if real_part(sum(b) - sum(a)) > 0:
     454                        return True
     455            return False
     456
     457        def terms(self, ex, a, b, z, n=None):
     458            """
     459            Generate the terms of self (optionally only n terms).
     460
     461            EXAMPLES::
     462
     463                sage: list(hypergeometric([-2, 1], [3, 4], x).terms())
     464                [1, -1/6*x, 1/120*x^2]
     465                sage: list(hypergeometric([-2, 1], [3, 4], x).terms(2))
     466                [1, -1/6*x]
     467                sage: list(hypergeometric([-2, 1], [3, 4], x).terms(0))
     468                []
     469
     470            """
     471            if n is None:
     472                n = Infinity
     473            t = Integer(1)
     474            k = 1
     475            while k <= n:
     476                yield t
     477                for aa in a:
     478                    t *= (aa + k - 1)
     479                for bb in b:
     480                    t /= (bb + k - 1)
     481                t *= z
     482                if t == 0:
     483                    break
     484                t /= k
     485                k += 1
     486
     487        def deflated(self, ex, a, b, z):
     488            """
     489            Rewrite as a linear combination of functions of strictly lower
     490            degree by eliminating all parameters a[i] and b[j] such that
     491            a[i] = b[i] + m for nonnegative integer m.
     492
     493            EXAMPLES::
     494
     495                sage: x = hypergeometric([5], [4], 3)
     496                sage: y = x.deflated()
     497                sage: y
     498                7/4*hypergeometric((), (), 3)
     499                sage: x.n(); y.n()
     500                35.1496896155784
     501                35.1496896155784
     502
     503                sage: x = hypergeometric([6, 1], [3, 4, 5], 10)
     504                sage: y = x.deflated()
     505                sage: y
     506                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)
     507                sage: x.n(); y.n()
     508                2.87893612686782
     509                2.87893612686782
     510
     511                sage: x = hypergeometric([6, 7], [3, 4, 5], 10)
     512                sage: y = x.deflated()
     513                sage: y
     514                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)
     515                sage: x.n(); y.n()
     516                63.0734110716969
     517                63.0734110716969
     518
     519            """
     520            return sum(map(prod, ex._deflated()))
     521
     522        def _deflated(self, ex, a, b, z):
     523            new = ex.eliminate_parameters()
     524            aa = new.operands()[0].operands()
     525            bb = new.operands()[1].operands()
     526            for i, aaa in enumerate(aa):
     527                for j, bbb in enumerate(bb):
     528                    m = aaa - bbb
     529                    if m in ZZ and m > 0:
     530                        aaaa = aa[:i] + aa[i + 1:]
     531                        bbbb = bb[:j] + bb[j + 1:]
     532                        terms = []
     533                        for k in xrange(m + 1):
     534                            # TODO: could rewrite prefactors as recurrence
     535                            term = binomial(m, k)
     536                            for c in aaaa:
     537                                term *= rising_factorial(c, k)
     538                            for c in bbbb:
     539                                term /= rising_factorial(c, k)
     540                            term *= z ** k
     541                            term /= rising_factorial(aaa - m, k)
     542                            F = hypergeometric([c + k for c in aaaa],
     543                                               [c + k for c in bbbb], z)
     544                            unique = []
     545                            counts = []
     546                            for c, f in F._deflated():
     547                                if f in unique:
     548                                    counts[unique.index(f)] += c
     549                                else:
     550                                    unique.append(f)
     551                                    counts.append(c)
     552                            Fterms = zip(counts, unique)
     553                            terms += [(term*termG, G) for (termG, G) in
     554                                      Fterms]
     555                        return terms
     556            return ((1, new),)
     557
     558hypergeometric = Hypergeometric()
     559
     560def closed_form(ex):
     561    """
     562    Try to evaluate self in closed form using elementary
     563    (and other simple) functions.
     564
     565    It may be necessary to call :meth:`deflated` first to
     566    find some closed forms.
     567
     568    EXAMPLES::
     569
     570        sage: from sage.functions.hypergeometric import closed_form
     571        sage: var('a b c z')
     572        (a, b, c, z)
     573        sage: closed_form(hypergeometric([1], [], 1 + z))
     574        -1/z
     575        sage: closed_form(hypergeometric([], [], 1 + z))
     576        e^(z + 1)
     577        sage: closed_form(hypergeometric([], [1/2], 4))
     578        cosh(4)
     579        sage: closed_form(hypergeometric([], [3/2], 4))
     580        1/4*sinh(4)
     581        sage: closed_form(hypergeometric([], [5/2], 4))
     582        -3/64*sinh(4) + 3/16*cosh(4)
     583        sage: closed_form(hypergeometric([], [-3/2], 4))
     584        -4*sinh(4) + 19/3*cosh(4)
     585        sage: closed_form(hypergeometric([-3, 1], [var('a')], z))
     586        -6*z^3/((a + 1)*(a + 2)*a) + 6*z^2/((a + 1)*a) - 3*z/a + 1
     587        sage: closed_form(hypergeometric([-3,1/3],[-4],z))
     588        7/162*z^3 + 1/9*z^2 + 1/4*z + 1
     589        sage: closed_form(hypergeometric([],[],z))
     590        e^z
     591        sage: closed_form(hypergeometric([a],[],z))
     592        (-z + 1)^(-a)
     593        sage: closed_form(hypergeometric([1,1,2],[1,1],z))
     594        (z - 1)^(-2)
     595        sage: closed_form(hypergeometric([2,3],[1],x))
     596        -1/(x - 1)^3 + 3*x/(x - 1)^4
     597        sage: closed_form(hypergeometric([1/2],[3/2],-5))
     598        1/10*sqrt(pi)*sqrt(5)*erf(sqrt(5))
     599        sage: closed_form(hypergeometric([2],[5],3))
     600        4
     601        sage: closed_form(hypergeometric([2],[5],5))
     602        48/625*e^5 + 612/625
     603        sage: closed_form(hypergeometric([1/2,7/2],[3/2],z))
     604        1/sqrt(-z + 1) + 2/3*z/(-z + 1)^(3/2) + 1/5*z^2/(-z + 1)^(5/2)
     605        sage: closed_form(hypergeometric([1/2,1],[2],z))
     606        -2*(sqrt(-z + 1) - 1)/z
     607        sage: closed_form(hypergeometric([1,1],[2],z))
     608        -log(-z + 1)/z
     609        sage: closed_form(hypergeometric([1,1],[3],z))
     610        -2*((z - 1)*log(-z + 1)/z - 1)/z
     611
     612    TODO: Try Python int
     613    """
     614    if ex.is_terminating():
     615        return sum(ex.terms())
     616
     617    # necessary?
     618    new = ex.eliminate_parameters()
     619    #new = ex.deflated()
     620    def _closed_form(ex):
     621        a, b, z = ex.operands()
     622        a, b = a.operands(), b.operands()
     623        p, q = len(a), len(b)
     624
     625        if z == 0:
     626            return Integer(1)
     627        if p == q == 0:
     628            return exp(z)
     629        if p == 1 and q == 0:
     630            return (1 - z) ** (-a[0])
     631
     632        if p == 0 and q == 1:
     633            # TODO: make this require only linear time
     634            def _0f1(b, z):
     635                F12 = cosh(2 * sqrt(z))
     636                F32 = sinh(2 * sqrt(z)) / (2 * sqrt(z))
     637                if 2 * b == 1:
     638                    return F12
     639                if 2 * b == 3:
     640                    return F32
     641                if 2 * b > 3:
     642                    return ((b - 2) * (b - 1) / z * (_0f1(b - 2, z) -
     643                            _0f1(b - 1, z)))
     644                if 2 * b < 1:
     645                    return (_0f1(b + 1, z) + z / (b * (b + 1)) *
     646                            _0f1(b + 2, z))
     647                raise ValueError
     648            # Can evaluate 0F1 in terms of elementary functions when
     649            # the parameter is a half-integer
     650            if 2*b[0] in ZZ and b[0] not in ZZ:
     651                return _0f1(b[0], z)
     652
     653        # Confluent hypergeometric function
     654        if p == 1 and q == 1:
     655            aa, bb = a[0], b[0]
     656            if aa * 2 == 1 and bb * 2 == 3:
     657                t = sqrt(-z)
     658                return sqrt(pi) / 2 * erf(t) / t
     659            if a == 1 and b == 2:
     660                return (exp(z) - 1) / z
     661            n, m = aa, bb
     662            if n in ZZ and m in ZZ and m > 0 and n > 0:
     663                rf = rising_factorial
     664                if m <= n:
     665                    return (exp(z) * sum(rf(m - n, k) * (-z) ** k /
     666                            factorial(k) / rf(m, k) for k in
     667                            xrange(n - m + 1)))
     668                else:
     669                    T = sum(rf(n - m + 1, k) * z ** k /
     670                            (factorial(k) * rf(2 - m, k)) for k in
     671                            xrange(m - n))
     672                    U = sum(rf(1 - n, k) * (-z) ** k /
     673                            (factorial(k) * rf(2 - m, k)) for k in
     674                            xrange(n))
     675                    return (factorial(m - 2) * rf(1 - m, n) *
     676                            z ** (1 - m) / factorial(n - 1) *
     677                            (T - exp(z) * U))
     678
     679        if p == 2 and q == 1:
     680            R12 = QQ('1/2')
     681            R32 = QQ('3/2')
     682            def _2f1(a, b, c, z):
     683                """
     684                Evaluation of 2F1(a, b, c, z), assuming a, b, c positive
     685                integers or half-integers
     686                """
     687                if b == c:
     688                    return (1 - z) ** (-a)
     689                if a == c:
     690                    return (1 - z) ** (-b)
     691                if a == 0 or b == 0:
     692                    return Integer(1)
     693                if a > b:
     694                    a, b = b, a
     695                if b >= 2:
     696                    F1 = _2f1(a, b - 1, c, z)
     697                    F2 = _2f1(a, b - 2, c, z)
     698                    q = (b - 1) * (z - 1)
     699                    return (((c - 2 * b + 2 + (b - a - 1) * z) * F1 +
     700                            (b - c - 1) * F2) / q)
     701                if c > 2:
     702                    # how to handle this case?
     703                    if a - c + 1 == 0 or b - c + 1 == 0:
     704                        raise NotImplementedError
     705                    F1 = _2f1(a, b, c - 1, z)
     706                    F2 = _2f1(a, b, c - 2, z)
     707                    r1 = (c - 1) * (2 - c - (a + b - 2 * c + 3) * z)
     708                    r2 = (c - 1) * (c - 2) * (1 - z)
     709                    q = (a - c + 1) * (b - c + 1) * z
     710                    return (r1 * F1 + r2 * F2) / q
     711
     712                if (a, b, c) == (R12, 1, 2):
     713                    return (2 - 2 * sqrt(1 - z)) / z
     714                if (a, b, c) == (1, 1, 2):
     715                    return -log(1 - z) / z
     716                if (a, b, c) == (1, R32, R12):
     717                    return (1 + z) / (1 - z) ** 2
     718                if (a, b, c) == (1, R32, 2):
     719                    return 2 * (1 / sqrt(1 - z) - 1) / z
     720                if (a, b, c) == (R32, 2, R12):
     721                    return (1 + 3 * z) / (1 - z) ** 3
     722                if (a, b, c) == (R32, 2, 1):
     723                    return (2 + z) / (2 * (sqrt(1 - z) * (1 - z) ** 2))
     724                if (a, b, c) == (2, 2, 1):
     725                    return (1 + z) / (1 - z) ** 3
     726                raise NotImplementedError
     727            aa, bb = a
     728            cc, = b
     729            if z == 1:
     730                return (gamma(cc) * gamma(cc - aa - bb) / gamma(cc - aa) /
     731                        gamma(cc - bb))
     732            if ((aa * 2) in ZZ and (bb * 2) in ZZ and (cc * 2) in ZZ and
     733                aa > 0 and bb > 0 and cc > 0):
     734                try:
     735                    return _2f1(aa, bb, cc, z)
     736                except NotImplementedError:
     737                    pass
     738    return sum([coeff * _closed_form(pfq) for coeff, pfq in new._deflated()])
     739#    return new
  • sage/interfaces/maxima_lib.py

    diff --git a/sage/interfaces/maxima_lib.py b/sage/interfaces/maxima_lib.py
    a b  
    11581158
    11591159
    11601160## Here we build dictionaries for operators needing special conversions.
    1161 ratdisrep=EclObject("ratdisrep")
    1162 mrat=EclObject("MRAT")
    1163 mqapply=EclObject("MQAPPLY")
    1164 max_li=EclObject("$LI")
    1165 max_psi=EclObject("$PSI")
    1166 max_array=EclObject("ARRAY")
    1167 mdiff=EclObject("%DERIVATIVE")
    1168 max_lambert_w=sage_op_dict[sage.functions.log.lambert_w]
     1161ratdisrep = EclObject("ratdisrep")
     1162mrat = EclObject("MRAT")
     1163mqapply = EclObject("MQAPPLY")
     1164max_li = EclObject("$LI")
     1165max_psi = EclObject("$PSI")
     1166max_hyper = EclObject("$%F")
     1167max_array = EclObject("ARRAY")
     1168mdiff = EclObject("%DERIVATIVE")
     1169max_lambert_w = sage_op_dict[sage.functions.log.lambert_w]
    11691170
    11701171def mrat_to_sage(expr):
    11711172    r"""
     
    12211222    """
    12221223    if caaadr(expr) == max_li:
    12231224        return sage.functions.log.polylog(max_to_sr(cadadr(expr)),
    1224                                            max_to_sr(caddr(expr)))
     1225                                          max_to_sr(caddr(expr)))
    12251226    if caaadr(expr) == max_psi:
    12261227        return sage.functions.other.psi(max_to_sr(cadadr(expr)),
    1227                                          max_to_sr(caddr(expr)))
     1228                                        max_to_sr(caddr(expr)))
     1229    if caaadr(expr) == max_hyper:
     1230        return sage.functions.hypergeometric.hypergeometric(mlist_to_sage(car(cdr(cdr(expr)))),
     1231                                                            mlist_to_sage(car(cdr(cdr(cdr(expr))))),
     1232                                                            max_to_sr(car(cdr(cdr(cdr(cdr(expr)))))))
    12281233    else:
    12291234        op=max_to_sr(cadr(expr))
    12301235        max_args=cddr(expr)
     
    12601265   
    12611266    - ``expr`` - ECL object; a Maxima MLIST expression (i.e., a list)
    12621267
    1263     OUTPUT: a python list of converted expressions.
     1268    OUTPUT: a Python list of converted expressions.
    12641269   
    12651270    EXAMPLES::
    12661271   
     
    14791484            return EclObject(l)
    14801485        elif (op in special_sage_to_max):
    14811486            return EclObject(special_sage_to_max[op](*[sr_to_max(o) for o in expr.operands()]))
     1487        elif op == tuple:
     1488            return maxima(expr.operands()).ecl()
    14821489        elif not (op in sage_op_dict):
    14831490            # Maxima does some simplifications automatically by default
    14841491            # so calling maxima(expr) can change the structure of expr
  • sage/symbolic/expression.pyx

    diff --git a/sage/symbolic/expression.pyx b/sage/symbolic/expression.pyx
    a b  
    74697469            factorial(n)
    74707470        """
    74717471        x = self
     7472        x = x.simplify_hypergeometric()
    74727473        x = x.simplify_factorial()
    74737474        x = x.simplify_trig()
    74747475        x = x.simplify_rational()
     
    74797480
    74807481    full_simplify = simplify_full
    74817482
     7483    def simplify_hypergeometric(self, algorithm='maxima'):
     7484        """
     7485        EXAMPLES::
     7486
     7487            sage: hypergeometric((5,4), (4,1,2,3),x).simplify_hypergeometric()         
     7488            1/144*x^2*hypergeometric((), (3, 4), x) + 1/3*x*hypergeometric((), (2, 3), x) + hypergeometric((), (1, 2), x)
     7489            sage: (2*hypergeometric((), (), x)).simplify_hypergeometric()
     7490            2*e^x
     7491        """
     7492        from sage.functions.hypergeometric import hypergeometric, closed_form
     7493        from sage.calculus.calculus import maxima
     7494        try:
     7495            op = self.operator()
     7496        except RuntimeError:
     7497            return self
     7498        ops = self.operands()
     7499        if op == hypergeometric:
     7500            if algorithm == 'maxima':
     7501                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)))
     7502            elif algorithm == 'sage':
     7503                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)))
     7504            else:
     7505                return NotImplementedError('unknown algorithm')
     7506        if not op:
     7507            return self
     7508        return op(*map(lambda o: o.simplify_hypergeometric(algorithm), ops))
     7509
     7510    hypergeometric_simplify = simplify_hypergeometric
     7511
    74827512    def simplify_trig(self,expand=True):
    74837513        r"""
    74847514        Optionally expands and then employs identities such as
  • sage/symbolic/expression_conversions.py

    diff --git a/sage/symbolic/expression_conversions.py b/sage/symbolic/expression_conversions.py
    a b  
    216216            return self.relation(ex, operator)
    217217        elif isinstance(operator, FDerivativeOperator):
    218218            return self.derivative(ex, operator)
     219        elif operator == tuple:
     220            return self.tuple(ex, operator)
    219221        else:
    220222            return self.composition(ex, operator)
    221223
     
    459461        return "%s %s %s"%(self(ex.lhs()), self.relation_symbols[operator],
    460462                           self(ex.rhs()))
    461463
     464    def tuple(self, ex, operator):
     465        x = map(self, ex.operands())
     466        X = ','.join(x)
     467        return '%s%s%s'%(self.interface._left_list_delim(), X, self.interface._right_list_delim())
     468
    462469    def derivative(self, ex, operator):
    463470        """
    464471        EXAMPLES::
     
    12381245            1.41421356237309...
    12391246        """
    12401247        f = operator
    1241         g = map(self, ex.operands())
     1248        try:
     1249            g = map(self, ex.operands())
     1250        except AttributeError, error:
     1251            from sage.functions.hypergeometric import hypergeometric
     1252            if f is hypergeometric:
     1253                ops = ex.operands()
     1254                return hypergeometric(*[map(lambda x: self.pyobject(None, x), ops[0]),
     1255                                        map(lambda x: self.pyobject(None, x), ops[1]),
     1256                                        self(ops[2])])
     1257            else:
     1258                raise AttributeError(error)
    12421259        try:
    12431260            return f(*g)
    12441261        except TypeError:
     
    12471264                return abs(*g) # special case
    12481265            else:
    12491266                return self.ff.fast_float_func(f, *g)
    1250             
     1267           
    12511268def fast_float(ex, *vars):
    12521269    """
    12531270    Returns an object which provides fast floating point evaluation of