Ticket #2516: trac_2516_2.patch

File trac_2516_2.patch, 41.7 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 2c773352ebc9bc079a27fa06e03ca858bf60e39a
    # Parent  000a0f82d23d6cb39a00637051ec77a8fec2a99f
    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 $\,_pF_q$ 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),...
     21    (2, -1/2*I*sqrt(3) + 1/2, 1/2*I*sqrt(3) + 1/2), 1/2*I)
     22    sage: sum((-1)^x/((2*x + 1)*factorial(2*x + 1)), x, 0, oo)
     23    hypergeometric((1/2,), (3/2, 3/2), -1/4)
     24
     25Simplification (note that ``simplify_full`` does not yet call
     26``simplify_hypergeometric``)::
     27
     28    sage: hypergeometric([-2], [], x).simplify_hypergeometric()
     29    x^2 - 2*x + 1
     30    sage: hypergeometric([], [], x).simplify_hypergeometric()
     31    e^x
     32    sage: a = hypergeometric((hypergeometric((), (), x),), (),
     33    ....:                    hypergeometric((), (), x))
     34    sage: a.simplify_hypergeometric()
     35    1/((-e^x + 1)^e^x)
     36    sage: a.simplify_hypergeometric(algorithm='sage')
     37    (-e^x + 1)^(-e^x)
     38
     39Equality testing::
     40
     41    sage: bool(hypergeometric([], [], x).derivative(x) ==
     42    ....:      hypergeometric([], [], x))  # diff(e^x, x) == e^x
     43    True
     44    sage: bool(hypergeometric([], [], x) == hypergeometric([], [1], x))
     45    False
     46
     47Computing terms and series::
     48
     49    sage: z = var('z')
     50    sage: hypergeometric([], [], z).series(z, 0)
     51    Order(1)
     52    sage: hypergeometric([], [], z).series(z, 1)
     53    1 + Order(z)
     54    sage: hypergeometric([], [], z).series(z, 2)
     55    1 + 1*z + Order(z^2)
     56    sage: hypergeometric([], [], z).series(z, 3)
     57    1 + 1*z + 1/2*z^2 + Order(z^3)
     58
     59    sage: hypergeometric([-2], [], z).series(z, 3)
     60    1 + (-2)*z + 1*z^2
     61    sage: hypergeometric([-2], [], z).series(z, 6)
     62    1 + (-2)*z + 1*z^2
     63    sage: hypergeometric([-2], [], z).series(z, 6).is_terminating_series()
     64    True
     65    sage: hypergeometric([-2], [], z).series(z, 2)
     66    1 + (-2)*z + Order(z^2)
     67    sage: hypergeometric([-2], [], z).series(z, 2).is_terminating_series()
     68    False
     69
     70    sage: hypergeometric([1], [], z).series(z, 6)
     71    1 + 1*z + 1*z^2 + 1*z^3 + 1*z^4 + 1*z^5 + Order(z^6)
     72    sage: hypergeometric([], [1/2], -z^2/4).series(z, 11)
     73    1 + (-1/2)*z^2 + 1/24*z^4 + (-1/720)*z^6 + 1/40320*z^8 +...
     74    (-1/3628800)*z^10 + Order(z^11)
     75
     76    sage: hypergeometric([1], [5], x).series(x, 5)
     77    1 + 1/5*x + 1/30*x^2 + 1/210*x^3 + 1/1680*x^4 + Order(x^5)
     78
     79    sage: sum(hypergeometric([1, 2], [3], 1/3).terms(6)).n()
     80    1.29788359788360
     81    sage: hypergeometric([1, 2], [3], 1/3).n()
     82    1.29837194594696
     83    sage: hypergeometric([], [], x).series(x, 20)(x=1).n() == e.n()
     84    True
     85
     86Plotting::
     87
     88    sage: plot(hypergeometric([1, 1], [3, 3, 3], x), x, -30, 30)
     89    sage: complex_plot(hypergeometric([x], [], 2), (-1, 1), (-1, 1))
     90
     91Numeric evaluation::
     92
     93    sage: hypergeometric([1], [], 1/10).n()  # geometric series
     94    1.11111111111111
     95    sage: hypergeometric([], [], 1).n()  # e
     96    2.71828182845905
     97    sage: hypergeometric([], [], 3., hold=True)
     98    hypergeometric((), (), 3.00000000000000)
     99    sage: hypergeometric([1, 2, 3], [4, 5, 6], 1/2).n()
     100    1.02573619590134
     101    sage: hypergeometric([1, 2, 3], [4, 5, 6], 1/2).n(digits=30)
     102    1.02573619590133865036584139535
     103    sage: hypergeometric([5 - 3*I], [3/2, 2 + I, sqrt(2)], 4 + I).n()
     104    5.52605111678805 - 7.86331357527544*I
     105    sage: hypergeometric((10, 10), (50,), 2.)
     106    -1705.75733163554 - 356.749986056024*I
     107
     108Conversions::
     109
     110    sage: maxima(hypergeometric([1, 1, 1], [3, 3, 3], x))
     111    hypergeometric([1,1,1],[3,3,3],x)
     112    sage: hypergeometric((5, 4), (4, 4), 3)._sympy_()
     113    hyper((5, 4), (4, 4), 3)
     114    sage: hypergeometric((5, 4), (4, 4), 3)._mathematica_init_()
     115    'HypergeometricPFQ[{5,4},{4,4},3]'
     116
     117Arbitrary level of nesting for conversions::
     118
     119    sage: maxima(nest(lambda y: hypergeometric([y], [], x), 3, 1))
     120    1/(1-x)^(1/(1-x)^(1/(1-x)))
     121    sage: maxima(nest(lambda y: hypergeometric([y], [3], x), 3, 1))._sage_()
     122    hypergeometric((hypergeometric((hypergeometric((1,), (3,), x),), (3,),...
     123    x),), (3,), x)
     124    sage: nest(lambda y: hypergeometric([y], [], x), 3, 1)._mathematica_init_()
     125    'HypergeometricPFQ[{HypergeometricPFQ[{HypergeometricPFQ[{1},{},x]},...
     126"""
     127#*****************************************************************************
     128#       Copyright (C) 2010 Fredrik Johansson <fredrik.johansson@gmail.com>
     129#       Copyright (C) 2013 Eviatar Bach <eviatarbach@gmail.com>
     130#
     131#  Distributed under the terms of the GNU General Public License (GPL)
     132#  as published by the Free Software Foundation; either version 2 of
     133#  the License, or (at your option) any later version.
     134#                  http://www.gnu.org/licenses/
     135#*****************************************************************************
     136
     137from sage.rings.all import (Integer, ZZ, QQ, Infinity, binomial,
     138                            rising_factorial, factorial)
     139from sage.functions.other import sqrt, gamma, real_part
     140from sage.functions.log import exp, log
     141from sage.functions.trig import cos, sin
     142from sage.functions.hyperbolic import cosh, sinh
     143from sage.functions.other import erf
     144from sage.symbolic.constants import pi
     145from sage.symbolic.all import I
     146from sage.symbolic.function import BuiltinFunction, is_inexact
     147from sage.symbolic.ring import SR
     148from sage.structure.element import get_coercion_model
     149from sage.misc.latex import latex
     150from sage.misc.misc_c import prod
     151from sage.libs.mpmath import utils as mpmath_utils
     152from sage.symbolic.expression import Expression
     153from sage.calculus.functional import derivative
     154
     155
     156def rational_param_as_tuple(x):
     157    """
     158    Utility function for converting rational pFq parameters to
     159    tuples (which mpmath handles more efficiently).
     160
     161    EXAMPLES::
     162
     163        sage: from sage.functions.hypergeometric import rational_param_as_tuple
     164        sage: rational_param_as_tuple(1/2)
     165        (1, 2)
     166        sage: rational_param_as_tuple(3)
     167        3
     168        sage: rational_param_as_tuple(pi)
     169        pi
     170
     171    """
     172    try:
     173        x = x.pyobject()
     174    except AttributeError:
     175        pass
     176    try:
     177        if x.parent() is QQ:
     178            p = int(x.numer())
     179            q = int(x.denom())
     180            return p, q
     181    except AttributeError:
     182        pass
     183    return x
     184
     185
     186class Hypergeometric(BuiltinFunction):
     187    r"""
     188    Represents a (formal) generalized infinite hypergeometric series. It is
     189    defined as `\,{}_pF_q(a_1,\ldots,a_p;b_1,\ldots,b_q;z) = \sum_{n=0}^\infty
     190    \frac{(a_1)_n\dots(a_p)_n}{(b_1)_n\dots(b_q)_n} \, \frac{z^n}{n!},`
     191
     192    where `(x)_n` is the rising factorial.
     193
     194    INPUT:
     195
     196    - ``a`` -- a list or tuple of parameters
     197    - ``b`` -- a list or tuple of parameters
     198    - ``z`` -- a number or symbolic expression
     199
     200    EXAMPLES::
     201
     202        sage: hypergeometric([], [], 1)
     203        hypergeometric((), (), 1)
     204        sage: hypergeometric([], [1], 1)
     205        hypergeometric((), (1,), 1)
     206        sage: hypergeometric([2, 3], [1], 1)
     207        hypergeometric((2, 3), (1,), 1)
     208        sage: hypergeometric([], [], x)
     209        hypergeometric((), (), x)
     210        sage: hypergeometric([x], [], x^2)
     211        hypergeometric((x,), (), x^2)
     212
     213    The only simplification that is done automatically is returning 1 if ``z``
     214    is 0::
     215
     216        sage: hypergeometric([], [], 0)
     217        1
     218
     219    For other simplifications use the ``simplify_hypergeometric`` method.
     220    """
     221    def __init__(self):
     222        BuiltinFunction.__init__(self, 'hypergeometric', nargs=3,
     223                                 conversions={'mathematica':
     224                                              'HypergeometricPFQ',
     225                                              'maxima': 'hypergeometric',
     226                                              'sympy': 'hyper'})
     227
     228    def __call__(self, a, b, z, **kwargs):
     229        return BuiltinFunction.__call__(self,
     230                                        SR._force_pyobject(a),
     231                                        SR._force_pyobject(b),
     232                                        z, **kwargs)
     233
     234    def _print_latex_(self, a, b, z):
     235        r"""
     236        TESTS::
     237
     238            sage: latex(hypergeometric([1, 1], [2], -1))
     239            \,_2F_1\left(\begin{matrix} 1,1 \\ 2 \end{matrix} ; -1 \right)
     240
     241        """
     242        aa = ",".join(latex(c) for c in a)
     243        bb = ",".join(latex(c) for c in b)
     244        z = latex(z)
     245        return (r"\,_{}F_{}\left(\begin{{matrix}} {} \\ {} \end{{matrix}} ; "
     246                r"{} \right)").format(len(a), len(b), aa, bb, z)
     247
     248    def _eval_(self, a, b, z, **kwargs):
     249        coercion_model = get_coercion_model()
     250        co = reduce(lambda x, y: coercion_model.canonical_coercion(x, y)[0],
     251                    a + b + (z,))
     252        if is_inexact(co) and not isinstance(co, Expression):
     253            from sage.structure.coerce import parent
     254            return self._evalf_(a, b, z, parent=parent(co))
     255        if not isinstance(z, Expression) and z == 0:  # Expression is excluded
     256            return Integer(1)                         # to avoid call to Maxima
     257        return
     258
     259    def _evalf_(self, a, b, z, parent):
     260        """
     261        TESTS::
     262
     263            sage: hypergeometric([1, 1], [2], -1).n()
     264            0.693147180559945
     265            sage: hypergeometric([], [], RealField(100)(1))
     266            2.7182818284590452353602874714
     267
     268        """
     269        from mpmath import hyper
     270        aa = [rational_param_as_tuple(c) for c in a]
     271        bb = [rational_param_as_tuple(c) for c in b]
     272        return mpmath_utils.call(hyper, aa, bb, z, parent=parent)
     273
     274    def _tderivative_(self, a, b, z, *args, **kwargs):
     275        """
     276        EXAMPLES::
     277
     278            sage: hypergeometric([1/3, 2/3], [5], x^2).diff(x)
     279            4/45*x*hypergeometric((4/3, 5/3), (6,), x^2)
     280            sage: hypergeometric([1, 2], [x], 2).diff(x)
     281            Traceback (most recent call last):
     282            ...
     283            NotImplementedError: derivative of hypergeometric function with...
     284            respect to parameters. Try calling .simplify_hypergeometric()...
     285            first.
     286            sage: hypergeometric([1/3, 2/3], [5], 2).diff(x)
     287            0
     288        """
     289        diff_param = kwargs['diff_param']
     290        if diff_param in hypergeometric(a, b, 1).variables():  # ignore z
     291            raise NotImplementedError("derivative of hypergeometric function "
     292                                      "with respect to parameters. Try calling"
     293                                      " .simplify_hypergeometric() first.")
     294        t = (reduce(lambda x, y: x * y, a, 1) *
     295             reduce(lambda x, y: x / y, b, Integer(1)))
     296        return (t * derivative(z, diff_param) *
     297                hypergeometric([c + 1 for c in a], [c + 1 for c in b], z))
     298
     299    class EvaluationMethods:
     300        def _fast_float_(cls, self, *args):
     301            """
     302            Do not support the old ``fast_float``
     303
     304            OUTPUT:
     305
     306            This method raises ``NotImplementedError``; use the newer
     307            ``fast_callable`` implementation
     308
     309            EXAMPLES::
     310
     311                sage: f = hypergeometric([], [], x)
     312                sage: f._fast_float_()
     313                Traceback (most recent call last):
     314                ...
     315                NotImplementedError
     316            """
     317            raise NotImplementedError
     318
     319        def _fast_callable_(cls, self, a, b, z, etb):
     320            """
     321            Override the ``fast_callable`` method
     322
     323            OUTPUT:
     324
     325            A :class:`~sage.ext.fast_callable.ExpressionCall` representing the
     326            hypergeometric function in the expression tree
     327
     328            EXAMPLES::
     329
     330                sage: h = hypergeometric([], [], x)
     331                sage: from sage.ext.fast_callable import ExpressionTreeBuilder
     332                sage: etb = ExpressionTreeBuilder(vars=['x'])
     333                sage: h._fast_callable_(etb)
     334                {CommutativeRings.element_class}(v_0)
     335
     336                sage: y = var('y')
     337                sage: fast_callable(hypergeometric([y], [], x),
     338                ....:               vars=[x, y])(3, 4)
     339                hypergeometric((4,), (), 3)
     340            """
     341            return etb.call(self, *map(etb.var, etb._vars))
     342
     343        def sorted_parameters(cls, self, a, b, z):
     344            """
     345            Return with parameters sorted in a canonical order
     346
     347            EXAMPLES::
     348
     349                sage: hypergeometric([2, 1, 3], [5, 4],
     350                ....:                1/2).sorted_parameters()
     351                hypergeometric((1, 2, 3), (4, 5), 1/2)
     352
     353            """
     354            return hypergeometric(sorted(a), sorted(b), z)
     355
     356        def eliminate_parameters(cls, self, a, b, z):
     357            """
     358            Eliminate repeated parameters by pairwise cancellation of identical
     359            terms in ``a`` and ``b``
     360
     361            EXAMPLES::
     362
     363                sage: hypergeometric([1, 1, 2, 5], [5, 1, 4],
     364                ....:                1/2).eliminate_parameters()
     365                hypergeometric((1, 2), (4,), 1/2)
     366                sage: hypergeometric([x], [x], x).eliminate_parameters()
     367                hypergeometric((), (), x)
     368                sage: hypergeometric((5, 4), (4, 4), 3).eliminate_parameters()
     369                hypergeometric((5,), (4,), 3)
     370
     371            """
     372            aa = list(a)  # tuples are immutable
     373            bb = list(b)
     374            p = pp = len(aa)
     375            q = qq = len(bb)
     376            i = 0
     377            while i < qq and aa:
     378                bbb = bb[i]
     379                if bbb in aa:
     380                    aa.remove(bbb)
     381                    bb.remove(bbb)
     382                    pp -= 1
     383                    qq -= 1
     384                else:
     385                    i += 1
     386            if (pp, qq) != (p, q):
     387                return hypergeometric(aa, bb, z)
     388            return self
     389
     390        def is_termwise_finite(cls, self, a, b, z):
     391            """
     392            Determine whether all terms of self are finite. Any infinite
     393            terms or ambiguous terms beyond the first zero, if one exists,
     394            are ignored.
     395
     396            Ambiguous cases (where a term is the product of both zero
     397            and an infinity) are not considered finite.
     398
     399            EXAMPLES::
     400
     401                sage: hypergeometric([2], [3, 4], 5).is_termwise_finite()
     402                True
     403                sage: hypergeometric([2], [-3, 4], 5).is_termwise_finite()
     404                False
     405                sage: hypergeometric([-2], [-3, 4], 5).is_termwise_finite()
     406                True
     407                sage: hypergeometric([-3], [-3, 4],
     408                ....:                5).is_termwise_finite()  # ambiguous
     409                False
     410
     411                sage: hypergeometric([0], [-1], 5).is_termwise_finite()
     412                True
     413                sage: hypergeometric([0], [0],
     414                ....:                5).is_termwise_finite()  # ambiguous
     415                False
     416                sage: hypergeometric([1], [2], Infinity).is_termwise_finite()
     417                False
     418                sage: (hypergeometric([0], [0], Infinity)
     419                ....:  .is_termwise_finite())  # ambiguous
     420                False
     421                sage: (hypergeometric([0], [], Infinity)
     422                ....:  .is_termwise_finite())  # ambiguous
     423                False
     424
     425            """
     426            if z == 0:
     427                return 0 not in b
     428            if abs(z) == Infinity:
     429                return False
     430            if abs(z) == Infinity:
     431                return False
     432            for bb in b:
     433                if bb in ZZ and bb <= 0:
     434                    if any((aa in ZZ) and (bb < aa <= 0) for aa in a):
     435                        continue
     436                    return False
     437            return True
     438
     439        def is_terminating(cls, self, a, b, z):
     440            """
     441            Determine whether the series represented by self terminates
     442            after a finite number of terms, i.e. whether any of the
     443            numerator parameters are nonnegative integers (with no
     444            preceding nonnegative denominator parameters), or z = 0.
     445
     446            If terminating, the series represents a polynomial of z
     447
     448            EXAMPLES::
     449
     450                sage: hypergeometric([1, 2], [3, 4], x).is_terminating()
     451                False
     452                sage: hypergeometric([1, -2], [3, 4], x).is_terminating()
     453                True
     454                sage: hypergeometric([1, -2], [], x).is_terminating()
     455                True
     456
     457            """
     458            if z == 0:
     459                return True
     460            for aa in a:
     461                if (aa in ZZ) and (aa <= 0):
     462                    return self.is_termwise_finite()
     463            return False
     464
     465        def is_absolutely_convergent(cls, self, a, b, z):
     466            """
     467            Determine whether self converges absolutely as an infinite series.
     468            False is returned if not all terms are finite.
     469
     470            EXAMPLES:
     471
     472            Degree giving infinite radius of convergence::
     473
     474                sage: hypergeometric([2, 3], [4, 5],
     475                ....:                6).is_absolutely_convergent()
     476                True
     477                sage: hypergeometric([2, 3], [-4, 5],
     478                ....:                6).is_absolutely_convergent()  # undefined
     479                False
     480                sage: (hypergeometric([2, 3], [-4, 5], Infinity)
     481                ....:  .is_absolutely_convergent())  # undefined
     482                False
     483
     484            Ordinary geometric series (unit radius of convergence)::
     485
     486                sage: hypergeometric([1], [], 1/2).is_absolutely_convergent()
     487                True
     488                sage: hypergeometric([1], [], 2).is_absolutely_convergent()
     489                False
     490                sage: hypergeometric([1], [], 1).is_absolutely_convergent()
     491                False
     492                sage: hypergeometric([1], [], -1).is_absolutely_convergent()
     493                False
     494                sage: hypergeometric([1], [], -1).n()  # Sum still exists
     495                0.500000000000000
     496
     497            Degree $p = q+1$ (unit radius of convergence)::
     498
     499                sage: hypergeometric([2, 3], [4], 6).is_absolutely_convergent()
     500                False
     501                sage: hypergeometric([2, 3], [4], 1).is_absolutely_convergent()
     502                False
     503                sage: hypergeometric([2, 3], [5], 1).is_absolutely_convergent()
     504                False
     505                sage: hypergeometric([2, 3], [6], 1).is_absolutely_convergent()
     506                True
     507                sage: hypergeometric([-2, 3], [4],
     508                ....:                5).is_absolutely_convergent()
     509                True
     510                sage: hypergeometric([2, -3], [4],
     511                ....:                5).is_absolutely_convergent()
     512                True
     513                sage: hypergeometric([2, -3], [-4],
     514                ....:                5).is_absolutely_convergent()
     515                True
     516                sage: hypergeometric([2, -3], [-1],
     517                ....:                5).is_absolutely_convergent()
     518                False
     519
     520            Degree giving zero radius of convergence::
     521
     522                sage: hypergeometric([1, 2, 3], [4],
     523                ....:                2).is_absolutely_convergent()
     524                False
     525                sage: hypergeometric([1, 2, 3], [4],
     526                ....:                1/2).is_absolutely_convergent()
     527                False
     528                sage: (hypergeometric([1, 2, -3], [4], 1/2)
     529                ....:  .is_absolutely_convergent())  # polynomial
     530                True
     531
     532            """
     533            p, q = len(a), len(b)
     534            if not self.is_termwise_finite():
     535                return False
     536            if p <= q:
     537                return True
     538            if self.is_terminating():
     539                return True
     540            if p == q + 1:
     541                if abs(z) < 1:
     542                    return True
     543                if abs(z) == 1:
     544                    if real_part(sum(b) - sum(a)) > 0:
     545                        return True
     546            return False
     547
     548        def terms(cls, self, a, b, z, n=None):
     549            """
     550            Generate the terms of self (optionally only n terms).
     551
     552            EXAMPLES::
     553
     554                sage: list(hypergeometric([-2, 1], [3, 4], x).terms())
     555                [1, -1/6*x, 1/120*x^2]
     556                sage: list(hypergeometric([-2, 1], [3, 4], x).terms(2))
     557                [1, -1/6*x]
     558                sage: list(hypergeometric([-2, 1], [3, 4], x).terms(0))
     559                []
     560
     561            """
     562            if n is None:
     563                n = Infinity
     564            t = Integer(1)
     565            k = 1
     566            while k <= n:
     567                yield t
     568                for aa in a:
     569                    t *= (aa + k - 1)
     570                for bb in b:
     571                    t /= (bb + k - 1)
     572                t *= z
     573                if t == 0:
     574                    break
     575                t /= k
     576                k += 1
     577
     578        def deflated(cls, self, a, b, z):
     579            """
     580            Rewrite as a linear combination of functions of strictly lower
     581            degree by eliminating all parameters ``a[i]`` and ``b[j]`` such
     582            that ``a[i]`` = ``b[i]`` + ``m`` for nonnegative integer ``m``.
     583
     584            EXAMPLES::
     585
     586                sage: x = hypergeometric([5], [4], 3)
     587                sage: y = x.deflated()
     588                sage: y
     589                7/4*hypergeometric((), (), 3)
     590                sage: x.n(); y.n()
     591                35.1496896155784
     592                35.1496896155784
     593
     594                sage: x = hypergeometric([6, 1], [3, 4, 5], 10)
     595                sage: y = x.deflated()
     596                sage: y
     597                hypergeometric((1,), (4, 5), 10) +...
     598                1/2*hypergeometric((2,), (5, 6), 10) +...
     599                1/12*hypergeometric((3,), (6, 7), 10) +...
     600                1/252*hypergeometric((4,), (7, 8), 10)
     601                sage: x.n(); y.n()
     602                2.87893612686782
     603                2.87893612686782
     604
     605                sage: x = hypergeometric([6, 7], [3, 4, 5], 10)
     606                sage: y = x.deflated()
     607                sage: y
     608                hypergeometric((), (5,), 10) +...
     609                5*hypergeometric((), (6,), 10) +...
     610                19/3*hypergeometric((), (7,), 10) +...
     611                181/63*hypergeometric((), (8,), 10) +...
     612                265/504*hypergeometric((), (9,), 10) +...
     613                25/648*hypergeometric((), (10,), 10) +...
     614                25/27216*hypergeometric((), (11,), 10)
     615                sage: x.n(); y.n()
     616                63.0734110716969
     617                63.0734110716969
     618
     619            """
     620            return sum(map(prod, self._deflated()))
     621
     622        def _deflated(cls, self, a, b, z):
     623            new = self.eliminate_parameters()
     624            aa = new.operands()[0].operands()
     625            bb = new.operands()[1].operands()
     626            for i, aaa in enumerate(aa):
     627                for j, bbb in enumerate(bb):
     628                    m = aaa - bbb
     629                    if m in ZZ and m > 0:
     630                        aaaa = aa[:i] + aa[i + 1:]
     631                        bbbb = bb[:j] + bb[j + 1:]
     632                        terms = []
     633                        for k in xrange(m + 1):
     634                            # TODO: could rewrite prefactors as recurrence
     635                            term = binomial(m, k)
     636                            for c in aaaa:
     637                                term *= rising_factorial(c, k)
     638                            for c in bbbb:
     639                                term /= rising_factorial(c, k)
     640                            term *= z ** k
     641                            term /= rising_factorial(aaa - m, k)
     642                            F = hypergeometric([c + k for c in aaaa],
     643                                               [c + k for c in bbbb], z)
     644                            unique = []
     645                            counts = []
     646                            for c, f in F._deflated():
     647                                if f in unique:
     648                                    counts[unique.index(f)] += c
     649                                else:
     650                                    unique.append(f)
     651                                    counts.append(c)
     652                            Fterms = zip(counts, unique)
     653                            terms += [(term * termG, G) for (termG, G) in
     654                                      Fterms]
     655                        return terms
     656            return ((1, new),)
     657
     658hypergeometric = Hypergeometric()
     659
     660
     661def closed_form(hyp):
     662    """
     663    Try to evaluate self in closed form using elementary
     664    (and other simple) functions.
     665
     666    It may be necessary to call :meth:`deflated` first to
     667    find some closed forms.
     668
     669    EXAMPLES::
     670
     671        sage: from sage.functions.hypergeometric import closed_form
     672        sage: var('a b c z')
     673        (a, b, c, z)
     674        sage: closed_form(hypergeometric([1], [], 1 + z))
     675        -1/z
     676        sage: closed_form(hypergeometric([], [], 1 + z))
     677        e^(z + 1)
     678        sage: closed_form(hypergeometric([], [1/2], 4))
     679        cosh(4)
     680        sage: closed_form(hypergeometric([], [3/2], 4))
     681        1/4*sinh(4)
     682        sage: closed_form(hypergeometric([], [5/2], 4))
     683        -3/64*sinh(4) + 3/16*cosh(4)
     684        sage: closed_form(hypergeometric([], [-3/2], 4))
     685        -4*sinh(4) + 19/3*cosh(4)
     686        sage: closed_form(hypergeometric([-3, 1], [var('a')], z))
     687        -6*z^3/((a + 1)*(a + 2)*a) + 6*z^2/((a + 1)*a) - 3*z/a + 1
     688        sage: closed_form(hypergeometric([-3, 1/3], [-4], z))
     689        7/162*z^3 + 1/9*z^2 + 1/4*z + 1
     690        sage: closed_form(hypergeometric([], [], z))
     691        e^z
     692        sage: closed_form(hypergeometric([a], [], z))
     693        (-z + 1)^(-a)
     694        sage: closed_form(hypergeometric([1, 1, 2], [1, 1], z))
     695        (z - 1)^(-2)
     696        sage: closed_form(hypergeometric([2, 3], [1], x))
     697        -1/(x - 1)^3 + 3*x/(x - 1)^4
     698        sage: closed_form(hypergeometric([1/2], [3/2], -5))
     699        1/10*sqrt(pi)*sqrt(5)*erf(sqrt(5))
     700        sage: closed_form(hypergeometric([2], [5], 3))
     701        4
     702        sage: closed_form(hypergeometric([2], [5], 5))
     703        48/625*e^5 + 612/625
     704        sage: closed_form(hypergeometric([1/2, 7/2], [3/2], z))
     705        1/sqrt(-z + 1) + 2/3*z/(-z + 1)^(3/2) + 1/5*z^2/(-z + 1)^(5/2)
     706        sage: closed_form(hypergeometric([1/2, 1], [2], z))
     707        -2*(sqrt(-z + 1) - 1)/z
     708        sage: closed_form(hypergeometric([1, 1], [2], z))
     709        -log(-z + 1)/z
     710        sage: closed_form(hypergeometric([1, 1], [3], z))
     711        -2*((z - 1)*log(-z + 1)/z - 1)/z
     712        sage: closed_form(hypergeometric([1, 1, 1], [2, 2], x))
     713        hypergeometric((1, 1, 1), (2, 2), x)
     714    """
     715    if hyp.is_terminating():
     716        return sum(hyp.terms())
     717
     718    new = hyp.eliminate_parameters()
     719
     720    def _closed_form(hyp):
     721        a, b, z = hyp.operands()
     722        a, b = a.operands(), b.operands()
     723        p, q = len(a), len(b)
     724
     725        if z == 0:
     726            return Integer(1)
     727        if p == q == 0:
     728            return exp(z)
     729        if p == 1 and q == 0:
     730            return (1 - z) ** (-a[0])
     731
     732        if p == 0 and q == 1:
     733            # TODO: make this require only linear time
     734            def _0f1(b, z):
     735                F12 = cosh(2 * sqrt(z))
     736                F32 = sinh(2 * sqrt(z)) / (2 * sqrt(z))
     737                if 2 * b == 1:
     738                    return F12
     739                if 2 * b == 3:
     740                    return F32
     741                if 2 * b > 3:
     742                    return ((b - 2) * (b - 1) / z * (_0f1(b - 2, z) -
     743                            _0f1(b - 1, z)))
     744                if 2 * b < 1:
     745                    return (_0f1(b + 1, z) + z / (b * (b + 1)) *
     746                            _0f1(b + 2, z))
     747                raise ValueError
     748            # Can evaluate 0F1 in terms of elementary functions when
     749            # the parameter is a half-integer
     750            if 2 * b[0] in ZZ and b[0] not in ZZ:
     751                return _0f1(b[0], z)
     752
     753        # Confluent hypergeometric function
     754        if p == 1 and q == 1:
     755            aa, bb = a[0], b[0]
     756            if aa * 2 == 1 and bb * 2 == 3:
     757                t = sqrt(-z)
     758                return sqrt(pi) / 2 * erf(t) / t
     759            if a == 1 and b == 2:
     760                return (exp(z) - 1) / z
     761            n, m = aa, bb
     762            if n in ZZ and m in ZZ and m > 0 and n > 0:
     763                rf = rising_factorial
     764                if m <= n:
     765                    return (exp(z) * sum(rf(m - n, k) * (-z) ** k /
     766                            factorial(k) / rf(m, k) for k in
     767                            xrange(n - m + 1)))
     768                else:
     769                    T = sum(rf(n - m + 1, k) * z ** k /
     770                            (factorial(k) * rf(2 - m, k)) for k in
     771                            xrange(m - n))
     772                    U = sum(rf(1 - n, k) * (-z) ** k /
     773                            (factorial(k) * rf(2 - m, k)) for k in
     774                            xrange(n))
     775                    return (factorial(m - 2) * rf(1 - m, n) *
     776                            z ** (1 - m) / factorial(n - 1) *
     777                            (T - exp(z) * U))
     778
     779        if p == 2 and q == 1:
     780            R12 = QQ('1/2')
     781            R32 = QQ('3/2')
     782
     783            def _2f1(a, b, c, z):
     784                """
     785                Evaluation of 2F1(a, b, c, z), assuming a, b, c positive
     786                integers or half-integers
     787                """
     788                if b == c:
     789                    return (1 - z) ** (-a)
     790                if a == c:
     791                    return (1 - z) ** (-b)
     792                if a == 0 or b == 0:
     793                    return Integer(1)
     794                if a > b:
     795                    a, b = b, a
     796                if b >= 2:
     797                    F1 = _2f1(a, b - 1, c, z)
     798                    F2 = _2f1(a, b - 2, c, z)
     799                    q = (b - 1) * (z - 1)
     800                    return (((c - 2 * b + 2 + (b - a - 1) * z) * F1 +
     801                            (b - c - 1) * F2) / q)
     802                if c > 2:
     803                    # how to handle this case?
     804                    if a - c + 1 == 0 or b - c + 1 == 0:
     805                        raise NotImplementedError
     806                    F1 = _2f1(a, b, c - 1, z)
     807                    F2 = _2f1(a, b, c - 2, z)
     808                    r1 = (c - 1) * (2 - c - (a + b - 2 * c + 3) * z)
     809                    r2 = (c - 1) * (c - 2) * (1 - z)
     810                    q = (a - c + 1) * (b - c + 1) * z
     811                    return (r1 * F1 + r2 * F2) / q
     812
     813                if (a, b, c) == (R12, 1, 2):
     814                    return (2 - 2 * sqrt(1 - z)) / z
     815                if (a, b, c) == (1, 1, 2):
     816                    return -log(1 - z) / z
     817                if (a, b, c) == (1, R32, R12):
     818                    return (1 + z) / (1 - z) ** 2
     819                if (a, b, c) == (1, R32, 2):
     820                    return 2 * (1 / sqrt(1 - z) - 1) / z
     821                if (a, b, c) == (R32, 2, R12):
     822                    return (1 + 3 * z) / (1 - z) ** 3
     823                if (a, b, c) == (R32, 2, 1):
     824                    return (2 + z) / (2 * (sqrt(1 - z) * (1 - z) ** 2))
     825                if (a, b, c) == (2, 2, 1):
     826                    return (1 + z) / (1 - z) ** 3
     827                raise NotImplementedError
     828            aa, bb = a
     829            cc, = b
     830            if z == 1:
     831                return (gamma(cc) * gamma(cc - aa - bb) / gamma(cc - aa) /
     832                        gamma(cc - bb))
     833            if ((aa * 2) in ZZ and (bb * 2) in ZZ and (cc * 2) in ZZ and
     834                aa > 0 and bb > 0 and cc > 0):
     835                try:
     836                    return _2f1(aa, bb, cc, z)
     837                except NotImplementedError:
     838                    pass
     839        return hyp
     840    return sum([coeff * _closed_form(pfq) for coeff, pfq in new._deflated()])
  • 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  
    74797479
    74807480    full_simplify = simplify_full
    74817481
     7482    def simplify_hypergeometric(self, algorithm='maxima'):
     7483        """
     7484        Simplify an expression containing hypergeometric functions
     7485
     7486        INPUT:
     7487
     7488        - ``self`` -- symbolic expression
     7489 
     7490        - ``algorithm`` -- (default: ``'maxima'``) the algorithm to use for
     7491          for simplification. Implemented are ``'maxima'``, which uses Maxima's
     7492          ``hgfred`` function, and ``'sage'``, which uses an algorithm
     7493          implemented in the hypergeometric module
     7494
     7495        ALIAS: :meth:`hypergeometric_simplify` and
     7496        :meth:`simplify_hypergeometric` are the same
     7497
     7498        EXAMPLES::
     7499
     7500            sage: hypergeometric((5, 4), (4, 1, 2, 3),
     7501            ....:                x).simplify_hypergeometric()
     7502            1/144*x^2*hypergeometric((), (3, 4), x) +...
     7503            1/3*x*hypergeometric((), (2, 3), x) + hypergeometric((), (1, 2), x)
     7504            sage: (2*hypergeometric((), (), x)).simplify_hypergeometric()
     7505            2*e^x
     7506            sage: (nest(lambda y: hypergeometric([y], [1], x), 3, 1)
     7507            ....:  .simplify_hypergeometric())
     7508            laguerre(-laguerre(-e^x, x), x)
     7509            sage: (nest(lambda y: hypergeometric([y], [1], x), 3, 1)
     7510            ....:  .simplify_hypergeometric(algorithm='sage'))
     7511            hypergeometric((hypergeometric((e^x,), (1,), x),), (1,), x)
     7512
     7513        """
     7514        from sage.functions.hypergeometric import hypergeometric, closed_form
     7515        from sage.calculus.calculus import maxima
     7516        try:
     7517            op = self.operator()
     7518        except RuntimeError:
     7519            return self
     7520        ops = self.operands()
     7521        if op == hypergeometric:
     7522            if algorithm == 'maxima':
     7523                return (self.parent()
     7524                        (maxima.hgfred(map(lambda o: o.simplify_hypergeometric(algorithm),
     7525                                           ops[0].operands()),
     7526                                       map(lambda o: o.simplify_hypergeometric(algorithm),
     7527                                           ops[1].operands()),
     7528                                       ops[2].simplify_hypergeometric(algorithm))))
     7529            elif algorithm == 'sage':
     7530                return (closed_form
     7531                        (hypergeometric(map(lambda o: o.simplify_hypergeometric(algorithm),
     7532                                            ops[0].operands()),
     7533                                        map(lambda o: o.simplify_hypergeometric(algorithm),
     7534                                            ops[1].operands()),
     7535                                        ops[2].simplify_hypergeometric(algorithm))))
     7536            else:
     7537                return NotImplementedError('unknown algorithm')
     7538        if not op:
     7539            return self
     7540        return op(*map(lambda o: o.simplify_hypergeometric(algorithm), ops))
     7541
     7542    hypergeometric_simplify = simplify_hypergeometric
     7543
    74827544    def simplify_trig(self,expand=True):
    74837545        r"""
    74847546        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)
    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):
     465        """
     466        EXAMPLES::
     467
     468            sage: from sage.symbolic.expression_conversions import InterfaceInit
     469            sage: m = InterfaceInit(maxima)
     470            sage: t = SR._force_pyobject((3, 4, e^x))
     471            sage: m.tuple(t)
     472            '[3,4,exp(x)]'
     473        """
     474        x = map(self, ex.operands())
     475        X = ','.join(x)
     476        return '%s%s%s'%(self.interface._left_list_delim(), X, self.interface._right_list_delim())
     477
    462478    def derivative(self, ex, operator):
    463479        """
    464480        EXAMPLES::
     
    12471263                return abs(*g) # special case
    12481264            else:
    12491265                return self.ff.fast_float_func(f, *g)
    1250            
     1266
    12511267def fast_float(ex, *vars):
    12521268    """
    12531269    Returns an object which provides fast floating point evaluation of
     
    14081424        """
    14091425        return self.etb.call(function, *ex.operands())
    14101426
     1427    def tuple(self, ex):
     1428        r"""
     1429        Given a symbolic tuple, return its elements as a Python list
     1430
     1431        EXAMPLES::
     1432
     1433            sage: from sage.ext.fast_callable import ExpressionTreeBuilder
     1434            sage: etb = ExpressionTreeBuilder(vars=['x'])
     1435            sage: SR._force_pyobject((2, 3, x^2))._fast_callable_(etb)
     1436            [2, 3, x^2]
     1437        """
     1438        return ex.operands()
    14111439
    14121440def fast_callable(ex, etb):
    14131441    """