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