2 Adaptive numerical evaluation of SymPy expressions, using mpmath
3 for mathematical functions.
6 from sympy
.mpmath
.lib
import (from_int
, from_rational
, fpi
, fzero
, fcmp
,
7 normalize
, bitcount
, round_nearest
, to_str
, fpow
, fone
, fpowi
, fe
,
8 fnone
, fhalf
, fcos
, fsin
, flog
, fatan
, fmul
, fneg
, to_float
, fshift
,
9 fabs
, fatan2
, fadd
, fdiv
, flt
, dps_to_prec
, prec_to_dps
, fpos
, from_float
,
12 import sympy
.mpmath
.libmpc
as libmpc
13 from sympy
.mpmath
import mpf
, mpc
, quadts
, quadosc
, mp
, make_mpf
14 from sympy
.mpmath
.specfun
import mpf_gamma
15 from sympy
.mpmath
.lib
import MP_BASE
, from_man_exp
16 from sympy
.mpmath
.calculus
import shanks_extrapolation
, richardson_extrapolation
21 from basic
import Basic
, C
, S
22 from function
import Function
23 from sympify
import sympify
27 # Used in a few places as placeholder values to denote exponents and
28 # precision levels, e.g. of exact numbers. Must be careful to avoid
29 # passing these to mpmath functions or returning them in final results.
33 # ~= 100 digits. Real men set this to INF.
36 class PrecisionExhausted(ArithmeticError):
39 #----------------------------------------------------------------------------#
41 # Helper functions for arithmetic and complex parts #
43 #----------------------------------------------------------------------------#
46 An mpf value tuple is a tuple of integers (sign, man, exp, bc)
47 representing a floating-point numbers.
49 A temporary result is a tuple (re, im, re_acc, im_acc) where
50 re and im are nonzero mpf value tuples representing approximate
51 numbers, or None to denote exact zeros.
53 re_acc, im_acc are integers denoting log2(e) where is the estimated
54 relative accuracy of the respective complex part, but many be anything
55 if the corresponding complex part is None.
60 """Fast approximation of log2(x) for an mpf value tuple x."""
61 if not x
or x
== fzero
:
63 # log2(x) ~= exponent + width of mantissa
64 # Note: this actually gives ceil(log2(x)), which is a useful
65 # feature for interval arithmetic.
68 def complex_accuracy(result
):
70 Returns relative accuracy of a complex number with given accuracies
71 for the real and imaginary parts. The relative accuracy is defined
72 in the complex norm sense as ||z|+|error|| / |z| where error
73 is equal to (real absolute error) + (imag absolute error)*i.
75 The full expression for the (logarithmic) error can be approximated
76 easily by using the max norm to approximate the complex norm.
78 In the worst case (re and im equal), this is wrong by a factor
79 sqrt(2), or by log2(sqrt(2)) = 0.5 bit.
81 re
, im
, re_acc
, im_acc
= result
90 absolute_error
= max(re_size
-re_acc
, im_size
-im_acc
)
91 relative_error
= absolute_error
- max(re_size
, im_size
)
92 return -relative_error
94 def get_abs(expr
, prec
, options
):
95 re
, im
, re_acc
, im_acc
= evalf(expr
, prec
+2, options
)
97 re
, re_acc
= im
, im_acc
99 return libmpc
.mpc_abs((re
, im
), prec
), None, re_acc
, None
101 return fabs(re
), None, re_acc
, None
103 def get_complex_part(expr
, no
, prec
, options
):
104 """no = 0 for real part, no = 1 for imaginary part"""
108 res
= evalf(expr
, workprec
, options
)
109 value
, accuracy
= res
[no
::2]
110 if (not value
) or accuracy
>= prec
:
111 return value
, None, accuracy
, None
112 workprec
+= max(30, 2**i
)
115 def evalf_abs(expr
, prec
, options
):
116 return get_abs(expr
.args
[0], prec
, options
)
118 def evalf_re(expr
, prec
, options
):
119 return get_complex_part(expr
.args
[0], 0, prec
, options
)
121 def evalf_im(expr
, prec
, options
):
122 return get_complex_part(expr
.args
[0], 1, prec
, options
)
124 def finalize_complex(re
, im
, prec
):
126 if re
== fzero
and im
== fzero
:
127 raise ValueError("got complex zero with unknown accuracy")
128 size_re
= fastlog(re
)
129 size_im
= fastlog(im
)
130 # Convert fzeros to scaled zeros
132 re
= fshift(fone
, size_im
-prec
)
133 size_re
= fastlog(re
)
135 im
= fshift(fone
, size_re
-prec
)
136 size_im
= fastlog(im
)
137 if size_re
> size_im
:
139 im_acc
= prec
+ min(-(size_re
- size_im
), 0)
142 re_acc
= prec
+ min(-(size_im
- size_re
), 0)
143 return re
, im
, re_acc
, im_acc
145 def chop_parts(value
, prec
):
147 Chop off tiny real or complex parts.
149 re
, im
, re_acc
, im_acc
= value
150 # Method 1: chop based on absolute value
151 if re
and (fastlog(re
) < -prec
+4):
152 re
, re_acc
= None, None
153 if im
and (fastlog(im
) < -prec
+4):
154 im
, im_acc
= None, None
155 # Method 2: chop if inaccurate and relatively small
157 delta
= fastlog(re
) - fastlog(im
)
158 if re_acc
< 2 and (delta
- re_acc
<= -prec
+4):
159 re
, re_acc
= None, None
160 if im_acc
< 2 and (delta
- im_acc
>= prec
-4):
161 im
, im_acc
= None, None
162 return re
, im
, re_acc
, im_acc
164 def check_target(expr
, result
, prec
):
165 a
= complex_accuracy(result
)
167 raise PrecisionExhausted("Failed to distinguish the expression: \n\n%s\n\n"
168 "from zero. Try simplifying the input, using chop=True, or providing "
169 "a higher maxprec for evalf" % (expr
))
171 def get_integer_part(expr
, no
, options
, return_ints
=False):
173 With no = 1, computes ceiling(expr)
174 With no = -1, computes floor(expr)
176 Note: this function either gives the exact result or signals failure.
178 ire
, iim
, ire_acc
, iim_acc
= evalf(expr
, 30, options
)
180 gap
= max(fastlog(ire
)-ire_acc
, fastlog(iim
)-iim_acc
)
182 gap
= fastlog(ire
)-ire_acc
184 gap
= fastlog(iim
)-iim_acc
186 return None, None, None, None
188 ire
, iim
, ire_acc
, iim_acc
= evalf(expr
, 30+gap
, options
)
189 if ire
: nint_re
= int(to_int(ire
, round_nearest
))
190 if iim
: nint_im
= int(to_int(iim
, round_nearest
))
191 re
, im
, re_acc
, im_acc
= None, None, None, None
193 e
= C
.Add(C
.re(expr
, evaluate
=False), -nint_re
, evaluate
=False)
194 re
, _
, re_acc
, _
= evalf(e
, 10, options
)
195 check_target(e
, (re
, None, re_acc
, None), 3)
196 #assert (re_acc - fastlog(re)) > 3
197 nint_re
+= int(no
*(fcmp(re
or fzero
, fzero
) == no
))
198 re
= from_int(nint_re
)
200 e
= C
.Add(C
.im(expr
, evaluate
=False), -nint_im
*S
.ImaginaryUnit
, evaluate
=False)
201 _
, im
, _
, im_acc
= evalf(e
, 10, options
)
202 check_target(e
, (im
, None, im_acc
, None), 3)
203 #assert (im_acc - fastlog(im)) > 3
204 nint_im
+= int(no
*(fcmp(im
or fzero
, fzero
) == no
))
205 im
= from_int(nint_im
)
207 return int(to_int(re
or fzero
)), int(to_int(im
or fzero
))
208 return re
, im
, re_acc
, im_acc
210 def evalf_ceiling(expr
, prec
, options
):
211 return get_integer_part(expr
.args
[0], 1, options
)
213 def evalf_floor(expr
, prec
, options
):
214 return get_integer_part(expr
.args
[0], -1, options
)
216 #----------------------------------------------------------------------------#
218 # Arithmetic operations #
220 #----------------------------------------------------------------------------#
222 def add_terms(terms
, prec
, target_prec
):
224 Helper for evalf_add. Adds a list of (mpfval, accuracy) terms.
228 # XXX: this is supposed to represent a scaled zero
229 return fshift(fone
, target_prec
), -1
231 sum_man
, sum_exp
, absolute_error
= 0, 0, MINUS_INF
232 for x
, accuracy
in terms
:
235 sign
, man
, exp
, bc
= x
238 absolute_error
= max(absolute_error
, bc
+exp
-accuracy
)
239 delta
= exp
- sum_exp
245 sum_man
+= man
<< delta
247 if (-delta
) > 4*prec
:
250 sum_man
= (sum_man
<< (-delta
)) + man
252 if absolute_error
== MINUS_INF
:
255 # XXX: this is supposed to represent a scaled zero
256 return fshift(fone
, absolute_error
), -1
262 sum_bc
= bitcount(sum_man
)
263 sum_accuracy
= sum_exp
+ sum_bc
- absolute_error
264 r
= normalize(sum_sign
, sum_man
, sum_exp
, sum_bc
, target_prec
,
265 round_nearest
), sum_accuracy
266 #print "returning", to_str(r[0],50), r[1]
269 def evalf_add(v
, prec
, options
):
274 oldmaxprec
= options
.get('maxprec', DEFAULT_MAXPREC
)
275 options
['maxprec'] = min(oldmaxprec
, 2*prec
)
279 terms
= [evalf(arg
, prec
+10, options
) for arg
in args
]
280 re
, re_accuracy
= add_terms([(a
[0],a
[2]) for a
in terms
if a
[0]], prec
, target_prec
)
281 im
, im_accuracy
= add_terms([(a
[1],a
[3]) for a
in terms
if a
[1]], prec
, target_prec
)
282 accuracy
= complex_accuracy((re
, im
, re_accuracy
, im_accuracy
))
283 if accuracy
>= target_prec
:
284 if options
.get('verbose'):
285 print "ADD: wanted", target_prec
, "accurate bits, got", re_accuracy
, im_accuracy
286 return re
, im
, re_accuracy
, im_accuracy
288 diff
= target_prec
- accuracy
289 if (prec
-target_prec
) > options
.get('maxprec', DEFAULT_MAXPREC
):
290 return re
, im
, re_accuracy
, im_accuracy
292 prec
= prec
+ max(10+2**i
, diff
)
293 options
['maxprec'] = min(oldmaxprec
, 2*prec
)
294 if options
.get('verbose'):
295 print "ADD: restarting with prec", prec
298 options
['maxprec'] = oldmaxprec
300 # Helper for complex multiplication
301 # XXX: should be able to multiply directly, and use complex_accuracy
302 # to obtain the final accuracy
303 def cmul((a
, aacc
), (b
, bacc
), (c
, cacc
), (d
, dacc
), prec
, target_prec
):
304 A
, Aacc
= fmul(a
,c
,prec
), min(aacc
, cacc
)
305 B
, Bacc
= fmul(fneg(b
),d
,prec
), min(bacc
, dacc
)
306 C
, Cacc
= fmul(a
,d
,prec
), min(aacc
, dacc
)
307 D
, Dacc
= fmul(b
,c
,prec
), min(bacc
, cacc
)
308 re
, re_accuracy
= add_terms([(A
, Aacc
), (B
, Bacc
)], prec
, target_prec
)
309 im
, im_accuracy
= add_terms([(C
, Cacc
), (D
, Cacc
)], prec
, target_prec
)
310 return re
, im
, re_accuracy
, im_accuracy
312 def evalf_mul(v
, prec
, options
):
314 # With guard digits, multiplication in the real case does not destroy
315 # accuracy. This is also true in the complex case when considering the
316 # total accuracy; however accuracy for the real or imaginary parts
317 # separately may be lower.
320 # XXX: big overestimate
321 prec
= prec
+ len(args
) + 5
324 man
, exp
, bc
= 1, 0, 1
327 # First, we multiply all pure real or pure imaginary numbers.
328 # direction tells us that the result should be multiplied by
331 re
, im
, a
, aim
= evalf(arg
, prec
, options
)
333 complex_factors
.append((re
, im
, a
, aim
))
342 return None, None, None, None
351 sign
= (direction
& 2) >> 1
352 v
= normalize(sign
, man
, exp
, bitcount(man
), prec
, round_nearest
)
354 # Multiply first complex number by the existing real scalar
355 re
, im
, re_acc
, im_acc
= complex_factors
[0]
356 re
= fmul(re
, v
, prec
)
357 im
= fmul(im
, v
, prec
)
358 re_acc
= min(re_acc
, acc
)
359 im_acc
= min(im_acc
, acc
)
360 # Multiply consecutive complex factors
361 complex_factors
= complex_factors
[1:]
362 for wre
, wim
, wre_acc
, wim_acc
in complex_factors
:
363 re
, im
, re_acc
, im_acc
= cmul((re
, re_acc
), (im
,im_acc
),
364 (wre
,wre_acc
), (wim
,wim_acc
), prec
, target_prec
)
365 if options
.get('verbose'):
366 print "MUL: obtained accuracy", re_acc
, im_acc
, "expected", target_prec
369 return fneg(im
), re
, re_acc
, im_acc
371 return re
, im
, re_acc
, im_acc
375 return None, v
, None, acc
377 return v
, None, acc
, None
379 def evalf_pow(v
, prec
, options
):
382 # We handle x**n separately. This has two purposes: 1) it is much
383 # faster, because we avoid calling evalf on the exponent, and 2) it
384 # allows better handling of real/imaginary parts that are exactly zero
389 return fone
, None, prec
, None
390 # Exponentiation by p magnifies relative error by |p|, so the
391 # base must be evaluated with increased precision if p is large
392 prec
+= int(math
.log(abs(p
),2))
393 re
, im
, re_acc
, im_acc
= evalf(base
, prec
+5, options
)
394 # Real to integer power
396 return fpowi(re
, p
, target_prec
), None, target_prec
, None
397 # (x*I)**n = I**n * x**n
399 z
= fpowi(im
, p
, target_prec
)
401 if case
== 0: return z
, None, target_prec
, None
402 if case
== 1: return None, z
, None, target_prec
403 if case
== 2: return fneg(z
), None, target_prec
, None
404 if case
== 3: return None, fneg(z
), None, target_prec
405 # General complex number to arbitrary integer power
406 re
, im
= libmpc
.mpc_pow_int((re
, im
), p
, prec
)
407 # Assumes full accuracy in input
408 return finalize_complex(re
, im
, target_prec
)
410 # TODO: optimize these cases
411 #pure_exp = (base is S.Exp1)
412 #pure_sqrt = (base is S.Half)
414 # Complex or real to a real power
415 # We first evaluate the exponent to find its magnitude
417 yre
, yim
, yre_acc
, yim_acc
= evalf(exp
, prec
, options
)
419 # Need to restart if too big
422 yre
, yim
, yre_acc
, yim_acc
= evalf(exp
, prec
, options
)
425 re
, im
= libmpc
.mpc_exp((yre
or fzero
, yim
or fzero
), prec
)
426 return finalize_complex(re
, im
, target_prec
)
427 raise NotImplementedError
428 xre
, xim
, xre_acc
, yim_acc
= evalf(base
, prec
, options
)
431 re
, im
= libmpc
.mpc_pow_mpf((xre
or fzero
, xim
), yre
, target_prec
)
432 return finalize_complex(re
, im
, target_prec
)
433 # Fractional power of negative real
434 elif flt(xre
, fzero
):
435 return None, fpow(fneg(xre
), yre
, target_prec
), None, target_prec
437 return fpow(xre
, yre
, target_prec
), None, target_prec
, None
442 #----------------------------------------------------------------------------#
444 # Special functions #
446 #----------------------------------------------------------------------------#
448 def evalf_trig(v
, prec
, options
):
450 This function handles sin and cos of real arguments.
452 TODO: should also handle tan and complex arguments.
456 elif v
.func
is C
.sin
:
459 raise NotImplementedError
461 # 20 extra bits is possibly overkill. It does make the need
462 # to restart very unlikely
464 re
, im
, re_accuracy
, im_accuracy
= evalf(arg
, xprec
, options
)
466 raise NotImplementedError
469 return fone
, None, prec
, None
470 elif v
.func
is C
.sin
:
471 return None, None, None, None
473 raise NotImplementedError
474 # For trigonometric functions, we are interested in the
475 # fixed-point (absolute) accuracy of the argument.
477 # Magnitude <= 1.0. OK to compute directly, because there is no
478 # danger of hitting the first root of cos (with sin, magnitude
479 # <= 2.0 would actually be ok)
481 return func(re
, prec
, round_nearest
), None, prec
, None
485 re
, im
, re_accuracy
, im_accuracy
= evalf(arg
, xprec
, options
)
486 # Need to repeat in case the argument is very close to a
487 # multiple of pi (or pi/2), hitting close to a root
489 y
= func(re
, prec
, round_nearest
)
492 accuracy
= (xprec
- xsize
) - gap
494 if options
.get('verbose'):
495 print "SIN/COS", accuracy
, "wanted", prec
, "gap", gap
497 if xprec
> options
.get('maxprec', DEFAULT_MAXPREC
):
498 return y
, None, accuracy
, None
500 re
, im
, re_accuracy
, im_accuracy
= evalf(arg
, xprec
, options
)
503 return y
, None, prec
, None
505 def evalf_log(expr
, prec
, options
):
508 xre
, xim
, xacc
, _
= evalf(arg
, workprec
, options
)
511 # XXX: use get_abs etc instead
512 re
= evalf_log(C
.log(C
.abs(arg
, evaluate
=False), evaluate
=False), prec
, options
)
513 im
= fatan2(xim
, xre
, prec
)
514 return re
[0], im
, re
[2], prec
516 imaginary_term
= (fcmp(xre
, fzero
) < 0)
518 re
= flog(fabs(xre
), prec
, round_nearest
)
520 if prec
- size
> workprec
:
521 # We actually need to compute 1+x accurately, not x
522 arg
= C
.Add(S
.NegativeOne
,arg
,evaluate
=False)
523 xre
, xim
, xre_acc
, xim_acc
= evalf_add(arg
, prec
, options
)
524 prec2
= workprec
- fastlog(xre
)
525 re
= flog(fadd(xre
, fone
, prec2
), prec
, round_nearest
)
530 return re
, fpi(prec
), re_acc
, prec
532 return re
, None, re_acc
, None
534 def evalf_atan(v
, prec
, options
):
536 xre
, xim
, reacc
, imacc
= evalf(arg
, prec
+5, options
)
538 raise NotImplementedError
539 return fatan(xre
, prec
, round_nearest
), None, prec
, None
542 #----------------------------------------------------------------------------#
544 # High-level operations #
546 #----------------------------------------------------------------------------#
548 def as_mpmath(x
, prec
, options
):
550 if isinstance(x
, C
.Zero
):
552 if isinstance(x
, C
.Infinity
):
554 if isinstance(x
, C
.NegativeInfinity
):
557 re
, im
, _
, _
= evalf(x
, prec
, options
)
559 return mpc(re
or fzero
, im
)
562 def do_integral(expr
, prec
, options
):
564 x
, (xlow
, xhigh
) = expr
.args
[1][0]
567 oldmaxprec
= options
.get('maxprec', DEFAULT_MAXPREC
)
568 options
['maxprec'] = min(oldmaxprec
, 2*prec
)
572 xlow
= as_mpmath(xlow
, prec
+15, options
)
573 xhigh
= as_mpmath(xhigh
, prec
+15, options
)
575 # Integration is like summation, and we can phone home from
576 # the integrand function to update accuracy summation style
577 # Note that this accuracy is inaccurate, since it fails
578 # to account for the variable quadrature weights,
579 # but it is better than nothing
581 have_part
= [False, False]
582 max_real_term
= [MINUS_INF
]
583 max_imag_term
= [MINUS_INF
]
586 re
, im
, re_acc
, im_acc
= evalf(func
, prec
+15, {'subs':{x
:t
}})
588 have_part
[0] = re
or have_part
[0]
589 have_part
[1] = im
or have_part
[1]
591 max_real_term
[0] = max(max_real_term
[0], fastlog(re
))
592 max_imag_term
[0] = max(max_imag_term
[0], fastlog(im
))
595 return mpc(re
or fzero
, im
)
596 return mpf(re
or fzero
)
598 if options
.get('quad') == 'osc':
599 A
= C
.Wild('A', exclude
=[x
])
600 B
= C
.Wild('B', exclude
=[x
])
602 m
= func
.match(C
.cos(A
*x
+B
)*D
)
604 m
= func
.match(C
.sin(A
*x
+B
)*D
)
606 raise ValueError("An integrand of the form sin(A*x+B)*f(x) "
607 "or cos(A*x+B)*f(x) is required for oscillatory quadrature")
608 period
= as_mpmath(2*S
.Pi
/m
[A
], prec
+15, options
)
609 result
= quadosc(f
, xlow
, xhigh
, period
=period
)
610 # XXX: quadosc does not do error detection yet
611 quadrature_error
= MINUS_INF
613 result
, quadrature_error
= quadts(f
, xlow
, xhigh
, error
=1)
614 quadrature_error
= fastlog(quadrature_error
._mpf
_)
617 options
['maxprec'] = oldmaxprec
621 re
= result
.real
._mpf
_
623 re
= fshift(fone
, min(-prec
,-max_real_term
[0],-quadrature_error
))
626 re_acc
= -max(max_real_term
[0]-fastlog(re
)-prec
, quadrature_error
)
628 re
, re_acc
= None, None
631 im
= result
.imag
._mpf
_
633 im
= fshift(fone
, min(-prec
,-max_imag_term
[0],-quadrature_error
))
636 im_acc
= -max(max_imag_term
[0]-fastlog(im
)-prec
, quadrature_error
)
638 im
, im_acc
= None, None
640 result
= re
, im
, re_acc
, im_acc
643 def evalf_integral(expr
, prec
, options
):
646 maxprec
= options
.get('maxprec', INF
)
648 result
= do_integral(expr
, workprec
, options
)
649 accuracy
= complex_accuracy(result
)
650 if accuracy
>= prec
or workprec
>= maxprec
:
652 workprec
+= prec
- max(-2**i
, accuracy
)
655 def check_convergence(numer
, denom
, n
):
657 Returns (h, g, p) where
659 > 0 for convergence of rate 1/factorial(n)**h
660 < 0 for divergence of rate factorial(n)**(-h)
661 = 0 for geometric or polynomial convergence or divergence
664 > 1 for geometric convergence of rate 1/h**n
665 < 1 for geometric divergence of rate h**n
666 = 1 for polynomial convergence or divergence
668 (g < 0 indicates an alternating series)
671 > 1 for polynomial convergence of rate 1/n**h
672 <= 1 for polynomial divergence of rate n**(-h)
675 npol
= C
.Poly(numer
, n
)
676 dpol
= C
.Poly(denom
, n
)
681 return rate
, None, None
682 constant
= dpol
.lead_term
[0] / npol
.lead_term
[0]
683 if abs(constant
) != 1:
684 return rate
, constant
, None
685 if npol
.degree
== dpol
.degree
== 0:
686 return rate
, constant
, 0
687 pc
= list(npol
.iter_all_terms())[1][0]
688 qc
= list(dpol
.iter_all_terms())[1][0]
689 return rate
, constant
, qc
-pc
691 def hypsum(expr
, n
, start
, prec
):
693 Sum a rapidly convergent infinite hypergeometric series with
694 given general term, e.g. e = hypsum(1/factorial(n), n). The
695 quotient between successive terms must be a quotient of integer
698 from sympy
import hypersimp
, lambdify
701 expr
= expr
.subs(n
, n
+start
)
702 hs
= hypersimp(expr
, n
)
704 raise NotImplementedError("a hypergeometric series is required")
705 num
, den
= hs
.as_numer_denom()
706 func1
= lambdify(n
, num
)
707 func2
= lambdify(n
, den
)
709 h
, g
, p
= check_convergence(num
, den
, n
)
712 raise ValueError("Sum diverges like (n!)^%i" % (-h
))
714 # Direct summation if geometric or faster
715 if h
> 0 or (h
== 0 and abs(g
) > 1):
716 one
= MP_BASE(1) << prec
717 term
= expr
.subs(n
, 0)
718 term
= (MP_BASE(term
.p
) << prec
) // term
.q
722 term
*= MP_BASE(func1(k
-1))
723 term
//= MP_BASE(func2(k
-1))
726 return from_man_exp(s
, -prec
)
730 raise ValueError("Sum diverges like (%i)^n" % abs(1/g
))
731 if p
< 1 or (p
== 1 and not alt
):
732 raise ValueError("Sum diverges like n^%i" % (-p
))
733 # We have polynomial convergence:
734 # Use Shanks extrapolation for alternating series,
735 # Richardson extrapolation for nonalternating series
737 # XXX: better parameters for Shanks transformation
738 # This tends to get bad somewhere > 50 digits
739 N
= 5 + int(prec
*0.36)
743 N
= 3 + int(prec
*0.15)
746 # Need to use at least double precision because a lot of cancellation
747 # might occur in the extrapolation process
749 one
= MP_BASE(1) << prec2
750 term
= expr
.subs(n
, 0)
751 term
= (MP_BASE(term
.p
) << prec2
) // term
.q
753 table
= [make_mpf(from_man_exp(s
, -prec2
))]
754 for k
in xrange(1, NTERMS
):
755 term
*= MP_BASE(func1(k
-1))
756 term
//= MP_BASE(func2(k
-1))
758 table
.append(make_mpf(from_man_exp(s
, -prec2
)))
764 v
= shanks_extrapolation(table
, N
, M
)
766 v
= richardson_extrapolation(table
, N
, M
)
771 def evalf_sum(expr
, prec
, options
):
774 if len(limits
) != 1 or not isinstance(limits
[0], tuple) or \
776 raise NotImplementedError
780 if b
!= S
.Infinity
or a
!= int(a
):
781 raise NotImplementedError
782 # Use fast hypergeometric summation if possible
783 v
= hypsum(func
, n
, int(a
), prec2
)
784 delta
= prec
- fastlog(v
)
786 v
= hypsum(func
, n
, int(a
), delta
)
787 return v
, None, min(prec
, delta
), None
788 except NotImplementedError:
789 # Euler-Maclaurin summation for general series
790 eps
= C
.Real(2.0)**(-prec
)
791 for i
in range(1, 5):
793 s
, err
= expr
.euler_maclaurin(m
=m
, n
=n
, eps
=eps
, \
798 err
= fastlog(evalf(abs(err
), 20, options
)[0])
799 re
, im
, re_acc
, im_acc
= evalf(s
, prec2
, options
)
800 re_acc
= max(re_acc
, -err
)
801 im_acc
= max(im_acc
, -err
)
802 return re
, im
, re_acc
, im_acc
805 #----------------------------------------------------------------------------#
807 # Symbolic interface #
809 #----------------------------------------------------------------------------#
811 def evalf_symbol(x
, prec
, options
):
812 val
= options
['subs'][x
]
813 if isinstance(val
, mpf
):
815 return None, None, None, None
816 return val
._mpf
_, None, prec
, None
818 if not '_cache' in options
:
819 options
['_cache'] = {}
820 cache
= options
['_cache']
821 cached
, cached_prec
= cache
.get(x
.name
, (None, MINUS_INF
))
822 if cached_prec
>= prec
:
824 v
= evalf(sympify(val
), prec
, options
)
825 cache
[x
.name
] = (v
, prec
)
830 def _create_evalf_table():
833 C
.Symbol
: evalf_symbol
,
834 C
.Dummy
: evalf_symbol
,
835 C
.Real
: lambda x
, prec
, options
: (x
._mpf
_, None, prec
, None),
836 C
.Rational
: lambda x
, prec
, options
: (from_rational(x
.p
, x
.q
, prec
), None, prec
, None),
837 C
.Integer
: lambda x
, prec
, options
: (from_int(x
.p
, prec
), None, prec
, None),
838 C
.Zero
: lambda x
, prec
, options
: (None, None, prec
, None),
839 C
.One
: lambda x
, prec
, options
: (fone
, None, prec
, None),
840 C
.Half
: lambda x
, prec
, options
: (fhalf
, None, prec
, None),
841 C
.Pi
: lambda x
, prec
, options
: (fpi(prec
), None, prec
, None),
842 C
.Exp1
: lambda x
, prec
, options
: (fe(prec
), None, prec
, None),
843 C
.ImaginaryUnit
: lambda x
, prec
, options
: (None, fone
, None, prec
),
844 C
.NegativeOne
: lambda x
, prec
, options
: (fnone
, None, prec
, None),
846 C
.exp
: lambda x
, prec
, options
: evalf_pow(C
.Pow(S
.Exp1
, x
.args
[0],
847 evaluate
=False), prec
, options
),
862 C
.floor
: evalf_floor
,
863 C
.ceiling
: evalf_ceiling
,
865 C
.Integral
: evalf_integral
,
869 def evalf(x
, prec
, options
):
871 r
= evalf_table
[x
.func
](x
, prec
, options
)
873 #r = finalize_complex(x._eval_evalf(prec)._mpf_, fzero, prec)
875 # Fall back to ordinary evalf if possible
876 if 'subs' in options
:
877 x
= x
.subs(options
['subs'])
878 r
= x
._eval
_evalf
(prec
)._mpf
_, None, prec
, None
879 except AttributeError:
880 raise NotImplementedError
881 if options
.get("verbose"):
883 print "### output", to_str(r
[0] or fzero
, 50)
884 #print "### raw", r[0], r[2]
886 if options
.get("chop"):
887 r
= chop_parts(r
, prec
)
888 if options
.get("strict"):
889 check_target(x
, r
, prec
)
892 def Basic_evalf(x
, n
=15, **options
):
894 Evaluate the given formula to an accuracy of n digits.
895 Optional keyword arguments:
898 Substitute numerical values for symbols, e.g.
902 Allow a maximum temporary working precision of N digits
906 Replace tiny real or imaginary parts in subresults
907 by exact zeros (default=False)
910 Raise PrecisionExhausted if any subresult fails to evaluate
911 to full accuracy, given the available maxprec
915 Choose algorithm for numerical quadrature. By default,
916 tanh-sinh quadrature is used. For oscillatory
917 integrals on an infinite interval, try quad='osc'.
920 Print debug information (default=False)
924 _create_evalf_table()
925 prec
= dps_to_prec(n
)
926 if 'maxprec' in options
:
927 options
['maxprec'] = int(options
['maxprec']*LG10
)
929 options
['maxprec'] = max(prec
, DEFAULT_MAXPREC
)
931 result
= evalf(x
, prec
+4, options
)
932 except NotImplementedError:
933 # Fall back to the ordinary evalf
934 v
= x
._eval
_evalf
(prec
)
938 # If the result is numerical, normalize it
939 result
= evalf(v
, prec
, options
)
941 # Probably contains symbols or unknown functions
943 re
, im
, re_acc
, im_acc
= result
945 p
= max(min(prec
, re_acc
), 1)
946 #re = fpos(re, p, round_nearest)
947 re
= C
.Real
._new
(re
, p
)
951 p
= max(min(prec
, im_acc
), 1)
952 #im = fpos(im, p, round_nearest)
953 im
= C
.Real
._new
(im
, p
)
954 return re
+ im
*S
.ImaginaryUnit
958 Basic
.evalf
= Basic
.n
= Basic_evalf
960 def N(x
, n
=15, **options
):
962 Calls x.evalf(n, **options).
964 Both .evalf() and N() are equivalent, use the one that you like better.
967 >>> from sympy import Sum, Symbol, oo
969 >>> Sum(1/k**k, (k, 1, oo))
970 Sum(k**(-k), (k, 1, oo))
971 >>> N(Sum(1/k**k, (k, 1, oo)), 4)
975 return sympify(x
).evalf(n
, **options
)