sympy/core/evalf.py

   1 """
   2 Adaptive numerical evaluation of SymPy expressions, using mpmath
   3 for mathematical functions.
   4 """
   5
   6 from sympy.mpmath.lib import (from_int, from_rational, fpi, fzero, fcmp,
   7     normalize, bitcount, round_nearest, to_str, fpow, fone, fpowi, fe,
   8     fnone, fhalf, fcos, fsin, flog, fatan, fmul, fneg, to_float, fshift,
   9     fabs, fatan2, fadd, fdiv, flt, dps_to_prec, prec_to_dps, fpos, from_float,
  10     fnone, to_int, flt, fexp, fsqrt)
  11
  12 import sympy.mpmath.libmpc as libmpc
  13 from sympy.mpmath import mpf, mpc, quadts, quadosc, mp, make_mpf
  14 from sympy.mpmath.specfun import mpf_gamma
  15 from sympy.mpmath.lib import MP_BASE, from_man_exp
  16 from sympy.mpmath.calculus import shanks_extrapolation, richardson_extrapolation
  17
  18
  19 import math
  20
  21 from basic import Basic, C, S
  22 from function import Function
  23 from sympify import sympify
  24
  25 LG10 = math.log(10,2)
  26
  27 # Used in a few places as placeholder values to denote exponents and
  28 # precision levels, e.g. of exact numbers. Must be careful to avoid
  29 # passing these to mpmath functions or returning them in final results.
  30 INF = 1e1000
  31 MINUS_INF = -1e1000
  32
  33 # ~= 100 digits. Real men set this to INF.
  34 DEFAULT_MAXPREC = 333
  35
  36 class PrecisionExhausted(ArithmeticError):
  37     pass
  38
  39 #----------------------------------------------------------------------------#
  40 #                                                                            #
  41 #              Helper functions for arithmetic and complex parts             #
  42 #                                                                            #
  43 #----------------------------------------------------------------------------#
  44
  45 """
  46 An mpf value tuple is a tuple of integers (sign, man, exp, bc)
  47 representing a floating-point numbers.
  48
  49 A temporary result is a tuple (re, im, re_acc, im_acc) where
  50 re and im are nonzero mpf value tuples representing approximate
  51 numbers, or None to denote exact zeros.
  52
  53 re_acc, im_acc are integers denoting log2(e) where e is the estimated
  54 relative accuracy of the respective complex part, but may be anything
  55 if the corresponding complex part is None.
  56
  57 """
  58
  59 def fastlog(x):
  60     """Fast approximation of log2(x) for an mpf value tuple x."""
  61     if not x or x == fzero:
  62         return MINUS_INF
  63     # log2(x) ~= exponent + width of mantissa
  64     # Note: this actually gives ceil(log2(x)), which is a useful
  65     # feature for interval arithmetic.
  66     return x[2] + x[3]
  67
  68 def complex_accuracy(result):
  69     """
  70     Returns relative accuracy of a complex number with given accuracies
  71     for the real and imaginary parts. The relative accuracy is defined
  72     in the complex norm sense as ||z|+|error|| / |z| where error
  73     is equal to (real absolute error) + (imag absolute error)*i.
  74
  75     The full expression for the (logarithmic) error can be approximated
  76     easily by using the max norm to approximate the complex norm.
  77
  78     In the worst case (re and im equal), this is wrong by a factor
  79     sqrt(2), or by log2(sqrt(2)) = 0.5 bit.
  80     """
  81     re, im, re_acc, im_acc = result
  82     if not im:
  83         if not re:
  84             return INF
  85         return re_acc
  86     if not re:
  87         return im_acc
  88     re_size = fastlog(re)
  89     im_size = fastlog(im)
  90     absolute_error = max(re_size-re_acc, im_size-im_acc)
  91     relative_error = absolute_error - max(re_size, im_size)
  92     return -relative_error
  93
  94 def get_abs(expr, prec, options):
  95     re, im, re_acc, im_acc = evalf(expr, prec+2, options)
  96     if not re:
  97         re, re_acc, im, im_acc = im, im_acc, re, re_acc
  98     if im:
  99         return libmpc.mpc_abs((re, im), prec), None, re_acc, None
 100     else:
 101         return fabs(re), None, re_acc, None
 102
 103 def get_complex_part(expr, no, prec, options):
 104     """no = 0 for real part, no = 1 for imaginary part"""
 105     workprec = prec
 106     i = 0
 107     while 1:
 108         res = evalf(expr, workprec, options)
 109         value, accuracy = res[no::2]
 110         if (not value) or accuracy >= prec:
 111             return value, None, accuracy, None
 112         workprec += max(30, 2**i)
 113         i += 1
 114
 115 def evalf_abs(expr, prec, options):
 116     return get_abs(expr.args[0], prec, options)
 117
 118 def evalf_re(expr, prec, options):
 119     return get_complex_part(expr.args[0], 0, prec, options)
 120
 121 def evalf_im(expr, prec, options):
 122     return get_complex_part(expr.args[0], 1, prec, options)
 123
 124 def finalize_complex(re, im, prec):
 125     assert re and im
 126     if re == fzero and im == fzero:
 127         raise ValueError("got complex zero with unknown accuracy")
 128     size_re = fastlog(re)
 129     size_im = fastlog(im)
 130     # Convert fzeros to scaled zeros
 131     if re == fzero:
 132         re = fshift(fone, size_im-prec)
 133         size_re = fastlog(re)
 134     elif im == fzero:
 135         im = fshift(fone, size_re-prec)
 136         size_im = fastlog(im)
 137     if size_re > size_im:
 138         re_acc = prec
 139         im_acc = prec + min(-(size_re - size_im), 0)
 140     else:
 141         im_acc = prec
 142         re_acc = prec + min(-(size_im - size_re), 0)
 143     return re, im, re_acc, im_acc
 144
 145 def chop_parts(value, prec):
 146     """
 147     Chop off tiny real or complex parts.
 148     """
 149     re, im, re_acc, im_acc = value
 150     # Method 1: chop based on absolute value
 151     if re and (fastlog(re) < -prec+4):
 152         re, re_acc = None, None
 153     if im and (fastlog(im) < -prec+4):
 154         im, im_acc = None, None
 155     # Method 2: chop if inaccurate and relatively small
 156     if re and im:
 157         delta = fastlog(re) - fastlog(im)
 158         if re_acc < 2 and (delta - re_acc <= -prec+4):
 159             re, re_acc = None, None
 160         if im_acc < 2 and (delta - im_acc >= prec-4):
 161             im, im_acc = None, None
 162     return re, im, re_acc, im_acc
 163
 164 def check_target(expr, result, prec):
 165     a = complex_accuracy(result)
 166     if a < prec:
 167         raise PrecisionExhausted("Failed to distinguish the expression: \n\n%s\n\n"
 168             "from zero. Try simplifying the input, using chop=True, or providing "
 169             "a higher maxprec for evalf" % (expr))
 170
 171 def get_integer_part(expr, no, options, return_ints=False):
 172     """
 173     With no = 1, computes ceiling(expr)
 174     With no = -1, computes floor(expr)
 175
 176     Note: this function either gives the exact result or signals failure.
 177     """
 178
 179     # The expression is likely less than 2^30 or so
 180     assumed_size = 30
 181     ire, iim, ire_acc, iim_acc = evalf(expr, assumed_size, options)
 182
 183     # We now know the size, so we can calculate how much extra precision
 184     # (if any) is needed to get within the nearest integer
 185     if ire and iim:
 186         gap = max(fastlog(ire)-ire_acc, fastlog(iim)-iim_acc)
 187     elif ire:
 188         gap = fastlog(ire)-ire_acc
 189     elif iim:
 190         gap = fastlog(iim)-iim_acc
 191     else:
 192         # ... or maybe the expression was exactly zero
 193         return None, None, None, None
 194
 195     margin = 10
 196
 197     if gap >= -margin:
 198         ire, iim, ire_acc, iim_acc = evalf(expr, margin+assumed_size+gap, options)
 199
 200     # We can now easily find the nearest integer, but to find floor/ceil, we
 201     # must also calculate whether the difference to the nearest integer is
 202     # positive or negative (which may fail if very close)
 203     def calc_part(expr, nexpr):
 204         nint = int(to_int(nexpr, round_nearest))
 205         expr = C.Add(expr, -nint, evaluate=False)
 206         x, _, x_acc, _ = evalf(expr, 10, options)
 207         check_target(expr, (x, None, x_acc, None), 3)
 208         nint += int(no*(fcmp(x or fzero, fzero) == no))
 209         nint = from_int(nint)
 210         return nint, fastlog(nint) + 10
 211
 212     re, im, re_acc, im_acc = None, None, None, None
 213
 214     if ire:
 215         re, re_acc = calc_part(C.re(expr, evaluate=False), ire)
 216     if iim:
 217         im, im_acc = calc_part(C.im(expr, evaluate=False), iim)
 218
 219     if return_ints:
 220         return int(to_int(re or fzero)), int(to_int(im or fzero))
 221     return re, im, re_acc, im_acc
 222
 223 def evalf_ceiling(expr, prec, options):
 224     return get_integer_part(expr.args[0], 1, options)
 225
 226 def evalf_floor(expr, prec, options):
 227     return get_integer_part(expr.args[0], -1, options)
 228
 229 #----------------------------------------------------------------------------#
 230 #                                                                            #
 231 #                            Arithmetic operations                           #
 232 #                                                                            #
 233 #----------------------------------------------------------------------------#
 234
 235 def add_terms(terms, prec, target_prec):
 236     """
 237     Helper for evalf_add. Adds a list of (mpfval, accuracy) terms.
 238     """
 239     if len(terms) == 1:
 240         if not terms[0]:
 241             # XXX: this is supposed to represent a scaled zero
 242             return fshift(fone, target_prec), -1
 243         return terms[0]
 244     sum_man, sum_exp, absolute_error = 0, 0, MINUS_INF
 245     for x, accuracy in terms:
 246         if not x:
 247             continue
 248         sign, man, exp, bc = x
 249         if sign:
 250             man = -man
 251         absolute_error = max(absolute_error, bc+exp-accuracy)
 252         delta = exp - sum_exp
 253         if exp >= sum_exp:
 254             if delta > 4*prec:
 255                 sum_man = man
 256                 sum_exp = exp
 257             else:
 258                 sum_man += man << delta
 259         else:
 260             if (-delta) > 4*prec:
 261                 pass
 262             else:
 263                 sum_man = (sum_man << (-delta)) + man
 264                 sum_exp = exp
 265     if absolute_error == MINUS_INF:
 266         return None, None
 267     if not sum_man:
 268         # XXX: this is supposed to represent a scaled zero
 269         return fshift(fone, absolute_error), -1
 270     if sum_man < 0:
 271         sum_sign = 1
 272         sum_man = -sum_man
 273     else:
 274         sum_sign = 0
 275     sum_bc = bitcount(sum_man)
 276     sum_accuracy = sum_exp + sum_bc - absolute_error
 277     r = normalize(sum_sign, sum_man, sum_exp, sum_bc, target_prec,
 278         round_nearest), sum_accuracy
 279     #print "returning", to_str(r[0],50), r[1]
 280     return r
 281
 282 def evalf_add(v, prec, options):
 283     args = v.args
 284     target_prec = prec
 285     i = 0
 286
 287     oldmaxprec = options.get('maxprec', DEFAULT_MAXPREC)
 288     options['maxprec'] = min(oldmaxprec, 2*prec)
 289
 290     try:
 291         while 1:
 292             terms = [evalf(arg, prec+10, options) for arg in args]
 293             re, re_accuracy = add_terms([(a[0],a[2]) for a in terms if a[0]], prec, target_prec)
 294             im, im_accuracy = add_terms([(a[1],a[3]) for a in terms if a[1]], prec, target_prec)
 295             accuracy = complex_accuracy((re, im, re_accuracy, im_accuracy))
 296             if accuracy >= target_prec:
 297                 if options.get('verbose'):
 298                     print "ADD: wanted", target_prec, "accurate bits, got", re_accuracy, im_accuracy
 299                 return re, im, re_accuracy, im_accuracy
 300             else:
 301                 diff = target_prec - accuracy
 302                 if (prec-target_prec) > options.get('maxprec', DEFAULT_MAXPREC):
 303                     return re, im, re_accuracy, im_accuracy
 304
 305                 prec = prec + max(10+2**i, diff)
 306                 options['maxprec'] = min(oldmaxprec, 2*prec)
 307                 if options.get('verbose'):
 308                     print "ADD: restarting with prec", prec
 309             i += 1
 310     finally:
 311         options['maxprec'] = oldmaxprec
 312
 313 # Helper for complex multiplication
 314 # XXX: should be able to multiply directly, and use complex_accuracy
 315 # to obtain the final accuracy
 316 def cmul((a, aacc), (b, bacc), (c, cacc), (d, dacc), prec, target_prec):
 317     A, Aacc = fmul(a,c,prec), min(aacc, cacc)
 318     B, Bacc = fmul(fneg(b),d,prec), min(bacc, dacc)
 319     C, Cacc = fmul(a,d,prec), min(aacc, dacc)
 320     D, Dacc = fmul(b,c,prec), min(bacc, cacc)
 321     re, re_accuracy = add_terms([(A, Aacc), (B, Bacc)], prec, target_prec)
 322     im, im_accuracy = add_terms([(C, Cacc), (D, Cacc)], prec, target_prec)
 323     return re, im, re_accuracy, im_accuracy
 324
 325 def evalf_mul(v, prec, options):
 326     args = v.args
 327     # With guard digits, multiplication in the real case does not destroy
 328     # accuracy. This is also true in the complex case when considering the
 329     # total accuracy; however accuracy for the real or imaginary parts
 330     # separately may be lower.
 331     acc = prec
 332     target_prec = prec
 333     # XXX: big overestimate
 334     prec = prec + len(args) + 5
 335     direction = 0
 336     # Empty product is 1
 337     man, exp, bc = 1, 0, 1
 338     direction = 0
 339     complex_factors = []
 340     # First, we multiply all pure real or pure imaginary numbers.
 341     # direction tells us that the result should be multiplied by
 342     # i**direction
 343     for arg in args:
 344         re, im, a, aim = evalf(arg, prec, options)
 345         if re and im:
 346             complex_factors.append((re, im, a, aim))
 347             continue
 348         elif re:
 349             s, m, e, b = re
 350         elif im:
 351             a = aim
 352             direction += 1
 353             s, m, e, b = im
 354         else:
 355             return None, None, None, None
 356         direction += 2*s
 357         man *= m
 358         exp += e
 359         bc += b
 360         if bc > 3*prec:
 361             man >>= prec
 362             exp += prec
 363         acc = min(acc, a)
 364     sign = (direction & 2) >> 1
 365     v = normalize(sign, man, exp, bitcount(man), prec, round_nearest)
 366     if complex_factors:
 367         # Multiply first complex number by the existing real scalar
 368         re, im, re_acc, im_acc = complex_factors[0]
 369         re = fmul(re, v, prec)
 370         im = fmul(im, v, prec)
 371         re_acc = min(re_acc, acc)
 372         im_acc = min(im_acc, acc)
 373         # Multiply consecutive complex factors
 374         complex_factors = complex_factors[1:]
 375         for wre, wim, wre_acc, wim_acc in complex_factors:
 376             re, im, re_acc, im_acc = cmul((re, re_acc), (im,im_acc),
 377                 (wre,wre_acc), (wim,wim_acc), prec, target_prec)
 378         if options.get('verbose'):
 379             print "MUL: obtained accuracy", re_acc, im_acc, "expected", target_prec
 380         # multiply by i
 381         if direction & 1:
 382             return fneg(im), re, re_acc, im_acc
 383         else:
 384             return re, im, re_acc, im_acc
 385     else:
 386         # multiply by i
 387         if direction & 1:
 388             return None, v, None, acc
 389         else:
 390             return v, None, acc, None
 391
 392 def evalf_pow(v, prec, options):
 393
 394     target_prec = prec
 395     base, exp = v.args
 396
 397     # We handle x**n separately. This has two purposes: 1) it is much
 398     # faster, because we avoid calling evalf on the exponent, and 2) it
 399     # allows better handling of real/imaginary parts that are exactly zero
 400     if exp.is_Integer:
 401         p = exp.p
 402         # Exact
 403         if not p:
 404             return fone, None, prec, None
 405         # Exponentiation by p magnifies relative error by |p|, so the
 406         # base must be evaluated with increased precision if p is large
 407         prec += int(math.log(abs(p),2))
 408         re, im, re_acc, im_acc = evalf(base, prec+5, options)
 409         # Real to integer power
 410         if re and not im:
 411             return fpowi(re, p, target_prec), None, target_prec, None
 412         # (x*I)**n = I**n * x**n
 413         if im and not re:
 414             z = fpowi(im, p, target_prec)
 415             case = p % 4
 416             if case == 0: return z, None, target_prec, None
 417             if case == 1: return None, z, None, target_prec
 418             if case == 2: return fneg(z), None, target_prec, None
 419             if case == 3: return None, fneg(z), None, target_prec
 420         # General complex number to arbitrary integer power
 421         re, im = libmpc.mpc_pow_int((re, im), p, prec)
 422         # Assumes full accuracy in input
 423         return finalize_complex(re, im, target_prec)
 424
 425     # Pure square root
 426     if exp is S.Half:
 427         xre, xim, xre_acc, yim_acc = evalf(base, prec+5, options)
 428         # General complex square root
 429         if xim:
 430             re, im = libmpc.mpc_sqrt((xre or fzero, xim), prec)
 431             return finalize_complex(re, im, prec)
 432         if not xre:
 433             return None, None, None, None
 434         # Square root of a negative real number
 435         if flt(xre, fzero):
 436             return None, fsqrt(fneg(xre), prec), None, prec
 437         # Positive square root
 438         return fsqrt(xre, prec), None, prec, None
 439
 440     # We first evaluate the exponent to find its magnitude
 441     # This determines the working precision that must be used
 442     prec += 10
 443     yre, yim, yre_acc, yim_acc = evalf(exp, prec, options)
 444     ysize = fastlog(yre)
 445     # Restart if too big
 446     # XXX: prec + ysize might exceed maxprec
 447     if ysize > 5:
 448         prec += ysize
 449         yre, yim, yre_acc, yim_acc = evalf(exp, prec, options)
 450
 451     # Pure exponential function; no need to evalf the base
 452     if base is S.Exp1:
 453         if yim:
 454             re, im = libmpc.mpc_exp((yre or fzero, yim), prec)
 455             return finalize_complex(re, im, target_prec)
 456         return fexp(yre, target_prec), None, target_prec, None
 457
 458     xre, xim, xre_acc, yim_acc = evalf(base, prec+5, options)
 459
 460     # (real ** complex) or (complex ** complex)
 461     if yim:
 462         re, im = libmpc.mpc_pow((xre or fzero, xim or fzero), (yre or fzero, yim),
 463             target_prec)
 464         return finalize_complex(re, im, target_prec)
 465     # complex ** real
 466     if xim:
 467         re, im = libmpc.mpc_pow_mpf((xre or fzero, xim), yre, target_prec)
 468         return finalize_complex(re, im, target_prec)
 469     # negative ** real
 470     elif flt(xre, fzero):
 471         re, im = libmpc.mpc_pow_mpf((xre, fzero), yre, target_prec)
 472         return finalize_complex(re, im, target_prec)
 473     # positive ** real
 474     else:
 475         return fpow(xre, yre, target_prec), None, target_prec, None
 476
 477
 478
 479
 480 #----------------------------------------------------------------------------#
 481 #                                                                            #
 482 #                            Special functions                               #
 483 #                                                                            #
 484 #----------------------------------------------------------------------------#
 485
 486 def evalf_trig(v, prec, options):
 487     """
 488     This function handles sin and cos of real arguments.
 489
 490     TODO: should also handle tan and complex arguments.
 491     """
 492     if v.func is C.cos:
 493         func = fcos
 494     elif v.func is C.sin:
 495         func = fsin
 496     else:
 497         raise NotImplementedError
 498     arg = v.args[0]
 499     # 20 extra bits is possibly overkill. It does make the need
 500     # to restart very unlikely
 501     xprec = prec + 20
 502     re, im, re_accuracy, im_accuracy = evalf(arg, xprec, options)
 503     if im:
 504         raise NotImplementedError
 505     if not re:
 506         if v.func is C.cos:
 507             return fone, None, prec, None
 508         elif v.func is C.sin:
 509             return None, None, None, None
 510         else:
 511             raise NotImplementedError
 512     # For trigonometric functions, we are interested in the
 513     # fixed-point (absolute) accuracy of the argument.
 514     xsize = fastlog(re)
 515     # Magnitude <= 1.0. OK to compute directly, because there is no
 516     # danger of hitting the first root of cos (with sin, magnitude
 517     # <= 2.0 would actually be ok)
 518     if xsize < 1:
 519         return func(re, prec, round_nearest), None, prec, None
 520     # Very large
 521     if xsize >= 10:
 522         xprec = prec + xsize
 523         re, im, re_accuracy, im_accuracy = evalf(arg, xprec, options)
 524     # Need to repeat in case the argument is very close to a
 525     # multiple of pi (or pi/2), hitting close to a root
 526     while 1:
 527         y = func(re, prec, round_nearest)
 528         ysize = fastlog(y)
 529         gap = -ysize
 530         accuracy = (xprec - xsize) - gap
 531         if accuracy < prec:
 532             if options.get('verbose'):
 533                 print "SIN/COS", accuracy, "wanted", prec, "gap", gap
 534                 print to_str(y,10)
 535             if xprec > options.get('maxprec', DEFAULT_MAXPREC):
 536                 return y, None, accuracy, None
 537             xprec += gap
 538             re, im, re_accuracy, im_accuracy = evalf(arg, xprec, options)
 539             continue
 540         else:
 541             return y, None, prec, None
 542
 543 def evalf_log(expr, prec, options):
 544     arg = expr.args[0]
 545     workprec = prec+10
 546     xre, xim, xacc, _ = evalf(arg, workprec, options)
 547
 548     if xim:
 549         # XXX: use get_abs etc instead
 550         re = evalf_log(C.log(C.abs(arg, evaluate=False), evaluate=False), prec, options)
 551         im = fatan2(xim, xre or fzero, prec)
 552         return re[0], im, re[2], prec
 553
 554     imaginary_term = (fcmp(xre, fzero) < 0)
 555
 556     re = flog(fabs(xre), prec, round_nearest)
 557     size = fastlog(re)
 558     if prec - size > workprec:
 559         # We actually need to compute 1+x accurately, not x
 560         arg = C.Add(S.NegativeOne,arg,evaluate=False)
 561         xre, xim, xre_acc, xim_acc = evalf_add(arg, prec, options)
 562         prec2 = workprec - fastlog(xre)
 563         re = flog(fadd(xre, fone, prec2), prec, round_nearest)
 564
 565     re_acc = prec
 566
 567     if imaginary_term:
 568         return re, fpi(prec), re_acc, prec
 569     else:
 570         return re, None, re_acc, None
 571
 572 def evalf_atan(v, prec, options):
 573     arg = v.args[0]
 574     xre, xim, reacc, imacc = evalf(arg, prec+5, options)
 575     if xim:
 576         raise NotImplementedError
 577     return fatan(xre, prec, round_nearest), None, prec, None
 578
 579
 580 #----------------------------------------------------------------------------#
 581 #                                                                            #
 582 #                            High-level operations                           #
 583 #                                                                            #
 584 #----------------------------------------------------------------------------#
 585
 586 def as_mpmath(x, prec, options):
 587     x = sympify(x)
 588     if isinstance(x, C.Zero):
 589         return mpf(0)
 590     if isinstance(x, C.Infinity):
 591         return mpf('inf')
 592     if isinstance(x, C.NegativeInfinity):
 593         return mpf('-inf')
 594     # XXX
 595     re, im, _, _ = evalf(x, prec, options)
 596     if im:
 597         return mpc(re or fzero, im)
 598     return mpf(re)
 599
 600 def do_integral(expr, prec, options):
 601     func = expr.args[0]
 602     x, (xlow, xhigh) = expr.args[1][0]
 603     orig = mp.prec
 604
 605     oldmaxprec = options.get('maxprec', DEFAULT_MAXPREC)
 606     options['maxprec'] = min(oldmaxprec, 2*prec)
 607
 608     try:
 609         mp.prec = prec+5
 610         xlow = as_mpmath(xlow, prec+15, options)
 611         xhigh = as_mpmath(xhigh, prec+15, options)
 612
 613         # Integration is like summation, and we can phone home from
 614         # the integrand function to update accuracy summation style
 615         # Note that this accuracy is inaccurate, since it fails
 616         # to account for the variable quadrature weights,
 617         # but it is better than nothing
 618
 619         have_part = [False, False]
 620         max_real_term = [MINUS_INF]
 621         max_imag_term = [MINUS_INF]
 622
 623         def f(t):
 624             re, im, re_acc, im_acc = evalf(func, prec+15, {'subs':{x:t}})
 625
 626             have_part[0] = re or have_part[0]
 627             have_part[1] = im or have_part[1]
 628
 629             max_real_term[0] = max(max_real_term[0], fastlog(re))
 630             max_imag_term[0] = max(max_imag_term[0], fastlog(im))
 631
 632             if im:
 633                 return mpc(re or fzero, im)
 634             return mpf(re or fzero)
 635
 636         if options.get('quad') == 'osc':
 637             A = C.Wild('A', exclude=[x])
 638             B = C.Wild('B', exclude=[x])
 639             D = C.Wild('D')
 640             m = func.match(C.cos(A*x+B)*D)
 641             if not m:
 642                 m = func.match(C.sin(A*x+B)*D)
 643             if not m:
 644                 raise ValueError("An integrand of the form sin(A*x+B)*f(x) "
 645                   "or cos(A*x+B)*f(x) is required for oscillatory quadrature")
 646             period = as_mpmath(2*S.Pi/m[A], prec+15, options)
 647             result = quadosc(f, xlow, xhigh, period=period)
 648             # XXX: quadosc does not do error detection yet
 649             quadrature_error = MINUS_INF
 650         else:
 651             result, quadrature_error = quadts(f, xlow, xhigh, error=1)
 652             quadrature_error = fastlog(quadrature_error._mpf_)
 653
 654     finally:
 655         options['maxprec'] = oldmaxprec
 656         mp.prec = orig
 657
 658     if have_part[0]:
 659         re = result.real._mpf_
 660         if re == fzero:
 661             re = fshift(fone, min(-prec,-max_real_term[0],-quadrature_error))
 662             re_acc = -1
 663         else:
 664             re_acc = -max(max_real_term[0]-fastlog(re)-prec, quadrature_error)
 665     else:
 666         re, re_acc = None, None
 667
 668     if have_part[1]:
 669         im = result.imag._mpf_
 670         if im == fzero:
 671             im = fshift(fone, min(-prec,-max_imag_term[0],-quadrature_error))
 672             im_acc = -1
 673         else:
 674             im_acc = -max(max_imag_term[0]-fastlog(im)-prec, quadrature_error)
 675     else:
 676         im, im_acc = None, None
 677
 678     result = re, im, re_acc, im_acc
 679     return result
 680
 681 def evalf_integral(expr, prec, options):
 682     workprec = prec
 683     i = 0
 684     maxprec = options.get('maxprec', INF)
 685     while 1:
 686         result = do_integral(expr, workprec, options)
 687         accuracy = complex_accuracy(result)
 688         if accuracy >= prec or workprec >= maxprec:
 689             return result
 690         workprec += prec - max(-2**i, accuracy)
 691         i += 1
 692
 693 def check_convergence(numer, denom, n):
 694     """
 695     Returns (h, g, p) where
 696     -- h is:
 697         > 0 for convergence of rate 1/factorial(n)**h
 698         < 0 for divergence of rate factorial(n)**(-h)
 699         = 0 for geometric or polynomial convergence or divergence
 700
 701     -- abs(g) is:
 702         > 1 for geometric convergence of rate 1/h**n
 703         < 1 for geometric divergence of rate h**n
 704         = 1 for polynomial convergence or divergence
 705
 706         (g < 0 indicates an alternating series)
 707
 708     -- p is:
 709         > 1 for polynomial convergence of rate 1/n**h
 710         <= 1 for polynomial divergence of rate n**(-h)
 711
 712     """
 713     npol = C.Poly(numer, n)
 714     dpol = C.Poly(denom, n)
 715     p = npol.degree
 716     q = dpol.degree
 717     rate = q - p
 718     if rate:
 719         return rate, None, None
 720     constant = dpol.lead_term[0] / npol.lead_term[0]
 721     if abs(constant) != 1:
 722         return rate, constant, None
 723     if npol.degree == dpol.degree == 0:
 724         return rate, constant, 0
 725     pc = list(npol.iter_all_terms())[1][0]
 726     qc = list(dpol.iter_all_terms())[1][0]
 727     return rate, constant, qc-pc
 728
 729 def hypsum(expr, n, start, prec):
 730     """
 731     Sum a rapidly convergent infinite hypergeometric series with
 732     given general term, e.g. e = hypsum(1/factorial(n), n). The
 733     quotient between successive terms must be a quotient of integer
 734     polynomials.
 735     """
 736     from sympy import hypersimp, lambdify
 737
 738     if start:
 739         expr = expr.subs(n, n+start)
 740     hs = hypersimp(expr, n)
 741     if hs is None:
 742         raise NotImplementedError("a hypergeometric series is required")
 743     num, den = hs.as_numer_denom()
 744     func1 = lambdify(n, num)
 745     func2 = lambdify(n, den)
 746
 747     h, g, p = check_convergence(num, den, n)
 748
 749     if h < 0:
 750         raise ValueError("Sum diverges like (n!)^%i" % (-h))
 751
 752     # Direct summation if geometric or faster
 753     if h > 0 or (h == 0 and abs(g) > 1):
 754         one = MP_BASE(1) << prec
 755         term = expr.subs(n, 0)
 756         term = (MP_BASE(term.p) << prec) // term.q
 757         s = term
 758         k = 1
 759         while abs(term) > 5:
 760             term *= MP_BASE(func1(k-1))
 761             term //= MP_BASE(func2(k-1))
 762             s += term
 763             k += 1
 764         return from_man_exp(s, -prec)
 765     else:
 766         alt = g < 0
 767         if abs(g) < 1:
 768             raise ValueError("Sum diverges like (%i)^n" % abs(1/g))
 769         if p < 1 or (p == 1 and not alt):
 770             raise ValueError("Sum diverges like n^%i" % (-p))
 771         # We have polynomial convergence:
 772         # Use Shanks extrapolation for alternating series,
 773         # Richardson extrapolation for nonalternating series
 774         if alt:
 775             # XXX: better parameters for Shanks transformation
 776             # This tends to get bad somewhere > 50 digits
 777             N = 5 + int(prec*0.36)
 778             M = 2 + N//3
 779             NTERMS = M + N + 2
 780         else:
 781             N = 3 + int(prec*0.15)
 782             M = 2*N
 783             NTERMS = M + N + 2
 784         # Need to use at least double precision because a lot of cancellation
 785         # might occur in the extrapolation process
 786         prec2 = 2*prec
 787         one = MP_BASE(1) << prec2
 788         term = expr.subs(n, 0)
 789         term = (MP_BASE(term.p) << prec2) // term.q
 790         s = term
 791         table = [make_mpf(from_man_exp(s, -prec2))]
 792         for k in xrange(1, NTERMS):
 793             term *= MP_BASE(func1(k-1))
 794             term //= MP_BASE(func2(k-1))
 795             s += term
 796             table.append(make_mpf(from_man_exp(s, -prec2)))
 797             k += 1
 798         orig = mp.prec
 799         try:
 800             mp.prec = prec
 801             if alt:
 802                 v = shanks_extrapolation(table, N, M)
 803             else:
 804                 v = richardson_extrapolation(table, N, M)
 805         finally:
 806             mp.prec = orig
 807         return v._mpf_
 808
 809 def evalf_sum(expr, prec, options):
 810     func = expr.function
 811     limits = expr.limits
 812     if len(limits) != 1 or not isinstance(limits[0], tuple) or \
 813         len(limits[0]) != 3:
 814         raise NotImplementedError
 815     prec2 = prec+10
 816     try:
 817         n, a, b = limits[0]
 818         if b != S.Infinity or a != int(a):
 819             raise NotImplementedError
 820         # Use fast hypergeometric summation if possible
 821         v = hypsum(func, n, int(a), prec2)
 822         delta = prec - fastlog(v)
 823         if fastlog(v) < -10:
 824             v = hypsum(func, n, int(a), delta)
 825         return v, None, min(prec, delta), None
 826     except NotImplementedError:
 827         # Euler-Maclaurin summation for general series
 828         eps = C.Real(2.0)**(-prec)
 829         for i in range(1, 5):
 830             m = n = 2**i * prec
 831             s, err = expr.euler_maclaurin(m=m, n=n, eps=eps, \
 832                 eval_integral=False)
 833             err = err.evalf()
 834             if err <= eps:
 835                 break
 836         err = fastlog(evalf(abs(err), 20, options)[0])
 837         re, im, re_acc, im_acc = evalf(s, prec2, options)
 838         re_acc = max(re_acc, -err)
 839         im_acc = max(im_acc, -err)
 840         return re, im, re_acc, im_acc
 841
 842
 843 #----------------------------------------------------------------------------#
 844 #                                                                            #
 845 #                            Symbolic interface                              #
 846 #                                                                            #
 847 #----------------------------------------------------------------------------#
 848
 849 def evalf_symbol(x, prec, options):
 850     val = options['subs'][x]
 851     if isinstance(val, mpf):
 852         if not val:
 853             return None, None, None, None
 854         return val._mpf_, None, prec, None
 855     else:
 856         if not '_cache' in options:
 857             options['_cache'] = {}
 858         cache = options['_cache']
 859         cached, cached_prec = cache.get(x.name, (None, MINUS_INF))
 860         if cached_prec >= prec:
 861             return cached
 862         v = evalf(sympify(val), prec, options)
 863         cache[x.name] = (v, prec)
 864         return v
 865
 866 evalf_table = None
 867
 868 def _create_evalf_table():
 869     global evalf_table
 870     evalf_table = {
 871     C.Symbol : evalf_symbol,
 872     C.Dummy : evalf_symbol,
 873     C.Real : lambda x, prec, options: (x._mpf_, None, prec, None),
 874     C.Rational : lambda x, prec, options: (from_rational(x.p, x.q, prec), None, prec, None),
 875     C.Integer : lambda x, prec, options: (from_int(x.p, prec), None, prec, None),
 876     C.Zero : lambda x, prec, options: (None, None, prec, None),
 877     C.One : lambda x, prec, options: (fone, None, prec, None),
 878     C.Half : lambda x, prec, options: (fhalf, None, prec, None),
 879     C.Pi : lambda x, prec, options: (fpi(prec), None, prec, None),
 880     C.Exp1 : lambda x, prec, options: (fe(prec), None, prec, None),
 881     C.ImaginaryUnit : lambda x, prec, options: (None, fone, None, prec),
 882     C.NegativeOne : lambda x, prec, options: (fnone, None, prec, None),
 883
 884     C.exp : lambda x, prec, options: evalf_pow(C.Pow(S.Exp1, x.args[0],
 885         evaluate=False), prec, options),
 886
 887     C.cos : evalf_trig,
 888     C.sin : evalf_trig,
 889
 890     C.Add : evalf_add,
 891     C.Mul : evalf_mul,
 892     C.Pow : evalf_pow,
 893
 894     C.log : evalf_log,
 895     C.atan : evalf_atan,
 896     C.abs : evalf_abs,
 897
 898     C.re : evalf_re,
 899     C.im : evalf_im,
 900     C.floor : evalf_floor,
 901     C.ceiling : evalf_ceiling,
 902
 903     C.Integral : evalf_integral,
 904     C.Sum : evalf_sum,
 905     }
 906
 907 def evalf(x, prec, options):
 908     try:
 909         r = evalf_table[x.func](x, prec, options)
 910     except KeyError:
 911         #r = finalize_complex(x._eval_evalf(prec)._mpf_, fzero, prec)
 912         try:
 913             # Fall back to ordinary evalf if possible
 914             if 'subs' in options:
 915                 x = x.subs(options['subs'])
 916             r = x._eval_evalf(prec)._mpf_, None, prec, None
 917         except AttributeError:
 918             raise NotImplementedError
 919     if options.get("verbose"):
 920         print "### input", x
 921         print "### output", to_str(r[0] or fzero, 50)
 922         #print "### raw", r[0], r[2]
 923         print
 924     if options.get("chop"):
 925         r = chop_parts(r, prec)
 926     if options.get("strict"):
 927         check_target(x, r, prec)
 928     return r
 929
 930 def Basic_evalf(x, n=15, **options):
 931     """
 932     Evaluate the given formula to an accuracy of n digits.
 933     Optional keyword arguments:
 934
 935         subs=<dict>
 936             Substitute numerical values for symbols, e.g.
 937             subs={x:3, y:1+pi}.
 938
 939         maxprec=N
 940             Allow a maximum temporary working precision of N digits
 941             (default=100)
 942
 943         chop=<bool>
 944             Replace tiny real or imaginary parts in subresults
 945             by exact zeros (default=False)
 946
 947         strict=<bool>
 948             Raise PrecisionExhausted if any subresult fails to evaluate
 949             to full accuracy, given the available maxprec
 950             (default=False)
 951
 952         quad=<str>
 953             Choose algorithm for numerical quadrature. By default,
 954             tanh-sinh quadrature is used. For oscillatory
 955             integrals on an infinite interval, try quad='osc'.
 956
 957         verbose=<bool>
 958             Print debug information (default=False)
 959
 960     """
 961     if not evalf_table:
 962         _create_evalf_table()
 963     prec = dps_to_prec(n)
 964     if 'maxprec' in options:
 965         options['maxprec'] = int(options['maxprec']*LG10)
 966     else:
 967         options['maxprec'] = max(prec, DEFAULT_MAXPREC)
 968     try:
 969         result = evalf(x, prec+4, options)
 970     except NotImplementedError:
 971         # Fall back to the ordinary evalf
 972         v = x._eval_evalf(prec)
 973         if v is None:
 974             return x
 975         try:
 976             # If the result is numerical, normalize it
 977             result = evalf(v, prec, options)
 978         except:
 979             # Probably contains symbols or unknown functions
 980             return v
 981     re, im, re_acc, im_acc = result
 982     if re:
 983         p = max(min(prec, re_acc), 1)
 984         #re = fpos(re, p, round_nearest)
 985         re = C.Real._new(re, p)
 986     else:
 987         re = S.Zero
 988     if im:
 989         p = max(min(prec, im_acc), 1)
 990         #im = fpos(im, p, round_nearest)
 991         im = C.Real._new(im, p)
 992         return re + im*S.ImaginaryUnit
 993     else:
 994         return re
 995
 996 Basic.evalf = Basic.n = Basic_evalf
 997
 998 def N(x, n=15, **options):
 999     """
1000     Calls x.evalf(n, **options).
1001
1002     Both .evalf() and N() are equivalent, use the one that you like better.
1003
1004     Example:
1005     >>> from sympy import Sum, Symbol, oo
1006     >>> k = Symbol("k")
1007     >>> Sum(1/k**k, (k, 1, oo))
1008     Sum(k**(-k), (k, 1, oo))
1009     >>> N(Sum(1/k**k, (k, 1, oo)), 4)
1010     1.291
1011
1012     """
1013     return sympify(x).evalf(n, **options)