| 1 | r""" |
| 2 | Hypergeometric Functions |
| 3 | |
| 4 | This module implements manipulation of infinite hypergeometric series |
| 5 | represented in standard parametric form (as $\,_pF_q$ functions). |
| 6 | |
| 7 | AUTHORS: |
| 8 | |
| 9 | - Fredrik Johansson (2010): initial version |
| 10 | |
| 11 | - Eviatar Bach (2013): major changes |
| 12 | |
| 13 | EXAMPLES: |
| 14 | |
| 15 | Examples 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 | |
| 25 | Simplification (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 | |
| 39 | Equality 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 | |
| 47 | Computing 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 | |
| 86 | Plotting:: |
| 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 | |
| 91 | Numeric 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 | |
| 106 | Sometimes it may be better to simplify before numeric evaluation (the error is |
| 107 | due 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 | |
| 116 | Conversions:: |
| 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 | |
| 125 | Arbitrary 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 | |
| 145 | from sage.rings.all import (Integer, ZZ, QQ, Infinity, binomial, |
| 146 | rising_factorial, factorial) |
| 147 | from sage.functions.other import sqrt, gamma, real_part |
| 148 | from sage.functions.log import exp, log |
| 149 | from sage.functions.trig import cos, sin |
| 150 | from sage.functions.hyperbolic import cosh, sinh |
| 151 | from sage.functions.other import erf |
| 152 | from sage.symbolic.constants import pi |
| 153 | from sage.symbolic.all import I |
| 154 | from sage.symbolic.function import BuiltinFunction, is_inexact |
| 155 | from sage.symbolic.ring import SR |
| 156 | from sage.structure.element import get_coercion_model |
| 157 | from sage.misc.latex import latex |
| 158 | from sage.misc.misc_c import prod |
| 159 | from mpmath import hyper |
| 160 | from sage.libs.mpmath import utils as mpmath_utils |
| 161 | from sage.symbolic.expression import Expression |
| 162 | from sage.calculus.functional import derivative |
| 163 | |
| 164 | |
| 165 | def 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 | |
| 195 | class 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 | |
| 667 | hypergeometric = Hypergeometric() |
| 668 | |
| 669 | |
| 670 | def 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()]) |