Coverage for sympy/core/evalf.py : 48%

""" 

Adaptive numerical evaluation of SymPy expressions, using mpmath 

for mathematical functions. 

""" 

from __future__ import print_function, division 

import math 

import mpmath.libmp as libmp 

from mpmath import ( 

make_mpc, make_mpf, mp, mpc, mpf, nsum, quadts, quadosc, workprec) 

from mpmath import inf as mpmath_inf 

from mpmath.libmp import (from_int, from_man_exp, from_rational, fhalf, 

fnan, fnone, fone, fzero, mpf_abs, mpf_add, 

mpf_atan, mpf_atan2, mpf_cmp, mpf_cos, mpf_e, mpf_exp, mpf_log, mpf_lt, 

mpf_mul, mpf_neg, mpf_pi, mpf_pow, mpf_pow_int, mpf_shift, mpf_sin, 

mpf_sqrt, normalize, round_nearest, to_int, to_str) 

from mpmath.libmp import bitcount as mpmath_bitcount 

from mpmath.libmp.backend import MPZ 

from mpmath.libmp.libmpc import _infs_nan 

from mpmath.libmp.libmpf import dps_to_prec, prec_to_dps 

from mpmath.libmp.gammazeta import mpf_bernoulli 

from .compatibility import SYMPY_INTS, range 

from .sympify import sympify 

from .singleton import S 

from sympy.utilities.iterables import is_sequence 

LG10 = math.log(10, 2) 

rnd = round_nearest 

def bitcount(n): 

return mpmath_bitcount(int(n)) 

# Used in a few places as placeholder values to denote exponents and 

# precision levels, e.g. of exact numbers. Must be careful to avoid 

# passing these to mpmath functions or returning them in final results. 

INF = float(mpmath_inf) 

MINUS_INF = float(-mpmath_inf) 

# ~= 100 digits. Real men set this to INF. 

DEFAULT_MAXPREC = 333 

class PrecisionExhausted(ArithmeticError): 

pass 

#----------------------------------------------------------------------------# 

# # 

# Helper functions for arithmetic and complex parts # 

# # 

#----------------------------------------------------------------------------# 

""" 

An mpf value tuple is a tuple of integers (sign, man, exp, bc) 

representing a floating-point number: [1, -1][sign]*man*2**exp where 

sign is 0 or 1 and bc should correspond to the number of bits used to 

represent the mantissa (man) in binary notation, e.g. 

>>> from sympy.core.evalf import bitcount 

>>> sign, man, exp, bc = 0, 5, 1, 3 

>>> n = [1, -1][sign]*man*2**exp 

>>> n, bitcount(man) 

(10, 3) 

A temporary result is a tuple (re, im, re_acc, im_acc) where 

re and im are nonzero mpf value tuples representing approximate 

numbers, or None to denote exact zeros. 

re_acc, im_acc are integers denoting log2(e) where e is the estimated 

relative accuracy of the respective complex part, but may be anything 

if the corresponding complex part is None. 

""" 

def fastlog(x): 

"""Fast approximation of log2(x) for an mpf value tuple x. 

Notes: Calculated as exponent + width of mantissa. This is an 

approximation for two reasons: 1) it gives the ceil(log2(abs(x))) 

value and 2) it is too high by 1 in the case that x is an exact 

power of 2. Although this is easy to remedy by testing to see if 

the odd mpf mantissa is 1 (indicating that one was dealing with 

an exact power of 2) that would decrease the speed and is not 

necessary as this is only being used as an approximation for the 

number of bits in x. The correct return value could be written as 

"x[2] + (x[3] if x[1] != 1 else 0)". 

Since mpf tuples always have an odd mantissa, no check is done 

to see if the mantissa is a multiple of 2 (in which case the 

result would be too large by 1). 

Examples 

======== 

>>> from sympy import log 

>>> from sympy.core.evalf import fastlog, bitcount 

>>> s, m, e = 0, 5, 1 

>>> bc = bitcount(m) 

>>> n = [1, -1][s]*m*2**e 

>>> n, (log(n)/log(2)).evalf(2), fastlog((s, m, e, bc)) 

(10, 3.3, 4) 

""" 

if not x or x == fzero: 

return MINUS_INF 

return x[2] + x[3] 

def pure_complex(v, or_real=False): 

"""Return a and b if v matches a + I*b where b is not zero and 

a and b are Numbers, else None. If `or_real` is True then 0 will 

be returned for `b` if `v` is a real number. 

>>> from sympy.core.evalf import pure_complex 

>>> from sympy import sqrt, I, S 

>>> a, b, surd = S(2), S(3), sqrt(2) 

>>> pure_complex(a) 

>>> pure_complex(a, or_real=True) 

(2, 0) 

>>> pure_complex(surd) 

>>> pure_complex(a + b*I) 

(2, 3) 

>>> pure_complex(I) 

(0, 1) 

""" 

h, t = v.as_coeff_Add() 

if not t: 

if or_real: 

return h, t 

return 

c, i = t.as_coeff_Mul() 

if i is S.ImaginaryUnit: 

return h, c 

def scaled_zero(mag, sign=1): 

"""Return an mpf representing a power of two with magnitude ``mag`` 

and -1 for precision. Or, if ``mag`` is a scaled_zero tuple, then just 

remove the sign from within the list that it was initially wrapped 

in. 

Examples 

======== 

>>> from sympy.core.evalf import scaled_zero 

>>> from sympy import Float 

>>> z, p = scaled_zero(100) 

>>> z, p 

(([0], 1, 100, 1), -1) 

>>> ok = scaled_zero(z) 

>>> ok 

(0, 1, 100, 1) 

>>> Float(ok) 

1.26765060022823e+30 

>>> Float(ok, p) 

0.e+30 

>>> ok, p = scaled_zero(100, -1) 

>>> Float(scaled_zero(ok), p) 

-0.e+30 

""" 

if type(mag) is tuple and len(mag) == 4 and iszero(mag, scaled=True): 

return (mag[0][0],) + mag[1:] 

elif isinstance(mag, SYMPY_INTS): 

if sign not in [-1, 1]: 

raise ValueError('sign must be +/-1') 

rv, p = mpf_shift(fone, mag), -1 

s = 0 if sign == 1 else 1 

rv = ([s],) + rv[1:] 

return rv, p 

else: 

raise ValueError('scaled zero expects int or scaled_zero tuple.') 

def iszero(mpf, scaled=False): 

if not scaled: 

return not mpf or not mpf[1] and not mpf[-1] 

return mpf and type(mpf[0]) is list and mpf[1] == mpf[-1] == 1 

def complex_accuracy(result): 

""" 

Returns relative accuracy of a complex number with given accuracies 

for the real and imaginary parts. The relative accuracy is defined 

in the complex norm sense as ||z|+|error|| / |z| where error 

is equal to (real absolute error) + (imag absolute error)*i. 

The full expression for the (logarithmic) error can be approximated 

easily by using the max norm to approximate the complex norm. 

In the worst case (re and im equal), this is wrong by a factor 

sqrt(2), or by log2(sqrt(2)) = 0.5 bit. 

""" 

re, im, re_acc, im_acc = result 

if not im: 

if not re: 

return INF 

return re_acc 

if not re: 

return im_acc 

re_size = fastlog(re) 

im_size = fastlog(im) 

absolute_error = max(re_size - re_acc, im_size - im_acc) 

relative_error = absolute_error - max(re_size, im_size) 

return -relative_error 

def get_abs(expr, prec, options): 

re, im, re_acc, im_acc = evalf(expr, prec + 2, options) 

if not re: 

re, re_acc, im, im_acc = im, im_acc, re, re_acc 

if im: 

return libmp.mpc_abs((re, im), prec), None, re_acc, None 

elif re: 

return mpf_abs(re), None, re_acc, None 

else: 

return None, None, None, None 

def get_complex_part(expr, no, prec, options): 

"""no = 0 for real part, no = 1 for imaginary part""" 

workprec = prec 

i = 0 

while 1: 

res = evalf(expr, workprec, options) 

value, accuracy = res[no::2] 

# XXX is the last one correct? Consider re((1+I)**2).n() 

if (not value) or accuracy >= prec or -value[2] > prec: 

return value, None, accuracy, None 

workprec += max(30, 2**i) 

i += 1 

def evalf_abs(expr, prec, options): 

return get_abs(expr.args[0], prec, options) 

def evalf_re(expr, prec, options): 

return get_complex_part(expr.args[0], 0, prec, options) 

def evalf_im(expr, prec, options): 

return get_complex_part(expr.args[0], 1, prec, options) 

def finalize_complex(re, im, prec): 

if re == fzero and im == fzero: 

raise ValueError("got complex zero with unknown accuracy") 

elif re == fzero: 

return None, im, None, prec 

elif im == fzero: 

return re, None, prec, None 

size_re = fastlog(re) 

size_im = fastlog(im) 

if size_re > size_im: 

re_acc = prec 

im_acc = prec + min(-(size_re - size_im), 0) 

else: 

im_acc = prec 

re_acc = prec + min(-(size_im - size_re), 0) 

return re, im, re_acc, im_acc 

def chop_parts(value, prec): 

""" 

Chop off tiny real or complex parts. 

""" 

re, im, re_acc, im_acc = value 

# Method 1: chop based on absolute value 

if re and re not in _infs_nan and (fastlog(re) < -prec + 4): 

re, re_acc = None, None 

if im and im not in _infs_nan and (fastlog(im) < -prec + 4): 

im, im_acc = None, None 

# Method 2: chop if inaccurate and relatively small 

if re and im: 

delta = fastlog(re) - fastlog(im) 

if re_acc < 2 and (delta - re_acc <= -prec + 4): 

re, re_acc = None, None 

if im_acc < 2 and (delta - im_acc >= prec - 4): 

im, im_acc = None, None 

return re, im, re_acc, im_acc 

def check_target(expr, result, prec): 

a = complex_accuracy(result) 

if a < prec: 

raise PrecisionExhausted("Failed to distinguish the expression: \n\n%s\n\n" 

"from zero. Try simplifying the input, using chop=True, or providing " 

"a higher maxn for evalf" % (expr)) 

def get_integer_part(expr, no, options, return_ints=False): 

""" 

With no = 1, computes ceiling(expr) 

With no = -1, computes floor(expr) 

Note: this function either gives the exact result or signals failure. 

""" 

from sympy.functions.elementary.complexes import re, im 

# The expression is likely less than 2^30 or so 

assumed_size = 30 

ire, iim, ire_acc, iim_acc = evalf(expr, assumed_size, options) 

# We now know the size, so we can calculate how much extra precision 

# (if any) is needed to get within the nearest integer 

if ire and iim: 

gap = max(fastlog(ire) - ire_acc, fastlog(iim) - iim_acc) 

elif ire: 

gap = fastlog(ire) - ire_acc 

elif iim: 

gap = fastlog(iim) - iim_acc 

else: 

# ... or maybe the expression was exactly zero 

return None, None, None, None 

margin = 10 

if gap >= -margin: 

ire, iim, ire_acc, iim_acc = \ 

evalf(expr, margin + assumed_size + gap, options) 

# We can now easily find the nearest integer, but to find floor/ceil, we 

# must also calculate whether the difference to the nearest integer is 

# positive or negative (which may fail if very close). 

def calc_part(expr, nexpr): 

from sympy.core.add import Add 

nint = int(to_int(nexpr, rnd)) 

n, c, p, b = nexpr 

is_int = (p == 0) 

if not is_int: 

# if there are subs and they all contain integer re/im parts 

# then we can (hopefully) safely substitute them into the 

# expression 

s = options.get('subs', False) 

if s: 

doit = True 

from sympy.core.compatibility import as_int 

for v in s.values(): 

try: 

as_int(v) 

except ValueError: 

try: 

[as_int(i) for i in v.as_real_imag()] 

continue 

except (ValueError, AttributeError): 

doit = False 

break 

if doit: 

expr = expr.subs(s) 

expr = Add(expr, -nint, evaluate=False) 

x, _, x_acc, _ = evalf(expr, 10, options) 

try: 

check_target(expr, (x, None, x_acc, None), 3) 

except PrecisionExhausted: 

if not expr.equals(0): 

raise PrecisionExhausted 

x = fzero 

nint += int(no*(mpf_cmp(x or fzero, fzero) == no)) 

nint = from_int(nint) 

return nint, fastlog(nint) + 10 

re_, im_, re_acc, im_acc = None, None, None, None 

if ire: 

re_, re_acc = calc_part(re(expr, evaluate=False), ire) 

if iim: 

im_, im_acc = calc_part(im(expr, evaluate=False), iim) 

if return_ints: 

return int(to_int(re_ or fzero)), int(to_int(im_ or fzero)) 

return re_, im_, re_acc, im_acc 

def evalf_ceiling(expr, prec, options): 

return get_integer_part(expr.args[0], 1, options) 

def evalf_floor(expr, prec, options): 

return get_integer_part(expr.args[0], -1, options) 

#----------------------------------------------------------------------------# 

# # 

# Arithmetic operations # 

# # 

#----------------------------------------------------------------------------# 

def add_terms(terms, prec, target_prec): 

""" 

Helper for evalf_add. Adds a list of (mpfval, accuracy) terms. 

Returns 

------- 

- None, None if there are no non-zero terms; 

- terms[0] if there is only 1 term; 

- scaled_zero if the sum of the terms produces a zero by cancellation 

e.g. mpfs representing 1 and -1 would produce a scaled zero which need 

special handling since they are not actually zero and they are purposely 

malformed to ensure that they can't be used in anything but accuracy 

calculations; 

- a tuple that is scaled to target_prec that corresponds to the 

sum of the terms. 

The returned mpf tuple will be normalized to target_prec; the input 

prec is used to define the working precision. 

XXX explain why this is needed and why one can't just loop using mpf_add 

""" 

terms = [t for t in terms if not iszero(t)] 

if not terms: 

return None, None 

elif len(terms) == 1: 

return terms[0] 

# see if any argument is NaN or oo and thus warrants a special return 

special = [] 

from sympy.core.numbers import Float 

for t in terms: 

arg = Float._new(t[0], 1) 

if arg is S.NaN or arg.is_infinite: 

special.append(arg) 

if special: 

from sympy.core.add import Add 

rv = evalf(Add(*special), prec + 4, {}) 

return rv[0], rv[2] 

working_prec = 2*prec 

sum_man, sum_exp, absolute_error = 0, 0, MINUS_INF 

for x, accuracy in terms: 

sign, man, exp, bc = x 

if sign: 

man = -man 

absolute_error = max(absolute_error, bc + exp - accuracy) 

delta = exp - sum_exp 

if exp >= sum_exp: 

# x much larger than existing sum? 

# first: quick test 

if ((delta > working_prec) and 

((not sum_man) or 

delta - bitcount(abs(sum_man)) > working_prec)): 

sum_man = man 

sum_exp = exp 

else: 

sum_man += (man << delta) 

else: 

delta = -delta 

# x much smaller than existing sum? 

if delta - bc > working_prec: 

if not sum_man: 

sum_man, sum_exp = man, exp 

else: 

sum_man = (sum_man << delta) + man 

sum_exp = exp 

if not sum_man: 

return scaled_zero(absolute_error) 

if sum_man < 0: 

sum_sign = 1 

sum_man = -sum_man 

else: 

sum_sign = 0 

sum_bc = bitcount(sum_man) 

sum_accuracy = sum_exp + sum_bc - absolute_error 

r = normalize(sum_sign, sum_man, sum_exp, sum_bc, target_prec, 

rnd), sum_accuracy 

return r 

def evalf_add(v, prec, options): 

res = pure_complex(v) 

if res: 

h, c = res 

re, _, re_acc, _ = evalf(h, prec, options) 

im, _, im_acc, _ = evalf(c, prec, options) 

return re, im, re_acc, im_acc 

oldmaxprec = options.get('maxprec', DEFAULT_MAXPREC) 

i = 0 

target_prec = prec 

while 1: 

options['maxprec'] = min(oldmaxprec, 2*prec) 

terms = [evalf(arg, prec + 10, options) for arg in v.args] 

re, re_acc = add_terms( 

[a[0::2] for a in terms if a[0]], prec, target_prec) 

im, im_acc = add_terms( 

[a[1::2] for a in terms if a[1]], prec, target_prec) 

acc = complex_accuracy((re, im, re_acc, im_acc)) 

if acc >= target_prec: 

if options.get('verbose'): 

print("ADD: wanted", target_prec, "accurate bits, got", re_acc, im_acc) 

break 

else: 

if (prec - target_prec) > options['maxprec']: 

break 

prec = prec + max(10 + 2**i, target_prec - acc) 

i += 1 

if options.get('verbose'): 

print("ADD: restarting with prec", prec) 

options['maxprec'] = oldmaxprec 

if iszero(re, scaled=True): 

re = scaled_zero(re) 

if iszero(im, scaled=True): 

im = scaled_zero(im) 

return re, im, re_acc, im_acc 

def evalf_mul(v, prec, options): 

res = pure_complex(v) 

if res: 

# the only pure complex that is a mul is h*I 

_, h = res 

im, _, im_acc, _ = evalf(h, prec, options) 

return None, im, None, im_acc 

args = list(v.args) 

# see if any argument is NaN or oo and thus warrants a special return 

special = [] 

from sympy.core.numbers import Float 

for arg in args: 

arg = evalf(arg, prec, options) 

if arg[0] is None: 

continue 

arg = Float._new(arg[0], 1) 

if arg is S.NaN or arg.is_infinite: 

special.append(arg) 

if special: 

from sympy.core.mul import Mul 

special = Mul(*special) 

return evalf(special, prec + 4, {}) 

# With guard digits, multiplication in the real case does not destroy 

# accuracy. This is also true in the complex case when considering the 

# total accuracy; however accuracy for the real or imaginary parts 

# separately may be lower. 

acc = prec 

# XXX: big overestimate 

working_prec = prec + len(args) + 5 

# Empty product is 1 

start = man, exp, bc = MPZ(1), 0, 1 

# First, we multiply all pure real or pure imaginary numbers. 

# direction tells us that the result should be multiplied by 

# I**direction; all other numbers get put into complex_factors 

# to be multiplied out after the first phase. 

last = len(args) 

direction = 0 

args.append(S.One) 

complex_factors = [] 

for i, arg in enumerate(args): 

if i != last and pure_complex(arg): 

args[-1] = (args[-1]*arg).expand() 

continue 

elif i == last and arg is S.One: 

continue 

re, im, re_acc, im_acc = evalf(arg, working_prec, options) 

if re and im: 

complex_factors.append((re, im, re_acc, im_acc)) 

continue 

elif re: 

(s, m, e, b), w_acc = re, re_acc 

elif im: 

(s, m, e, b), w_acc = im, im_acc 

direction += 1 

else: 

return None, None, None, None 

direction += 2*s 

man *= m 

exp += e 

bc += b 

if bc > 3*working_prec: 

man >>= working_prec 

exp += working_prec 

acc = min(acc, w_acc) 

sign = (direction & 2) >> 1 

if not complex_factors: 

v = normalize(sign, man, exp, bitcount(man), prec, rnd) 

# multiply by i 

if direction & 1: 

return None, v, None, acc 

else: 

return v, None, acc, None 

else: 

# initialize with the first term 

if (man, exp, bc) != start: 

# there was a real part; give it an imaginary part 

re, im = (sign, man, exp, bitcount(man)), (0, MPZ(0), 0, 0) 

i0 = 0 

else: 

# there is no real part to start (other than the starting 1) 

wre, wim, wre_acc, wim_acc = complex_factors[0] 

acc = min(acc, 

complex_accuracy((wre, wim, wre_acc, wim_acc))) 

re = wre 

im = wim 

i0 = 1 

for wre, wim, wre_acc, wim_acc in complex_factors[i0:]: 

# acc is the overall accuracy of the product; we aren't 

# computing exact accuracies of the product. 

acc = min(acc, 

complex_accuracy((wre, wim, wre_acc, wim_acc))) 

use_prec = working_prec 

A = mpf_mul(re, wre, use_prec) 

B = mpf_mul(mpf_neg(im), wim, use_prec) 

C = mpf_mul(re, wim, use_prec) 

D = mpf_mul(im, wre, use_prec) 

re = mpf_add(A, B, use_prec) 

im = mpf_add(C, D, use_prec) 

if options.get('verbose'): 

print("MUL: wanted", prec, "accurate bits, got", acc) 

# multiply by I 

if direction & 1: 

re, im = mpf_neg(im), re 

return re, im, acc, acc 

def evalf_pow(v, prec, options): 

target_prec = prec 

base, exp = v.args 

# We handle x**n separately. This has two purposes: 1) it is much 

# faster, because we avoid calling evalf on the exponent, and 2) it 

# allows better handling of real/imaginary parts that are exactly zero 

if exp.is_Integer: 

p = exp.p 

# Exact 

if not p: 

return fone, None, prec, None 

# Exponentiation by p magnifies relative error by |p|, so the 

# base must be evaluated with increased precision if p is large 

prec += int(math.log(abs(p), 2)) 

re, im, re_acc, im_acc = evalf(base, prec + 5, options) 

# Real to integer power 

if re and not im: 

return mpf_pow_int(re, p, target_prec), None, target_prec, None 

# (x*I)**n = I**n * x**n 

if im and not re: 

z = mpf_pow_int(im, p, target_prec) 

case = p % 4 

if case == 0: 

return z, None, target_prec, None 

if case == 1: 

return None, z, None, target_prec 

if case == 2: 

return mpf_neg(z), None, target_prec, None 

if case == 3: 

return None, mpf_neg(z), None, target_prec 

# Zero raised to an integer power 

if not re: 

return None, None, None, None 

# General complex number to arbitrary integer power 

re, im = libmp.mpc_pow_int((re, im), p, prec) 

# Assumes full accuracy in input 

return finalize_complex(re, im, target_prec) 

# Pure square root 

if exp is S.Half: 

xre, xim, _, _ = evalf(base, prec + 5, options) 

# General complex square root 

if xim: 

re, im = libmp.mpc_sqrt((xre or fzero, xim), prec) 

return finalize_complex(re, im, prec) 

if not xre: 

return None, None, None, None 

# Square root of a negative real number 

if mpf_lt(xre, fzero): 

return None, mpf_sqrt(mpf_neg(xre), prec), None, prec 

# Positive square root 

return mpf_sqrt(xre, prec), None, prec, None 

# We first evaluate the exponent to find its magnitude 

# This determines the working precision that must be used 

prec += 10 

yre, yim, _, _ = evalf(exp, prec, options) 

# Special cases: x**0 

if not (yre or yim): 

return fone, None, prec, None 

ysize = fastlog(yre) 

# Restart if too big 

# XXX: prec + ysize might exceed maxprec 

if ysize > 5: 

prec += ysize 

yre, yim, _, _ = evalf(exp, prec, options) 

# Pure exponential function; no need to evalf the base 

if base is S.Exp1: 

if yim: 

re, im = libmp.mpc_exp((yre or fzero, yim), prec) 

return finalize_complex(re, im, target_prec) 

return mpf_exp(yre, target_prec), None, target_prec, None 

xre, xim, _, _ = evalf(base, prec + 5, options) 

# 0**y 

if not (xre or xim): 

return None, None, None, None 

# (real ** complex) or (complex ** complex) 

if yim: 

re, im = libmp.mpc_pow( 

(xre or fzero, xim or fzero), (yre or fzero, yim), 

target_prec) 

return finalize_complex(re, im, target_prec)