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