Ticket #2516: trac_2516.patch

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