1""" 2Adaptive numerical evaluation of SymPy expressions, using mpmath 3for mathematical functions. 4""" 5 6from typing import Tuple 7 8import math 9 10import mpmath.libmp as libmp 11from mpmath import ( 12 make_mpc, make_mpf, mp, mpc, mpf, nsum, quadts, quadosc, workprec) 13from mpmath import inf as mpmath_inf 14from mpmath.libmp import (from_int, from_man_exp, from_rational, fhalf, 15 fnan, fnone, fone, fzero, mpf_abs, mpf_add, 16 mpf_atan, mpf_atan2, mpf_cmp, mpf_cos, mpf_e, mpf_exp, mpf_log, mpf_lt, 17 mpf_mul, mpf_neg, mpf_pi, mpf_pow, mpf_pow_int, mpf_shift, mpf_sin, 18 mpf_sqrt, normalize, round_nearest, to_int, to_str) 19from mpmath.libmp import bitcount as mpmath_bitcount 20from mpmath.libmp.backend import MPZ 21from mpmath.libmp.libmpc import _infs_nan 22from mpmath.libmp.libmpf import dps_to_prec, prec_to_dps 23from mpmath.libmp.gammazeta import mpf_bernoulli 24 25from .compatibility import SYMPY_INTS 26from .sympify import sympify 27from .singleton import S 28 29from sympy.utilities.iterables import is_sequence 30 31LG10 = math.log(10, 2) 32rnd = round_nearest 33 34 35def bitcount(n): 36 """Return smallest integer, b, such that |n|/2**b < 1. 37 """ 38 return mpmath_bitcount(abs(int(n))) 39 40# Used in a few places as placeholder values to denote exponents and 41# precision levels, e.g. of exact numbers. Must be careful to avoid 42# passing these to mpmath functions or returning them in final results. 43INF = float(mpmath_inf) 44MINUS_INF = float(-mpmath_inf) 45 46# ~= 100 digits. Real men set this to INF. 47DEFAULT_MAXPREC = 333 48 49 50class PrecisionExhausted(ArithmeticError): 51 pass 52 53#----------------------------------------------------------------------------# 54# # 55# Helper functions for arithmetic and complex parts # 56# # 57#----------------------------------------------------------------------------# 58 59""" 60An mpf value tuple is a tuple of integers (sign, man, exp, bc) 61representing a floating-point number: [1, -1][sign]*man*2**exp where 62sign is 0 or 1 and bc should correspond to the number of bits used to 63represent the mantissa (man) in binary notation, e.g. 64 65Explanation 66=========== 67 68>>> from sympy.core.evalf import bitcount 69>>> sign, man, exp, bc = 0, 5, 1, 3 70>>> n = [1, -1][sign]*man*2**exp 71>>> n, bitcount(man) 72(10, 3) 73 74A temporary result is a tuple (re, im, re_acc, im_acc) where 75re and im are nonzero mpf value tuples representing approximate 76numbers, or None to denote exact zeros. 77 78re_acc, im_acc are integers denoting log2(e) where e is the estimated 79relative accuracy of the respective complex part, but may be anything 80if the corresponding complex part is None. 81 82""" 83 84 85def fastlog(x): 86 """Fast approximation of log2(x) for an mpf value tuple x. 87 88 Explanation 89 =========== 90 91 Calculated as exponent + width of mantissa. This is an 92 approximation for two reasons: 1) it gives the ceil(log2(abs(x))) 93 value and 2) it is too high by 1 in the case that x is an exact 94 power of 2. Although this is easy to remedy by testing to see if 95 the odd mpf mantissa is 1 (indicating that one was dealing with 96 an exact power of 2) that would decrease the speed and is not 97 necessary as this is only being used as an approximation for the 98 number of bits in x. The correct return value could be written as 99 "x[2] + (x[3] if x[1] != 1 else 0)". 100 Since mpf tuples always have an odd mantissa, no check is done 101 to see if the mantissa is a multiple of 2 (in which case the 102 result would be too large by 1). 103 104 Examples 105 ======== 106 107 >>> from sympy import log 108 >>> from sympy.core.evalf import fastlog, bitcount 109 >>> s, m, e = 0, 5, 1 110 >>> bc = bitcount(m) 111 >>> n = [1, -1][s]*m*2**e 112 >>> n, (log(n)/log(2)).evalf(2), fastlog((s, m, e, bc)) 113 (10, 3.3, 4) 114 """ 115 116 if not x or x == fzero: 117 return MINUS_INF 118 return x[2] + x[3] 119 120 121def pure_complex(v, or_real=False): 122 """Return a and b if v matches a + I*b where b is not zero and 123 a and b are Numbers, else None. If `or_real` is True then 0 will 124 be returned for `b` if `v` is a real number. 125 126 Examples 127 ======== 128 129 >>> from sympy.core.evalf import pure_complex 130 >>> from sympy import sqrt, I, S 131 >>> a, b, surd = S(2), S(3), sqrt(2) 132 >>> pure_complex(a) 133 >>> pure_complex(a, or_real=True) 134 (2, 0) 135 >>> pure_complex(surd) 136 >>> pure_complex(a + b*I) 137 (2, 3) 138 >>> pure_complex(I) 139 (0, 1) 140 """ 141 h, t = v.as_coeff_Add() 142 if not t: 143 if or_real: 144 return h, t 145 return 146 c, i = t.as_coeff_Mul() 147 if i is S.ImaginaryUnit: 148 return h, c 149 150 151def scaled_zero(mag, sign=1): 152 """Return an mpf representing a power of two with magnitude ``mag`` 153 and -1 for precision. Or, if ``mag`` is a scaled_zero tuple, then just 154 remove the sign from within the list that it was initially wrapped 155 in. 156 157 Examples 158 ======== 159 160 >>> from sympy.core.evalf import scaled_zero 161 >>> from sympy import Float 162 >>> z, p = scaled_zero(100) 163 >>> z, p 164 (([0], 1, 100, 1), -1) 165 >>> ok = scaled_zero(z) 166 >>> ok 167 (0, 1, 100, 1) 168 >>> Float(ok) 169 1.26765060022823e+30 170 >>> Float(ok, p) 171 0.e+30 172 >>> ok, p = scaled_zero(100, -1) 173 >>> Float(scaled_zero(ok), p) 174 -0.e+30 175 """ 176 if type(mag) is tuple and len(mag) == 4 and iszero(mag, scaled=True): 177 return (mag[0][0],) + mag[1:] 178 elif isinstance(mag, SYMPY_INTS): 179 if sign not in [-1, 1]: 180 raise ValueError('sign must be +/-1') 181 rv, p = mpf_shift(fone, mag), -1 182 s = 0 if sign == 1 else 1 183 rv = ([s],) + rv[1:] 184 return rv, p 185 else: 186 raise ValueError('scaled zero expects int or scaled_zero tuple.') 187 188 189def iszero(mpf, scaled=False): 190 if not scaled: 191 return not mpf or not mpf[1] and not mpf[-1] 192 return mpf and type(mpf[0]) is list and mpf[1] == mpf[-1] == 1 193 194 195def complex_accuracy(result): 196 """ 197 Returns relative accuracy of a complex number with given accuracies 198 for the real and imaginary parts. The relative accuracy is defined 199 in the complex norm sense as ||z|+|error|| / |z| where error 200 is equal to (real absolute error) + (imag absolute error)*i. 201 202 The full expression for the (logarithmic) error can be approximated 203 easily by using the max norm to approximate the complex norm. 204 205 In the worst case (re and im equal), this is wrong by a factor 206 sqrt(2), or by log2(sqrt(2)) = 0.5 bit. 207 """ 208 re, im, re_acc, im_acc = result 209 if not im: 210 if not re: 211 return INF 212 return re_acc 213 if not re: 214 return im_acc 215 re_size = fastlog(re) 216 im_size = fastlog(im) 217 absolute_error = max(re_size - re_acc, im_size - im_acc) 218 relative_error = absolute_error - max(re_size, im_size) 219 return -relative_error 220 221 222def get_abs(expr, prec, options): 223 re, im, re_acc, im_acc = evalf(expr, prec + 2, options) 224 225 if not re: 226 re, re_acc, im, im_acc = im, im_acc, re, re_acc 227 if im: 228 if expr.is_number: 229 abs_expr, _, acc, _ = evalf(abs(N(expr, prec + 2)), 230 prec + 2, options) 231 return abs_expr, None, acc, None 232 else: 233 if 'subs' in options: 234 return libmp.mpc_abs((re, im), prec), None, re_acc, None 235 return abs(expr), None, prec, None 236 elif re: 237 return mpf_abs(re), None, re_acc, None 238 else: 239 return None, None, None, None 240 241 242def get_complex_part(expr, no, prec, options): 243 """no = 0 for real part, no = 1 for imaginary part""" 244 workprec = prec 245 i = 0 246 while 1: 247 res = evalf(expr, workprec, options) 248 value, accuracy = res[no::2] 249 # XXX is the last one correct? Consider re((1+I)**2).n() 250 if (not value) or accuracy >= prec or -value[2] > prec: 251 return value, None, accuracy, None 252 workprec += max(30, 2**i) 253 i += 1 254 255 256def evalf_abs(expr, prec, options): 257 return get_abs(expr.args[0], prec, options) 258 259 260def evalf_re(expr, prec, options): 261 return get_complex_part(expr.args[0], 0, prec, options) 262 263 264def evalf_im(expr, prec, options): 265 return get_complex_part(expr.args[0], 1, prec, options) 266 267 268def finalize_complex(re, im, prec): 269 if re == fzero and im == fzero: 270 raise ValueError("got complex zero with unknown accuracy") 271 elif re == fzero: 272 return None, im, None, prec 273 elif im == fzero: 274 return re, None, prec, None 275 276 size_re = fastlog(re) 277 size_im = fastlog(im) 278 if size_re > size_im: 279 re_acc = prec 280 im_acc = prec + min(-(size_re - size_im), 0) 281 else: 282 im_acc = prec 283 re_acc = prec + min(-(size_im - size_re), 0) 284 return re, im, re_acc, im_acc 285 286 287def chop_parts(value, prec): 288 """ 289 Chop off tiny real or complex parts. 290 """ 291 re, im, re_acc, im_acc = value 292 # Method 1: chop based on absolute value 293 if re and re not in _infs_nan and (fastlog(re) < -prec + 4): 294 re, re_acc = None, None 295 if im and im not in _infs_nan and (fastlog(im) < -prec + 4): 296 im, im_acc = None, None 297 # Method 2: chop if inaccurate and relatively small 298 if re and im: 299 delta = fastlog(re) - fastlog(im) 300 if re_acc < 2 and (delta - re_acc <= -prec + 4): 301 re, re_acc = None, None 302 if im_acc < 2 and (delta - im_acc >= prec - 4): 303 im, im_acc = None, None 304 return re, im, re_acc, im_acc 305 306 307def check_target(expr, result, prec): 308 a = complex_accuracy(result) 309 if a < prec: 310 raise PrecisionExhausted("Failed to distinguish the expression: \n\n%s\n\n" 311 "from zero. Try simplifying the input, using chop=True, or providing " 312 "a higher maxn for evalf" % (expr)) 313 314 315def get_integer_part(expr, no, options, return_ints=False): 316 """ 317 With no = 1, computes ceiling(expr) 318 With no = -1, computes floor(expr) 319 320 Note: this function either gives the exact result or signals failure. 321 """ 322 from sympy.functions.elementary.complexes import re, im 323 # The expression is likely less than 2^30 or so 324 assumed_size = 30 325 ire, iim, ire_acc, iim_acc = evalf(expr, assumed_size, options) 326 327 # We now know the size, so we can calculate how much extra precision 328 # (if any) is needed to get within the nearest integer 329 if ire and iim: 330 gap = max(fastlog(ire) - ire_acc, fastlog(iim) - iim_acc) 331 elif ire: 332 gap = fastlog(ire) - ire_acc 333 elif iim: 334 gap = fastlog(iim) - iim_acc 335 else: 336 # ... or maybe the expression was exactly zero 337 if return_ints: 338 return 0, 0 339 else: 340 return None, None, None, None 341 342 margin = 10 343 344 if gap >= -margin: 345 prec = margin + assumed_size + gap 346 ire, iim, ire_acc, iim_acc = evalf( 347 expr, prec, options) 348 else: 349 prec = assumed_size 350 351 # We can now easily find the nearest integer, but to find floor/ceil, we 352 # must also calculate whether the difference to the nearest integer is 353 # positive or negative (which may fail if very close). 354 def calc_part(re_im, nexpr): 355 from sympy.core.add import Add 356 _, _, exponent, _ = nexpr 357 is_int = exponent == 0 358 nint = int(to_int(nexpr, rnd)) 359 if is_int: 360 # make sure that we had enough precision to distinguish 361 # between nint and the re or im part (re_im) of expr that 362 # was passed to calc_part 363 ire, iim, ire_acc, iim_acc = evalf( 364 re_im - nint, 10, options) # don't need much precision 365 assert not iim 366 size = -fastlog(ire) + 2 # -ve b/c ire is less than 1 367 if size > prec: 368 ire, iim, ire_acc, iim_acc = evalf( 369 re_im, size, options) 370 assert not iim 371 nexpr = ire 372 nint = int(to_int(nexpr, rnd)) 373 _, _, new_exp, _ = ire 374 is_int = new_exp == 0 375 if not is_int: 376 # if there are subs and they all contain integer re/im parts 377 # then we can (hopefully) safely substitute them into the 378 # expression 379 s = options.get('subs', False) 380 if s: 381 doit = True 382 from sympy.core.compatibility import as_int 383 # use strict=False with as_int because we take 384 # 2.0 == 2 385 for v in s.values(): 386 try: 387 as_int(v, strict=False) 388 except ValueError: 389 try: 390 [as_int(i, strict=False) for i in v.as_real_imag()] 391 continue 392 except (ValueError, AttributeError): 393 doit = False 394 break 395 if doit: 396 re_im = re_im.subs(s) 397 398 re_im = Add(re_im, -nint, evaluate=False) 399 x, _, x_acc, _ = evalf(re_im, 10, options) 400 try: 401 check_target(re_im, (x, None, x_acc, None), 3) 402 except PrecisionExhausted: 403 if not re_im.equals(0): 404 raise PrecisionExhausted 405 x = fzero 406 nint += int(no*(mpf_cmp(x or fzero, fzero) == no)) 407 nint = from_int(nint) 408 return nint, INF 409 410 re_, im_, re_acc, im_acc = None, None, None, None 411 412 if ire: 413 re_, re_acc = calc_part(re(expr, evaluate=False), ire) 414 if iim: 415 im_, im_acc = calc_part(im(expr, evaluate=False), iim) 416 417 if return_ints: 418 return int(to_int(re_ or fzero)), int(to_int(im_ or fzero)) 419 return re_, im_, re_acc, im_acc 420 421 422def evalf_ceiling(expr, prec, options): 423 return get_integer_part(expr.args[0], 1, options) 424 425 426def evalf_floor(expr, prec, options): 427 return get_integer_part(expr.args[0], -1, options) 428 429#----------------------------------------------------------------------------# 430# # 431# Arithmetic operations # 432# # 433#----------------------------------------------------------------------------# 434 435 436def add_terms(terms, prec, target_prec): 437 """ 438 Helper for evalf_add. Adds a list of (mpfval, accuracy) terms. 439 440 Returns 441 ======= 442 443 - None, None if there are no non-zero terms; 444 - terms[0] if there is only 1 term; 445 - scaled_zero if the sum of the terms produces a zero by cancellation 446 e.g. mpfs representing 1 and -1 would produce a scaled zero which need 447 special handling since they are not actually zero and they are purposely 448 malformed to ensure that they can't be used in anything but accuracy 449 calculations; 450 - a tuple that is scaled to target_prec that corresponds to the 451 sum of the terms. 452 453 The returned mpf tuple will be normalized to target_prec; the input 454 prec is used to define the working precision. 455 456 XXX explain why this is needed and why one can't just loop using mpf_add 457 """ 458 459 terms = [t for t in terms if not iszero(t[0])] 460 if not terms: 461 return None, None 462 elif len(terms) == 1: 463 return terms[0] 464 465 # see if any argument is NaN or oo and thus warrants a special return 466 special = [] 467 from sympy.core.numbers import Float 468 for t in terms: 469 arg = Float._new(t[0], 1) 470 if arg is S.NaN or arg.is_infinite: 471 special.append(arg) 472 if special: 473 from sympy.core.add import Add 474 rv = evalf(Add(*special), prec + 4, {}) 475 return rv[0], rv[2] 476 477 working_prec = 2*prec 478 sum_man, sum_exp, absolute_error = 0, 0, MINUS_INF 479 480 for x, accuracy in terms: 481 sign, man, exp, bc = x 482 if sign: 483 man = -man 484 absolute_error = max(absolute_error, bc + exp - accuracy) 485 delta = exp - sum_exp 486 if exp >= sum_exp: 487 # x much larger than existing sum? 488 # first: quick test 489 if ((delta > working_prec) and 490 ((not sum_man) or 491 delta - bitcount(abs(sum_man)) > working_prec)): 492 sum_man = man 493 sum_exp = exp 494 else: 495 sum_man += (man << delta) 496 else: 497 delta = -delta 498 # x much smaller than existing sum? 499 if delta - bc > working_prec: 500 if not sum_man: 501 sum_man, sum_exp = man, exp 502 else: 503 sum_man = (sum_man << delta) + man 504 sum_exp = exp 505 if not sum_man: 506 return scaled_zero(absolute_error) 507 if sum_man < 0: 508 sum_sign = 1 509 sum_man = -sum_man 510 else: 511 sum_sign = 0 512 sum_bc = bitcount(sum_man) 513 sum_accuracy = sum_exp + sum_bc - absolute_error 514 r = normalize(sum_sign, sum_man, sum_exp, sum_bc, target_prec, 515 rnd), sum_accuracy 516 return r 517 518 519def evalf_add(v, prec, options): 520 res = pure_complex(v) 521 if res: 522 h, c = res 523 re, _, re_acc, _ = evalf(h, prec, options) 524 im, _, im_acc, _ = evalf(c, prec, options) 525 return re, im, re_acc, im_acc 526 527 oldmaxprec = options.get('maxprec', DEFAULT_MAXPREC) 528 529 i = 0 530 target_prec = prec 531 while 1: 532 options['maxprec'] = min(oldmaxprec, 2*prec) 533 534 terms = [evalf(arg, prec + 10, options) for arg in v.args] 535 re, re_acc = add_terms( 536 [a[0::2] for a in terms if a[0]], prec, target_prec) 537 im, im_acc = add_terms( 538 [a[1::2] for a in terms if a[1]], prec, target_prec) 539 acc = complex_accuracy((re, im, re_acc, im_acc)) 540 if acc >= target_prec: 541 if options.get('verbose'): 542 print("ADD: wanted", target_prec, "accurate bits, got", re_acc, im_acc) 543 break 544 else: 545 if (prec - target_prec) > options['maxprec']: 546 break 547 548 prec = prec + max(10 + 2**i, target_prec - acc) 549 i += 1 550 if options.get('verbose'): 551 print("ADD: restarting with prec", prec) 552 553 options['maxprec'] = oldmaxprec 554 if iszero(re, scaled=True): 555 re = scaled_zero(re) 556 if iszero(im, scaled=True): 557 im = scaled_zero(im) 558 return re, im, re_acc, im_acc 559 560 561def evalf_mul(v, prec, options): 562 res = pure_complex(v) 563 if res: 564 # the only pure complex that is a mul is h*I 565 _, h = res 566 im, _, im_acc, _ = evalf(h, prec, options) 567 return None, im, None, im_acc 568 args = list(v.args) 569 570 # see if any argument is NaN or oo and thus warrants a special return 571 special = [] 572 from sympy.core.numbers import Float 573 for arg in args: 574 arg = evalf(arg, prec, options) 575 if arg[0] is None: 576 continue 577 arg = Float._new(arg[0], 1) 578 if arg is S.NaN or arg.is_infinite: 579 special.append(arg) 580 if special: 581 from sympy.core.mul import Mul 582 special = Mul(*special) 583 return evalf(special, prec + 4, {}) 584 585 # With guard digits, multiplication in the real case does not destroy 586 # accuracy. This is also true in the complex case when considering the 587 # total accuracy; however accuracy for the real or imaginary parts 588 # separately may be lower. 589 acc = prec 590 591 # XXX: big overestimate 592 working_prec = prec + len(args) + 5 593 594 # Empty product is 1 595 start = man, exp, bc = MPZ(1), 0, 1 596 597 # First, we multiply all pure real or pure imaginary numbers. 598 # direction tells us that the result should be multiplied by 599 # I**direction; all other numbers get put into complex_factors 600 # to be multiplied out after the first phase. 601 last = len(args) 602 direction = 0 603 args.append(S.One) 604 complex_factors = [] 605 606 for i, arg in enumerate(args): 607 if i != last and pure_complex(arg): 608 args[-1] = (args[-1]*arg).expand() 609 continue 610 elif i == last and arg is S.One: 611 continue 612 re, im, re_acc, im_acc = evalf(arg, working_prec, options) 613 if re and im: 614 complex_factors.append((re, im, re_acc, im_acc)) 615 continue 616 elif re: 617 (s, m, e, b), w_acc = re, re_acc 618 elif im: 619 (s, m, e, b), w_acc = im, im_acc 620 direction += 1 621 else: 622 return None, None, None, None 623 direction += 2*s 624 man *= m 625 exp += e 626 bc += b 627 if bc > 3*working_prec: 628 man >>= working_prec 629 exp += working_prec 630 acc = min(acc, w_acc) 631 sign = (direction & 2) >> 1 632 if not complex_factors: 633 v = normalize(sign, man, exp, bitcount(man), prec, rnd) 634 # multiply by i 635 if direction & 1: 636 return None, v, None, acc 637 else: 638 return v, None, acc, None 639 else: 640 # initialize with the first term 641 if (man, exp, bc) != start: 642 # there was a real part; give it an imaginary part 643 re, im = (sign, man, exp, bitcount(man)), (0, MPZ(0), 0, 0) 644 i0 = 0 645 else: 646 # there is no real part to start (other than the starting 1) 647 wre, wim, wre_acc, wim_acc = complex_factors[0] 648 acc = min(acc, 649 complex_accuracy((wre, wim, wre_acc, wim_acc))) 650 re = wre 651 im = wim 652 i0 = 1 653 654 for wre, wim, wre_acc, wim_acc in complex_factors[i0:]: 655 # acc is the overall accuracy of the product; we aren't 656 # computing exact accuracies of the product. 657 acc = min(acc, 658 complex_accuracy((wre, wim, wre_acc, wim_acc))) 659 660 use_prec = working_prec 661 A = mpf_mul(re, wre, use_prec) 662 B = mpf_mul(mpf_neg(im), wim, use_prec) 663 C = mpf_mul(re, wim, use_prec) 664 D = mpf_mul(im, wre, use_prec) 665 re = mpf_add(A, B, use_prec) 666 im = mpf_add(C, D, use_prec) 667 if options.get('verbose'): 668 print("MUL: wanted", prec, "accurate bits, got", acc) 669 # multiply by I 670 if direction & 1: 671 re, im = mpf_neg(im), re 672 return re, im, acc, acc 673 674 675def evalf_pow(v, prec, options): 676 677 target_prec = prec 678 base, exp = v.args 679 680 # We handle x**n separately. This has two purposes: 1) it is much 681 # faster, because we avoid calling evalf on the exponent, and 2) it 682 # allows better handling of real/imaginary parts that are exactly zero 683 if exp.is_Integer: 684 p = exp.p 685 # Exact 686 if not p: 687 return fone, None, prec, None 688 # Exponentiation by p magnifies relative error by |p|, so the 689 # base must be evaluated with increased precision if p is large 690 prec += int(math.log(abs(p), 2)) 691 re, im, re_acc, im_acc = evalf(base, prec + 5, options) 692 # Real to integer power 693 if re and not im: 694 return mpf_pow_int(re, p, target_prec), None, target_prec, None 695 # (x*I)**n = I**n * x**n 696 if im and not re: 697 z = mpf_pow_int(im, p, target_prec) 698 case = p % 4 699 if case == 0: 700 return z, None, target_prec, None 701 if case == 1: 702 return None, z, None, target_prec 703 if case == 2: 704 return mpf_neg(z), None, target_prec, None 705 if case == 3: 706 return None, mpf_neg(z), None, target_prec 707 # Zero raised to an integer power 708 if not re: 709 return None, None, None, None 710 # General complex number to arbitrary integer power 711 re, im = libmp.mpc_pow_int((re, im), p, prec) 712 # Assumes full accuracy in input 713 return finalize_complex(re, im, target_prec) 714 715 # Pure square root 716 if exp is S.Half: 717 xre, xim, _, _ = evalf(base, prec + 5, options) 718 # General complex square root 719 if xim: 720 re, im = libmp.mpc_sqrt((xre or fzero, xim), prec) 721 return finalize_complex(re, im, prec) 722 if not xre: 723 return None, None, None, None 724 # Square root of a negative real number 725 if mpf_lt(xre, fzero): 726 return None, mpf_sqrt(mpf_neg(xre), prec), None, prec 727 # Positive square root 728 return mpf_sqrt(xre, prec), None, prec, None 729 730 # We first evaluate the exponent to find its magnitude 731 # This determines the working precision that must be used 732 prec += 10 733 yre, yim, _, _ = evalf(exp, prec, options) 734 # Special cases: x**0 735 if not (yre or yim): 736 return fone, None, prec, None 737 738 ysize = fastlog(yre) 739 # Restart if too big 740 # XXX: prec + ysize might exceed maxprec 741 if ysize > 5: 742 prec += ysize 743 yre, yim, _, _ = evalf(exp, prec, options) 744 745 # Pure exponential function; no need to evalf the base 746 if base is S.Exp1: 747 if yim: 748 re, im = libmp.mpc_exp((yre or fzero, yim), prec) 749 return finalize_complex(re, im, target_prec) 750 return mpf_exp(yre, target_prec), None, target_prec, None 751 752 xre, xim, _, _ = evalf(base, prec + 5, options) 753 # 0**y 754 if not (xre or xim): 755 return None, None, None, None 756 757 # (real ** complex) or (complex ** complex) 758 if yim: 759 re, im = libmp.mpc_pow( 760 (xre or fzero, xim or fzero), (yre or fzero, yim), 761 target_prec) 762 return finalize_complex(re, im, target_prec) 763 # complex ** real 764 if xim: 765 re, im = libmp.mpc_pow_mpf((xre or fzero, xim), yre, target_prec) 766 return finalize_complex(re, im, target_prec) 767 # negative ** real 768 elif mpf_lt(xre, fzero): 769 re, im = libmp.mpc_pow_mpf((xre, fzero), yre, target_prec) 770 return finalize_complex(re, im, target_prec) 771 # positive ** real 772 else: 773 return mpf_pow(xre, yre, target_prec), None, target_prec, None 774 775 776#----------------------------------------------------------------------------# 777# # 778# Special functions # 779# # 780#----------------------------------------------------------------------------# 781def evalf_trig(v, prec, options): 782 """ 783 This function handles sin and cos of complex arguments. 784 785 TODO: should also handle tan of complex arguments. 786 """ 787 from sympy import cos, sin 788 if isinstance(v, cos): 789 func = mpf_cos 790 elif isinstance(v, sin): 791 func = mpf_sin 792 else: 793 raise NotImplementedError 794 arg = v.args[0] 795 # 20 extra bits is possibly overkill. It does make the need 796 # to restart very unlikely 797 xprec = prec + 20 798 re, im, re_acc, im_acc = evalf(arg, xprec, options) 799 if im: 800 if 'subs' in options: 801 v = v.subs(options['subs']) 802 return evalf(v._eval_evalf(prec), prec, options) 803 if not re: 804 if isinstance(v, cos): 805 return fone, None, prec, None 806 elif isinstance(v, sin): 807 return None, None, None, None 808 else: 809 raise NotImplementedError 810 # For trigonometric functions, we are interested in the 811 # fixed-point (absolute) accuracy of the argument. 812 xsize = fastlog(re) 813 # Magnitude <= 1.0. OK to compute directly, because there is no 814 # danger of hitting the first root of cos (with sin, magnitude 815 # <= 2.0 would actually be ok) 816 if xsize < 1: 817 return func(re, prec, rnd), None, prec, None 818 # Very large 819 if xsize >= 10: 820 xprec = prec + xsize 821 re, im, re_acc, im_acc = evalf(arg, xprec, options) 822 # Need to repeat in case the argument is very close to a 823 # multiple of pi (or pi/2), hitting close to a root 824 while 1: 825 y = func(re, prec, rnd) 826 ysize = fastlog(y) 827 gap = -ysize 828 accuracy = (xprec - xsize) - gap 829 if accuracy < prec: 830 if options.get('verbose'): 831 print("SIN/COS", accuracy, "wanted", prec, "gap", gap) 832 print(to_str(y, 10)) 833 if xprec > options.get('maxprec', DEFAULT_MAXPREC): 834 return y, None, accuracy, None 835 xprec += gap 836 re, im, re_acc, im_acc = evalf(arg, xprec, options) 837 continue 838 else: 839 return y, None, prec, None 840 841 842def evalf_log(expr, prec, options): 843 from sympy import Abs, Add, log 844 if len(expr.args)>1: 845 expr = expr.doit() 846 return evalf(expr, prec, options) 847 arg = expr.args[0] 848 workprec = prec + 10 849 xre, xim, xacc, _ = evalf(arg, workprec, options) 850 851 # evalf can return NoneTypes if chop=True 852 # issue 18516, 19623 853 if xre is xim is None: 854 # Dear reviewer, I do not know what -inf is; 855 # it looks to be (1, 0, -789, -3) 856 # but I'm not sure in general, 857 # so we just let mpmath figure 858 # it out by taking log of 0 directly. 859 # It would be better to return -inf instead. 860 xre = fzero 861 862 if xim: 863 # XXX: use get_abs etc instead 864 re = evalf_log( 865 log(Abs(arg, evaluate=False), evaluate=False), prec, options) 866 im = mpf_atan2(xim, xre or fzero, prec) 867 return re[0], im, re[2], prec 868 869 imaginary_term = (mpf_cmp(xre, fzero) < 0) 870 871 re = mpf_log(mpf_abs(xre), prec, rnd) 872 size = fastlog(re) 873 if prec - size > workprec and re != fzero: 874 # We actually need to compute 1+x accurately, not x 875 arg = Add(S.NegativeOne, arg, evaluate=False) 876 xre, xim, _, _ = evalf_add(arg, prec, options) 877 prec2 = workprec - fastlog(xre) 878 # xre is now x - 1 so we add 1 back here to calculate x 879 re = mpf_log(mpf_abs(mpf_add(xre, fone, prec2)), prec, rnd) 880 881 re_acc = prec 882 883 if imaginary_term: 884 return re, mpf_pi(prec), re_acc, prec 885 else: 886 return re, None, re_acc, None 887 888 889def evalf_atan(v, prec, options): 890 arg = v.args[0] 891 xre, xim, reacc, imacc = evalf(arg, prec + 5, options) 892 if xre is xim is None: 893 return (None,)*4 894 if xim: 895 raise NotImplementedError 896 return mpf_atan(xre, prec, rnd), None, prec, None 897 898 899def evalf_subs(prec, subs): 900 """ Change all Float entries in `subs` to have precision prec. """ 901 newsubs = {} 902 for a, b in subs.items(): 903 b = S(b) 904 if b.is_Float: 905 b = b._eval_evalf(prec) 906 newsubs[a] = b 907 return newsubs 908 909 910def evalf_piecewise(expr, prec, options): 911 from sympy import Float, Integer 912 if 'subs' in options: 913 expr = expr.subs(evalf_subs(prec, options['subs'])) 914 newopts = options.copy() 915 del newopts['subs'] 916 if hasattr(expr, 'func'): 917 return evalf(expr, prec, newopts) 918 if type(expr) == float: 919 return evalf(Float(expr), prec, newopts) 920 if type(expr) == int: 921 return evalf(Integer(expr), prec, newopts) 922 923 # We still have undefined symbols 924 raise NotImplementedError 925 926 927def evalf_bernoulli(expr, prec, options): 928 arg = expr.args[0] 929 if not arg.is_Integer: 930 raise ValueError("Bernoulli number index must be an integer") 931 n = int(arg) 932 b = mpf_bernoulli(n, prec, rnd) 933 if b == fzero: 934 return None, None, None, None 935 return b, None, prec, None 936 937#----------------------------------------------------------------------------# 938# # 939# High-level operations # 940# # 941#----------------------------------------------------------------------------# 942 943 944def as_mpmath(x, prec, options): 945 from sympy.core.numbers import Infinity, NegativeInfinity, Zero 946 x = sympify(x) 947 if isinstance(x, Zero) or x == 0: 948 return mpf(0) 949 if isinstance(x, Infinity): 950 return mpf('inf') 951 if isinstance(x, NegativeInfinity): 952 return mpf('-inf') 953 # XXX 954 re, im, _, _ = evalf(x, prec, options) 955 if im: 956 return mpc(re or fzero, im) 957 return mpf(re) 958 959 960def do_integral(expr, prec, options): 961 func = expr.args[0] 962 x, xlow, xhigh = expr.args[1] 963 if xlow == xhigh: 964 xlow = xhigh = 0 965 elif x not in func.free_symbols: 966 # only the difference in limits matters in this case 967 # so if there is a symbol in common that will cancel 968 # out when taking the difference, then use that 969 # difference 970 if xhigh.free_symbols & xlow.free_symbols: 971 diff = xhigh - xlow 972 if diff.is_number: 973 xlow, xhigh = 0, diff 974 975 oldmaxprec = options.get('maxprec', DEFAULT_MAXPREC) 976 options['maxprec'] = min(oldmaxprec, 2*prec) 977 978 with workprec(prec + 5): 979 xlow = as_mpmath(xlow, prec + 15, options) 980 xhigh = as_mpmath(xhigh, prec + 15, options) 981 982 # Integration is like summation, and we can phone home from 983 # the integrand function to update accuracy summation style 984 # Note that this accuracy is inaccurate, since it fails 985 # to account for the variable quadrature weights, 986 # but it is better than nothing 987 988 from sympy import cos, sin, Wild 989 990 have_part = [False, False] 991 max_real_term = [MINUS_INF] 992 max_imag_term = [MINUS_INF] 993 994 def f(t): 995 re, im, re_acc, im_acc = evalf(func, mp.prec, {'subs': {x: t}}) 996 997 have_part[0] = re or have_part[0] 998 have_part[1] = im or have_part[1] 999 1000 max_real_term[0] = max(max_real_term[0], fastlog(re)) 1001 max_imag_term[0] = max(max_imag_term[0], fastlog(im)) 1002 1003 if im: 1004 return mpc(re or fzero, im) 1005 return mpf(re or fzero) 1006 1007 if options.get('quad') == 'osc': 1008 A = Wild('A', exclude=[x]) 1009 B = Wild('B', exclude=[x]) 1010 D = Wild('D') 1011 m = func.match(cos(A*x + B)*D) 1012 if not m: 1013 m = func.match(sin(A*x + B)*D) 1014 if not m: 1015 raise ValueError("An integrand of the form sin(A*x+B)*f(x) " 1016 "or cos(A*x+B)*f(x) is required for oscillatory quadrature") 1017 period = as_mpmath(2*S.Pi/m[A], prec + 15, options) 1018 result = quadosc(f, [xlow, xhigh], period=period) 1019 # XXX: quadosc does not do error detection yet 1020 quadrature_error = MINUS_INF 1021 else: 1022 result, quadrature_error = quadts(f, [xlow, xhigh], error=1) 1023 quadrature_error = fastlog(quadrature_error._mpf_) 1024 1025 options['maxprec'] = oldmaxprec 1026 1027 if have_part[0]: 1028 re = result.real._mpf_ 1029 if re == fzero: 1030 re, re_acc = scaled_zero( 1031 min(-prec, -max_real_term[0], -quadrature_error)) 1032 re = scaled_zero(re) # handled ok in evalf_integral 1033 else: 1034 re_acc = -max(max_real_term[0] - fastlog(re) - 1035 prec, quadrature_error) 1036 else: 1037 re, re_acc = None, None 1038 1039 if have_part[1]: 1040 im = result.imag._mpf_ 1041 if im == fzero: 1042 im, im_acc = scaled_zero( 1043 min(-prec, -max_imag_term[0], -quadrature_error)) 1044 im = scaled_zero(im) # handled ok in evalf_integral 1045 else: 1046 im_acc = -max(max_imag_term[0] - fastlog(im) - 1047 prec, quadrature_error) 1048 else: 1049 im, im_acc = None, None 1050 1051 result = re, im, re_acc, im_acc 1052 return result 1053 1054 1055def evalf_integral(expr, prec, options): 1056 limits = expr.limits 1057 if len(limits) != 1 or len(limits[0]) != 3: 1058 raise NotImplementedError 1059 workprec = prec 1060 i = 0 1061 maxprec = options.get('maxprec', INF) 1062 while 1: 1063 result = do_integral(expr, workprec, options) 1064 accuracy = complex_accuracy(result) 1065 if accuracy >= prec: # achieved desired precision 1066 break 1067 if workprec >= maxprec: # can't increase accuracy any more 1068 break 1069 if accuracy == -1: 1070 # maybe the answer really is zero and maybe we just haven't increased 1071 # the precision enough. So increase by doubling to not take too long 1072 # to get to maxprec. 1073 workprec *= 2 1074 else: 1075 workprec += max(prec, 2**i) 1076 workprec = min(workprec, maxprec) 1077 i += 1 1078 return result 1079 1080 1081def check_convergence(numer, denom, n): 1082 """ 1083 Returns 1084 ======= 1085 1086 (h, g, p) where 1087 -- h is: 1088 > 0 for convergence of rate 1/factorial(n)**h 1089 < 0 for divergence of rate factorial(n)**(-h) 1090 = 0 for geometric or polynomial convergence or divergence 1091 1092 -- abs(g) is: 1093 > 1 for geometric convergence of rate 1/h**n 1094 < 1 for geometric divergence of rate h**n 1095 = 1 for polynomial convergence or divergence 1096 1097 (g < 0 indicates an alternating series) 1098 1099 -- p is: 1100 > 1 for polynomial convergence of rate 1/n**h 1101 <= 1 for polynomial divergence of rate n**(-h) 1102 1103 """ 1104 from sympy import Poly 1105 npol = Poly(numer, n) 1106 dpol = Poly(denom, n) 1107 p = npol.degree() 1108 q = dpol.degree() 1109 rate = q - p 1110 if rate: 1111 return rate, None, None 1112 constant = dpol.LC() / npol.LC() 1113 if abs(constant) != 1: 1114 return rate, constant, None 1115 if npol.degree() == dpol.degree() == 0: 1116 return rate, constant, 0 1117 pc = npol.all_coeffs()[1] 1118 qc = dpol.all_coeffs()[1] 1119 return rate, constant, (qc - pc)/dpol.LC() 1120 1121 1122def hypsum(expr, n, start, prec): 1123 """ 1124 Sum a rapidly convergent infinite hypergeometric series with 1125 given general term, e.g. e = hypsum(1/factorial(n), n). The 1126 quotient between successive terms must be a quotient of integer 1127 polynomials. 1128 """ 1129 from sympy import Float, hypersimp, lambdify 1130 1131 if prec == float('inf'): 1132 raise NotImplementedError('does not support inf prec') 1133 1134 if start: 1135 expr = expr.subs(n, n + start) 1136 hs = hypersimp(expr, n) 1137 if hs is None: 1138 raise NotImplementedError("a hypergeometric series is required") 1139 num, den = hs.as_numer_denom() 1140 1141 func1 = lambdify(n, num) 1142 func2 = lambdify(n, den) 1143 1144 h, g, p = check_convergence(num, den, n) 1145 1146 if h < 0: 1147 raise ValueError("Sum diverges like (n!)^%i" % (-h)) 1148 1149 term = expr.subs(n, 0) 1150 if not term.is_Rational: 1151 raise NotImplementedError("Non rational term functionality is not implemented.") 1152 1153 # Direct summation if geometric or faster 1154 if h > 0 or (h == 0 and abs(g) > 1): 1155 term = (MPZ(term.p) << prec) // term.q 1156 s = term 1157 k = 1 1158 while abs(term) > 5: 1159 term *= MPZ(func1(k - 1)) 1160 term //= MPZ(func2(k - 1)) 1161 s += term 1162 k += 1 1163 return from_man_exp(s, -prec) 1164 else: 1165 alt = g < 0 1166 if abs(g) < 1: 1167 raise ValueError("Sum diverges like (%i)^n" % abs(1/g)) 1168 if p < 1 or (p == 1 and not alt): 1169 raise ValueError("Sum diverges like n^%i" % (-p)) 1170 # We have polynomial convergence: use Richardson extrapolation 1171 vold = None 1172 ndig = prec_to_dps(prec) 1173 while True: 1174 # Need to use at least quad precision because a lot of cancellation 1175 # might occur in the extrapolation process; we check the answer to 1176 # make sure that the desired precision has been reached, too. 1177 prec2 = 4*prec 1178 term0 = (MPZ(term.p) << prec2) // term.q 1179 1180 def summand(k, _term=[term0]): 1181 if k: 1182 k = int(k) 1183 _term[0] *= MPZ(func1(k - 1)) 1184 _term[0] //= MPZ(func2(k - 1)) 1185 return make_mpf(from_man_exp(_term[0], -prec2)) 1186 1187 with workprec(prec): 1188 v = nsum(summand, [0, mpmath_inf], method='richardson') 1189 vf = Float(v, ndig) 1190 if vold is not None and vold == vf: 1191 break 1192 prec += prec # double precision each time 1193 vold = vf 1194 1195 return v._mpf_ 1196 1197 1198def evalf_prod(expr, prec, options): 1199 from sympy import Sum 1200 if all((l[1] - l[2]).is_Integer for l in expr.limits): 1201 re, im, re_acc, im_acc = evalf(expr.doit(), prec=prec, options=options) 1202 else: 1203 re, im, re_acc, im_acc = evalf(expr.rewrite(Sum), prec=prec, options=options) 1204 return re, im, re_acc, im_acc 1205 1206 1207def evalf_sum(expr, prec, options): 1208 from sympy import Float 1209 if 'subs' in options: 1210 expr = expr.subs(options['subs']) 1211 func = expr.function 1212 limits = expr.limits 1213 if len(limits) != 1 or len(limits[0]) != 3: 1214 raise NotImplementedError 1215 if func.is_zero: 1216 return None, None, prec, None 1217 prec2 = prec + 10 1218 try: 1219 n, a, b = limits[0] 1220 if b != S.Infinity or a != int(a): 1221 raise NotImplementedError 1222 # Use fast hypergeometric summation if possible 1223 v = hypsum(func, n, int(a), prec2) 1224 delta = prec - fastlog(v) 1225 if fastlog(v) < -10: 1226 v = hypsum(func, n, int(a), delta) 1227 return v, None, min(prec, delta), None 1228 except NotImplementedError: 1229 # Euler-Maclaurin summation for general series 1230 eps = Float(2.0)**(-prec) 1231 for i in range(1, 5): 1232 m = n = 2**i * prec 1233 s, err = expr.euler_maclaurin(m=m, n=n, eps=eps, 1234 eval_integral=False) 1235 err = err.evalf() 1236 if err <= eps: 1237 break 1238 err = fastlog(evalf(abs(err), 20, options)[0]) 1239 re, im, re_acc, im_acc = evalf(s, prec2, options) 1240 if re_acc is None: 1241 re_acc = -err 1242 if im_acc is None: 1243 im_acc = -err 1244 return re, im, re_acc, im_acc 1245 1246 1247#----------------------------------------------------------------------------# 1248# # 1249# Symbolic interface # 1250# # 1251#----------------------------------------------------------------------------# 1252 1253def evalf_symbol(x, prec, options): 1254 val = options['subs'][x] 1255 if isinstance(val, mpf): 1256 if not val: 1257 return None, None, None, None 1258 return val._mpf_, None, prec, None 1259 else: 1260 if not '_cache' in options: 1261 options['_cache'] = {} 1262 cache = options['_cache'] 1263 cached, cached_prec = cache.get(x, (None, MINUS_INF)) 1264 if cached_prec >= prec: 1265 return cached 1266 v = evalf(sympify(val), prec, options) 1267 cache[x] = (v, prec) 1268 return v 1269 1270evalf_table = None 1271 1272 1273def _create_evalf_table(): 1274 global evalf_table 1275 from sympy.functions.combinatorial.numbers import bernoulli 1276 from sympy.concrete.products import Product 1277 from sympy.concrete.summations import Sum 1278 from sympy.core.add import Add 1279 from sympy.core.mul import Mul 1280 from sympy.core.numbers import Exp1, Float, Half, ImaginaryUnit, Integer, NaN, NegativeOne, One, Pi, Rational, Zero 1281 from sympy.core.power import Pow 1282 from sympy.core.symbol import Dummy, Symbol 1283 from sympy.functions.elementary.complexes import Abs, im, re 1284 from sympy.functions.elementary.exponential import exp, log 1285 from sympy.functions.elementary.integers import ceiling, floor 1286 from sympy.functions.elementary.piecewise import Piecewise 1287 from sympy.functions.elementary.trigonometric import atan, cos, sin 1288 from sympy.integrals.integrals import Integral 1289 evalf_table = { 1290 Symbol: evalf_symbol, 1291 Dummy: evalf_symbol, 1292 Float: lambda x, prec, options: (x._mpf_, None, prec, None), 1293 Rational: lambda x, prec, options: (from_rational(x.p, x.q, prec), None, prec, None), 1294 Integer: lambda x, prec, options: (from_int(x.p, prec), None, prec, None), 1295 Zero: lambda x, prec, options: (None, None, prec, None), 1296 One: lambda x, prec, options: (fone, None, prec, None), 1297 Half: lambda x, prec, options: (fhalf, None, prec, None), 1298 Pi: lambda x, prec, options: (mpf_pi(prec), None, prec, None), 1299 Exp1: lambda x, prec, options: (mpf_e(prec), None, prec, None), 1300 ImaginaryUnit: lambda x, prec, options: (None, fone, None, prec), 1301 NegativeOne: lambda x, prec, options: (fnone, None, prec, None), 1302 NaN: lambda x, prec, options: (fnan, None, prec, None), 1303 1304 exp: lambda x, prec, options: evalf_pow( 1305 Pow(S.Exp1, x.exp, evaluate=False), prec, options), 1306 1307 cos: evalf_trig, 1308 sin: evalf_trig, 1309 1310 Add: evalf_add, 1311 Mul: evalf_mul, 1312 Pow: evalf_pow, 1313 1314 log: evalf_log, 1315 atan: evalf_atan, 1316 Abs: evalf_abs, 1317 1318 re: evalf_re, 1319 im: evalf_im, 1320 floor: evalf_floor, 1321 ceiling: evalf_ceiling, 1322 1323 Integral: evalf_integral, 1324 Sum: evalf_sum, 1325 Product: evalf_prod, 1326 Piecewise: evalf_piecewise, 1327 1328 bernoulli: evalf_bernoulli, 1329 } 1330 1331 1332def evalf(x, prec, options): 1333 """ 1334 Evaluate the ``Basic`` instance, ``x`` 1335 to a binary precision of ``prec``. This 1336 function is supposed to be used internally. 1337 1338 Parameters 1339 ========== 1340 1341 x : Basic 1342 The formula to evaluate to a float. 1343 prec : int 1344 The binary precision that the output should have. 1345 options : dict 1346 A dictionary with the same entries as 1347 ``EvalfMixin.evalf`` and in addition, 1348 ``maxprec`` which is the maximum working precision. 1349 1350 Returns 1351 ======= 1352 1353 An optional tuple, ``(re, im, re_acc, im_acc)`` 1354 which are the real, imaginary, real accuracy 1355 and imaginary accuracy respectively. ``re`` is 1356 an mpf value tuple and so is ``im``. ``re_acc`` 1357 and ``im_acc`` are ints. 1358 1359 NB: all these return values can be ``None``. 1360 If all values are ``None``, then that represents 0. 1361 Note that 0 is also represented as ``fzero = (0, 0, 0, 0)``. 1362 """ 1363 from sympy import re as re_, im as im_ 1364 try: 1365 rf = evalf_table[x.func] 1366 r = rf(x, prec, options) 1367 except KeyError: 1368 # Fall back to ordinary evalf if possible 1369 if 'subs' in options: 1370 x = x.subs(evalf_subs(prec, options['subs'])) 1371 xe = x._eval_evalf(prec) 1372 if xe is None: 1373 raise NotImplementedError 1374 as_real_imag = getattr(xe, "as_real_imag", None) 1375 if as_real_imag is None: 1376 raise NotImplementedError # e.g. FiniteSet(-1.0, 1.0).evalf() 1377 re, im = as_real_imag() 1378 if re.has(re_) or im.has(im_): 1379 raise NotImplementedError 1380 if re == 0: 1381 re = None 1382 reprec = None 1383 elif re.is_number: 1384 re = re._to_mpmath(prec, allow_ints=False)._mpf_ 1385 reprec = prec 1386 else: 1387 raise NotImplementedError 1388 if im == 0: 1389 im = None 1390 imprec = None 1391 elif im.is_number: 1392 im = im._to_mpmath(prec, allow_ints=False)._mpf_ 1393 imprec = prec 1394 else: 1395 raise NotImplementedError 1396 r = re, im, reprec, imprec 1397 1398 if options.get("verbose"): 1399 print("### input", x) 1400 print("### output", to_str(r[0] or fzero, 50)) 1401 print("### raw", r) # r[0], r[2] 1402 print() 1403 chop = options.get('chop', False) 1404 if chop: 1405 if chop is True: 1406 chop_prec = prec 1407 else: 1408 # convert (approximately) from given tolerance; 1409 # the formula here will will make 1e-i rounds to 0 for 1410 # i in the range +/-27 while 2e-i will not be chopped 1411 chop_prec = int(round(-3.321*math.log10(chop) + 2.5)) 1412 if chop_prec == 3: 1413 chop_prec -= 1 1414 r = chop_parts(r, chop_prec) 1415 if options.get("strict"): 1416 check_target(x, r, prec) 1417 return r 1418 1419 1420class EvalfMixin: 1421 """Mixin class adding evalf capabililty.""" 1422 1423 __slots__ = () # type: Tuple[str, ...] 1424 1425 def evalf(self, n=15, subs=None, maxn=100, chop=False, strict=False, quad=None, verbose=False): 1426 """ 1427 Evaluate the given formula to an accuracy of *n* digits. 1428 1429 Parameters 1430 ========== 1431 1432 subs : dict, optional 1433 Substitute numerical values for symbols, e.g. 1434 ``subs={x:3, y:1+pi}``. The substitutions must be given as a 1435 dictionary. 1436 1437 maxn : int, optional 1438 Allow a maximum temporary working precision of maxn digits. 1439 1440 chop : bool or number, optional 1441 Specifies how to replace tiny real or imaginary parts in 1442 subresults by exact zeros. 1443 1444 When ``True`` the chop value defaults to standard precision. 1445 1446 Otherwise the chop value is used to determine the 1447 magnitude of "small" for purposes of chopping. 1448 1449 >>> from sympy import N 1450 >>> x = 1e-4 1451 >>> N(x, chop=True) 1452 0.000100000000000000 1453 >>> N(x, chop=1e-5) 1454 0.000100000000000000 1455 >>> N(x, chop=1e-4) 1456 0 1457 1458 strict : bool, optional 1459 Raise ``PrecisionExhausted`` if any subresult fails to 1460 evaluate to full accuracy, given the available maxprec. 1461 1462 quad : str, optional 1463 Choose algorithm for numerical quadrature. By default, 1464 tanh-sinh quadrature is used. For oscillatory 1465 integrals on an infinite interval, try ``quad='osc'``. 1466 1467 verbose : bool, optional 1468 Print debug information. 1469 1470 Notes 1471 ===== 1472 1473 When Floats are naively substituted into an expression, 1474 precision errors may adversely affect the result. For example, 1475 adding 1e16 (a Float) to 1 will truncate to 1e16; if 1e16 is 1476 then subtracted, the result will be 0. 1477 That is exactly what happens in the following: 1478 1479 >>> from sympy.abc import x, y, z 1480 >>> values = {x: 1e16, y: 1, z: 1e16} 1481 >>> (x + y - z).subs(values) 1482 0 1483 1484 Using the subs argument for evalf is the accurate way to 1485 evaluate such an expression: 1486 1487 >>> (x + y - z).evalf(subs=values) 1488 1.00000000000000 1489 """ 1490 from sympy import Float, Number 1491 n = n if n is not None else 15 1492 1493 if subs and is_sequence(subs): 1494 raise TypeError('subs must be given as a dictionary') 1495 1496 # for sake of sage that doesn't like evalf(1) 1497 if n == 1 and isinstance(self, Number): 1498 from sympy.core.expr import _mag 1499 rv = self.evalf(2, subs, maxn, chop, strict, quad, verbose) 1500 m = _mag(rv) 1501 rv = rv.round(1 - m) 1502 return rv 1503 1504 if not evalf_table: 1505 _create_evalf_table() 1506 prec = dps_to_prec(n) 1507 options = {'maxprec': max(prec, int(maxn*LG10)), 'chop': chop, 1508 'strict': strict, 'verbose': verbose} 1509 if subs is not None: 1510 options['subs'] = subs 1511 if quad is not None: 1512 options['quad'] = quad 1513 try: 1514 result = evalf(self, prec + 4, options) 1515 except NotImplementedError: 1516 # Fall back to the ordinary evalf 1517 if hasattr(self, 'subs') and subs is not None: # issue 20291 1518 v = self.subs(subs)._eval_evalf(prec) 1519 else: 1520 v = self._eval_evalf(prec) 1521 if v is None: 1522 return self 1523 elif not v.is_number: 1524 return v 1525 try: 1526 # If the result is numerical, normalize it 1527 result = evalf(v, prec, options) 1528 except NotImplementedError: 1529 # Probably contains symbols or unknown functions 1530 return v 1531 re, im, re_acc, im_acc = result 1532 if re: 1533 p = max(min(prec, re_acc), 1) 1534 re = Float._new(re, p) 1535 else: 1536 re = S.Zero 1537 if im: 1538 p = max(min(prec, im_acc), 1) 1539 im = Float._new(im, p) 1540 return re + im*S.ImaginaryUnit 1541 else: 1542 return re 1543 1544 n = evalf 1545 1546 def _evalf(self, prec): 1547 """Helper for evalf. Does the same thing but takes binary precision""" 1548 r = self._eval_evalf(prec) 1549 if r is None: 1550 r = self 1551 return r 1552 1553 def _eval_evalf(self, prec): 1554 return 1555 1556 def _to_mpmath(self, prec, allow_ints=True): 1557 # mpmath functions accept ints as input 1558 errmsg = "cannot convert to mpmath number" 1559 if allow_ints and self.is_Integer: 1560 return self.p 1561 if hasattr(self, '_as_mpf_val'): 1562 return make_mpf(self._as_mpf_val(prec)) 1563 try: 1564 re, im, _, _ = evalf(self, prec, {}) 1565 if im: 1566 if not re: 1567 re = fzero 1568 return make_mpc((re, im)) 1569 elif re: 1570 return make_mpf(re) 1571 else: 1572 return make_mpf(fzero) 1573 except NotImplementedError: 1574 v = self._eval_evalf(prec) 1575 if v is None: 1576 raise ValueError(errmsg) 1577 if v.is_Float: 1578 return make_mpf(v._mpf_) 1579 # Number + Number*I is also fine 1580 re, im = v.as_real_imag() 1581 if allow_ints and re.is_Integer: 1582 re = from_int(re.p) 1583 elif re.is_Float: 1584 re = re._mpf_ 1585 else: 1586 raise ValueError(errmsg) 1587 if allow_ints and im.is_Integer: 1588 im = from_int(im.p) 1589 elif im.is_Float: 1590 im = im._mpf_ 1591 else: 1592 raise ValueError(errmsg) 1593 return make_mpc((re, im)) 1594 1595 1596def N(x, n=15, **options): 1597 r""" 1598 Calls x.evalf(n, \*\*options). 1599 1600 Explanations 1601 ============ 1602 1603 Both .n() and N() are equivalent to .evalf(); use the one that you like better. 1604 See also the docstring of .evalf() for information on the options. 1605 1606 Examples 1607 ======== 1608 1609 >>> from sympy import Sum, oo, N 1610 >>> from sympy.abc import k 1611 >>> Sum(1/k**k, (k, 1, oo)) 1612 Sum(k**(-k), (k, 1, oo)) 1613 >>> N(_, 4) 1614 1.291 1615 1616 """ 1617 # by using rational=True, any evaluation of a string 1618 # will be done using exact values for the Floats 1619 return sympify(x, rational=True).evalf(n, **options) 1620