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