""" This module defines the mpf, mpc classes, and standard functions for operating with them. """ __docformat__ = 'plaintext' import functools import re from .ctx_base import StandardBaseContext from .libmp.backend import basestring, BACKEND from . import libmp from .libmp import (MPZ, MPZ_ZERO, MPZ_ONE, int_types, repr_dps, round_floor, round_ceiling, dps_to_prec, round_nearest, prec_to_dps, ComplexResult, to_pickable, from_pickable, normalize, from_int, from_float, from_str, to_int, to_float, to_str, from_rational, from_man_exp, fone, fzero, finf, fninf, fnan, mpf_abs, mpf_pos, mpf_neg, mpf_add, mpf_sub, mpf_mul, mpf_mul_int, mpf_div, mpf_rdiv_int, mpf_pow_int, mpf_mod, mpf_eq, mpf_cmp, mpf_lt, mpf_gt, mpf_le, mpf_ge, mpf_hash, mpf_rand, mpf_sum, bitcount, to_fixed, mpc_to_str, mpc_to_complex, mpc_hash, mpc_pos, mpc_is_nonzero, mpc_neg, mpc_conjugate, mpc_abs, mpc_add, mpc_add_mpf, mpc_sub, mpc_sub_mpf, mpc_mul, mpc_mul_mpf, mpc_mul_int, mpc_div, mpc_div_mpf, mpc_pow, mpc_pow_mpf, mpc_pow_int, mpc_mpf_div, mpf_pow, mpf_pi, mpf_degree, mpf_e, mpf_phi, mpf_ln2, mpf_ln10, mpf_euler, mpf_catalan, mpf_apery, mpf_khinchin, mpf_glaisher, mpf_twinprime, mpf_mertens, int_types) from . import function_docs from . import rational new = object.__new__ get_complex = re.compile(r'^\(?(?P[\+\-]?\d*\.?\d*(e[\+\-]?\d+)?)??' r'(?P[\+\-]?\d*\.?\d*(e[\+\-]?\d+)?j)?\)?$') if BACKEND == 'sage': from sage.libs.mpmath.ext_main import Context as BaseMPContext # pickle hack import sage.libs.mpmath.ext_main as _mpf_module else: from .ctx_mp_python import PythonMPContext as BaseMPContext from . import ctx_mp_python as _mpf_module from .ctx_mp_python import _mpf, _mpc, mpnumeric class MPContext(BaseMPContext, StandardBaseContext): """ Context for multiprecision arithmetic with a global precision. """ def __init__(ctx): BaseMPContext.__init__(ctx) ctx.trap_complex = False ctx.pretty = False ctx.types = [ctx.mpf, ctx.mpc, ctx.constant] ctx._mpq = rational.mpq ctx.default() StandardBaseContext.__init__(ctx) ctx.mpq = rational.mpq ctx.init_builtins() ctx.hyp_summators = {} ctx._init_aliases() # XXX: automate try: ctx.bernoulli.im_func.func_doc = function_docs.bernoulli ctx.primepi.im_func.func_doc = function_docs.primepi ctx.psi.im_func.func_doc = function_docs.psi ctx.atan2.im_func.func_doc = function_docs.atan2 except AttributeError: # python 3 ctx.bernoulli.__func__.func_doc = function_docs.bernoulli ctx.primepi.__func__.func_doc = function_docs.primepi ctx.psi.__func__.func_doc = function_docs.psi ctx.atan2.__func__.func_doc = function_docs.atan2 ctx.digamma.func_doc = function_docs.digamma ctx.cospi.func_doc = function_docs.cospi ctx.sinpi.func_doc = function_docs.sinpi def init_builtins(ctx): mpf = ctx.mpf mpc = ctx.mpc # Exact constants ctx.one = ctx.make_mpf(fone) ctx.zero = ctx.make_mpf(fzero) ctx.j = ctx.make_mpc((fzero,fone)) ctx.inf = ctx.make_mpf(finf) ctx.ninf = ctx.make_mpf(fninf) ctx.nan = ctx.make_mpf(fnan) eps = ctx.constant(lambda prec, rnd: (0, MPZ_ONE, 1-prec, 1), "epsilon of working precision", "eps") ctx.eps = eps # Approximate constants ctx.pi = ctx.constant(mpf_pi, "pi", "pi") ctx.ln2 = ctx.constant(mpf_ln2, "ln(2)", "ln2") ctx.ln10 = ctx.constant(mpf_ln10, "ln(10)", "ln10") ctx.phi = ctx.constant(mpf_phi, "Golden ratio phi", "phi") ctx.e = ctx.constant(mpf_e, "e = exp(1)", "e") ctx.euler = ctx.constant(mpf_euler, "Euler's constant", "euler") ctx.catalan = ctx.constant(mpf_catalan, "Catalan's constant", "catalan") ctx.khinchin = ctx.constant(mpf_khinchin, "Khinchin's constant", "khinchin") ctx.glaisher = ctx.constant(mpf_glaisher, "Glaisher's constant", "glaisher") ctx.apery = ctx.constant(mpf_apery, "Apery's constant", "apery") ctx.degree = ctx.constant(mpf_degree, "1 deg = pi / 180", "degree") ctx.twinprime = ctx.constant(mpf_twinprime, "Twin prime constant", "twinprime") ctx.mertens = ctx.constant(mpf_mertens, "Mertens' constant", "mertens") # Standard functions ctx.sqrt = ctx._wrap_libmp_function(libmp.mpf_sqrt, libmp.mpc_sqrt) ctx.cbrt = ctx._wrap_libmp_function(libmp.mpf_cbrt, libmp.mpc_cbrt) ctx.ln = ctx._wrap_libmp_function(libmp.mpf_log, libmp.mpc_log) ctx.atan = ctx._wrap_libmp_function(libmp.mpf_atan, libmp.mpc_atan) ctx.exp = ctx._wrap_libmp_function(libmp.mpf_exp, libmp.mpc_exp) ctx.expj = ctx._wrap_libmp_function(libmp.mpf_expj, libmp.mpc_expj) ctx.expjpi = ctx._wrap_libmp_function(libmp.mpf_expjpi, libmp.mpc_expjpi) ctx.sin = ctx._wrap_libmp_function(libmp.mpf_sin, libmp.mpc_sin) ctx.cos = ctx._wrap_libmp_function(libmp.mpf_cos, libmp.mpc_cos) ctx.tan = ctx._wrap_libmp_function(libmp.mpf_tan, libmp.mpc_tan) ctx.sinh = ctx._wrap_libmp_function(libmp.mpf_sinh, libmp.mpc_sinh) ctx.cosh = ctx._wrap_libmp_function(libmp.mpf_cosh, libmp.mpc_cosh) ctx.tanh = ctx._wrap_libmp_function(libmp.mpf_tanh, libmp.mpc_tanh) ctx.asin = ctx._wrap_libmp_function(libmp.mpf_asin, libmp.mpc_asin) ctx.acos = ctx._wrap_libmp_function(libmp.mpf_acos, libmp.mpc_acos) ctx.atan = ctx._wrap_libmp_function(libmp.mpf_atan, libmp.mpc_atan) ctx.asinh = ctx._wrap_libmp_function(libmp.mpf_asinh, libmp.mpc_asinh) ctx.acosh = ctx._wrap_libmp_function(libmp.mpf_acosh, libmp.mpc_acosh) ctx.atanh = ctx._wrap_libmp_function(libmp.mpf_atanh, libmp.mpc_atanh) ctx.sinpi = ctx._wrap_libmp_function(libmp.mpf_sin_pi, libmp.mpc_sin_pi) ctx.cospi = ctx._wrap_libmp_function(libmp.mpf_cos_pi, libmp.mpc_cos_pi) ctx.floor = ctx._wrap_libmp_function(libmp.mpf_floor, libmp.mpc_floor) ctx.ceil = ctx._wrap_libmp_function(libmp.mpf_ceil, libmp.mpc_ceil) ctx.nint = ctx._wrap_libmp_function(libmp.mpf_nint, libmp.mpc_nint) ctx.frac = ctx._wrap_libmp_function(libmp.mpf_frac, libmp.mpc_frac) ctx.fib = ctx.fibonacci = ctx._wrap_libmp_function(libmp.mpf_fibonacci, libmp.mpc_fibonacci) ctx.gamma = ctx._wrap_libmp_function(libmp.mpf_gamma, libmp.mpc_gamma) ctx.rgamma = ctx._wrap_libmp_function(libmp.mpf_rgamma, libmp.mpc_rgamma) ctx.loggamma = ctx._wrap_libmp_function(libmp.mpf_loggamma, libmp.mpc_loggamma) ctx.fac = ctx.factorial = ctx._wrap_libmp_function(libmp.mpf_factorial, libmp.mpc_factorial) ctx.digamma = ctx._wrap_libmp_function(libmp.mpf_psi0, libmp.mpc_psi0) ctx.harmonic = ctx._wrap_libmp_function(libmp.mpf_harmonic, libmp.mpc_harmonic) ctx.ei = ctx._wrap_libmp_function(libmp.mpf_ei, libmp.mpc_ei) ctx.e1 = ctx._wrap_libmp_function(libmp.mpf_e1, libmp.mpc_e1) ctx._ci = ctx._wrap_libmp_function(libmp.mpf_ci, libmp.mpc_ci) ctx._si = ctx._wrap_libmp_function(libmp.mpf_si, libmp.mpc_si) ctx.ellipk = ctx._wrap_libmp_function(libmp.mpf_ellipk, libmp.mpc_ellipk) ctx._ellipe = ctx._wrap_libmp_function(libmp.mpf_ellipe, libmp.mpc_ellipe) ctx.agm1 = ctx._wrap_libmp_function(libmp.mpf_agm1, libmp.mpc_agm1) ctx._erf = ctx._wrap_libmp_function(libmp.mpf_erf, None) ctx._erfc = ctx._wrap_libmp_function(libmp.mpf_erfc, None) ctx._zeta = ctx._wrap_libmp_function(libmp.mpf_zeta, libmp.mpc_zeta) ctx._altzeta = ctx._wrap_libmp_function(libmp.mpf_altzeta, libmp.mpc_altzeta) # Faster versions ctx.sqrt = getattr(ctx, "_sage_sqrt", ctx.sqrt) ctx.exp = getattr(ctx, "_sage_exp", ctx.exp) ctx.ln = getattr(ctx, "_sage_ln", ctx.ln) ctx.cos = getattr(ctx, "_sage_cos", ctx.cos) ctx.sin = getattr(ctx, "_sage_sin", ctx.sin) def to_fixed(ctx, x, prec): return x.to_fixed(prec) def hypot(ctx, x, y): r""" Computes the Euclidean norm of the vector `(x, y)`, equal to `\sqrt{x^2 + y^2}`. Both `x` and `y` must be real.""" x = ctx.convert(x) y = ctx.convert(y) return ctx.make_mpf(libmp.mpf_hypot(x._mpf_, y._mpf_, *ctx._prec_rounding)) def _gamma_upper_int(ctx, n, z): n = int(ctx._re(n)) if n == 0: return ctx.e1(z) if not hasattr(z, '_mpf_'): raise NotImplementedError prec, rounding = ctx._prec_rounding real, imag = libmp.mpf_expint(n, z._mpf_, prec, rounding, gamma=True) if imag is None: return ctx.make_mpf(real) else: return ctx.make_mpc((real, imag)) def _expint_int(ctx, n, z): n = int(n) if n == 1: return ctx.e1(z) if not hasattr(z, '_mpf_'): raise NotImplementedError prec, rounding = ctx._prec_rounding real, imag = libmp.mpf_expint(n, z._mpf_, prec, rounding) if imag is None: return ctx.make_mpf(real) else: return ctx.make_mpc((real, imag)) def _nthroot(ctx, x, n): if hasattr(x, '_mpf_'): try: return ctx.make_mpf(libmp.mpf_nthroot(x._mpf_, n, *ctx._prec_rounding)) except ComplexResult: if ctx.trap_complex: raise x = (x._mpf_, libmp.fzero) else: x = x._mpc_ return ctx.make_mpc(libmp.mpc_nthroot(x, n, *ctx._prec_rounding)) def _besselj(ctx, n, z): prec, rounding = ctx._prec_rounding if hasattr(z, '_mpf_'): return ctx.make_mpf(libmp.mpf_besseljn(n, z._mpf_, prec, rounding)) elif hasattr(z, '_mpc_'): return ctx.make_mpc(libmp.mpc_besseljn(n, z._mpc_, prec, rounding)) def _agm(ctx, a, b=1): prec, rounding = ctx._prec_rounding if hasattr(a, '_mpf_') and hasattr(b, '_mpf_'): try: v = libmp.mpf_agm(a._mpf_, b._mpf_, prec, rounding) return ctx.make_mpf(v) except ComplexResult: pass if hasattr(a, '_mpf_'): a = (a._mpf_, libmp.fzero) else: a = a._mpc_ if hasattr(b, '_mpf_'): b = (b._mpf_, libmp.fzero) else: b = b._mpc_ return ctx.make_mpc(libmp.mpc_agm(a, b, prec, rounding)) def bernoulli(ctx, n): return ctx.make_mpf(libmp.mpf_bernoulli(int(n), *ctx._prec_rounding)) def _zeta_int(ctx, n): return ctx.make_mpf(libmp.mpf_zeta_int(int(n), *ctx._prec_rounding)) def atan2(ctx, y, x): x = ctx.convert(x) y = ctx.convert(y) return ctx.make_mpf(libmp.mpf_atan2(y._mpf_, x._mpf_, *ctx._prec_rounding)) def psi(ctx, m, z): z = ctx.convert(z) m = int(m) if ctx._is_real_type(z): return ctx.make_mpf(libmp.mpf_psi(m, z._mpf_, *ctx._prec_rounding)) else: return ctx.make_mpc(libmp.mpc_psi(m, z._mpc_, *ctx._prec_rounding)) def cos_sin(ctx, x, **kwargs): if type(x) not in ctx.types: x = ctx.convert(x) prec, rounding = ctx._parse_prec(kwargs) if hasattr(x, '_mpf_'): c, s = libmp.mpf_cos_sin(x._mpf_, prec, rounding) return ctx.make_mpf(c), ctx.make_mpf(s) elif hasattr(x, '_mpc_'): c, s = libmp.mpc_cos_sin(x._mpc_, prec, rounding) return ctx.make_mpc(c), ctx.make_mpc(s) else: return ctx.cos(x, **kwargs), ctx.sin(x, **kwargs) def cospi_sinpi(ctx, x, **kwargs): if type(x) not in ctx.types: x = ctx.convert(x) prec, rounding = ctx._parse_prec(kwargs) if hasattr(x, '_mpf_'): c, s = libmp.mpf_cos_sin_pi(x._mpf_, prec, rounding) return ctx.make_mpf(c), ctx.make_mpf(s) elif hasattr(x, '_mpc_'): c, s = libmp.mpc_cos_sin_pi(x._mpc_, prec, rounding) return ctx.make_mpc(c), ctx.make_mpc(s) else: return ctx.cos(x, **kwargs), ctx.sin(x, **kwargs) def clone(ctx): """ Create a copy of the context, with the same working precision. """ a = ctx.__class__() a.prec = ctx.prec return a # Several helper methods # TODO: add more of these, make consistent, write docstrings, ... def _is_real_type(ctx, x): if hasattr(x, '_mpc_') or type(x) is complex: return False return True def _is_complex_type(ctx, x): if hasattr(x, '_mpc_') or type(x) is complex: return True return False def isnan(ctx, x): """ Return *True* if *x* is a NaN (not-a-number), or for a complex number, whether either the real or complex part is NaN; otherwise return *False*:: >>> from mpmath import * >>> isnan(3.14) False >>> isnan(nan) True >>> isnan(mpc(3.14,2.72)) False >>> isnan(mpc(3.14,nan)) True """ if hasattr(x, "_mpf_"): return x._mpf_ == fnan if hasattr(x, "_mpc_"): return fnan in x._mpc_ if isinstance(x, int_types) or isinstance(x, rational.mpq): return False x = ctx.convert(x) if hasattr(x, '_mpf_') or hasattr(x, '_mpc_'): return ctx.isnan(x) raise TypeError("isnan() needs a number as input") def isfinite(ctx, x): """ Return *True* if *x* is a finite number, i.e. neither an infinity or a NaN. >>> from mpmath import * >>> isfinite(inf) False >>> isfinite(-inf) False >>> isfinite(3) True >>> isfinite(nan) False >>> isfinite(3+4j) True >>> isfinite(mpc(3,inf)) False >>> isfinite(mpc(nan,3)) False """ if ctx.isinf(x) or ctx.isnan(x): return False return True def isnpint(ctx, x): """ Determine if *x* is a nonpositive integer. """ if not x: return True if hasattr(x, '_mpf_'): sign, man, exp, bc = x._mpf_ return sign and exp >= 0 if hasattr(x, '_mpc_'): return not x.imag and ctx.isnpint(x.real) if type(x) in int_types: return x <= 0 if isinstance(x, ctx.mpq): p, q = x._mpq_ if not p: return True return q == 1 and p <= 0 return ctx.isnpint(ctx.convert(x)) def __str__(ctx): lines = ["Mpmath settings:", (" mp.prec = %s" % ctx.prec).ljust(30) + "[default: 53]", (" mp.dps = %s" % ctx.dps).ljust(30) + "[default: 15]", (" mp.trap_complex = %s" % ctx.trap_complex).ljust(30) + "[default: False]", ] return "\n".join(lines) @property def _repr_digits(ctx): return repr_dps(ctx._prec) @property def _str_digits(ctx): return ctx._dps def extraprec(ctx, n, normalize_output=False): """ The block with extraprec(n): increases the precision n bits, executes , and then restores the precision. extraprec(n)(f) returns a decorated version of the function f that increases the working precision by n bits before execution, and restores the parent precision afterwards. With normalize_output=True, it rounds the return value to the parent precision. """ return PrecisionManager(ctx, lambda p: p + n, None, normalize_output) def extradps(ctx, n, normalize_output=False): """ This function is analogous to extraprec (see documentation) but changes the decimal precision instead of the number of bits. """ return PrecisionManager(ctx, None, lambda d: d + n, normalize_output) def workprec(ctx, n, normalize_output=False): """ The block with workprec(n): sets the precision to n bits, executes , and then restores the precision. workprec(n)(f) returns a decorated version of the function f that sets the precision to n bits before execution, and restores the precision afterwards. With normalize_output=True, it rounds the return value to the parent precision. """ return PrecisionManager(ctx, lambda p: n, None, normalize_output) def workdps(ctx, n, normalize_output=False): """ This function is analogous to workprec (see documentation) but changes the decimal precision instead of the number of bits. """ return PrecisionManager(ctx, None, lambda d: n, normalize_output) def autoprec(ctx, f, maxprec=None, catch=(), verbose=False): r""" Return a wrapped copy of *f* that repeatedly evaluates *f* with increasing precision until the result converges to the full precision used at the point of the call. This heuristically protects against rounding errors, at the cost of roughly a 2x slowdown compared to manually setting the optimal precision. This method can, however, easily be fooled if the results from *f* depend "discontinuously" on the precision, for instance if catastrophic cancellation can occur. Therefore, :func:`~mpmath.autoprec` should be used judiciously. **Examples** Many functions are sensitive to perturbations of the input arguments. If the arguments are decimal numbers, they may have to be converted to binary at a much higher precision. If the amount of required extra precision is unknown, :func:`~mpmath.autoprec` is convenient:: >>> from mpmath import * >>> mp.dps = 15 >>> mp.pretty = True >>> besselj(5, 125 * 10**28) # Exact input -8.03284785591801e-17 >>> besselj(5, '1.25e30') # Bad 7.12954868316652e-16 >>> autoprec(besselj)(5, '1.25e30') # Good -8.03284785591801e-17 The following fails to converge because `\sin(\pi) = 0` whereas all finite-precision approximations of `\pi` give nonzero values:: >>> autoprec(sin)(pi) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): ... NoConvergence: autoprec: prec increased to 2910 without convergence As the following example shows, :func:`~mpmath.autoprec` can protect against cancellation, but is fooled by too severe cancellation:: >>> x = 1e-10 >>> exp(x)-1; expm1(x); autoprec(lambda t: exp(t)-1)(x) 1.00000008274037e-10 1.00000000005e-10 1.00000000005e-10 >>> x = 1e-50 >>> exp(x)-1; expm1(x); autoprec(lambda t: exp(t)-1)(x) 0.0 1.0e-50 0.0 With *catch*, an exception or list of exceptions to intercept may be specified. The raised exception is interpreted as signaling insufficient precision. This permits, for example, evaluating a function where a too low precision results in a division by zero:: >>> f = lambda x: 1/(exp(x)-1) >>> f(1e-30) Traceback (most recent call last): ... ZeroDivisionError >>> autoprec(f, catch=ZeroDivisionError)(1e-30) 1.0e+30 """ def f_autoprec_wrapped(*args, **kwargs): prec = ctx.prec if maxprec is None: maxprec2 = ctx._default_hyper_maxprec(prec) else: maxprec2 = maxprec try: ctx.prec = prec + 10 try: v1 = f(*args, **kwargs) except catch: v1 = ctx.nan prec2 = prec + 20 while 1: ctx.prec = prec2 try: v2 = f(*args, **kwargs) except catch: v2 = ctx.nan if v1 == v2: break err = ctx.mag(v2-v1) - ctx.mag(v2) if err < (-prec): break if verbose: print("autoprec: target=%s, prec=%s, accuracy=%s" \ % (prec, prec2, -err)) v1 = v2 if prec2 >= maxprec2: raise ctx.NoConvergence(\ "autoprec: prec increased to %i without convergence"\ % prec2) prec2 += int(prec2*2) prec2 = min(prec2, maxprec2) finally: ctx.prec = prec return +v2 return f_autoprec_wrapped def nstr(ctx, x, n=6, **kwargs): """ Convert an ``mpf`` or ``mpc`` to a decimal string literal with *n* significant digits. The small default value for *n* is chosen to make this function useful for printing collections of numbers (lists, matrices, etc). If *x* is a list or tuple, :func:`~mpmath.nstr` is applied recursively to each element. For unrecognized classes, :func:`~mpmath.nstr` simply returns ``str(x)``. The companion function :func:`~mpmath.nprint` prints the result instead of returning it. The keyword arguments *strip_zeros*, *min_fixed*, *max_fixed* and *show_zero_exponent* are forwarded to :func:`~mpmath.libmp.to_str`. The number will be printed in fixed-point format if the position of the leading digit is strictly between min_fixed (default = min(-dps/3,-5)) and max_fixed (default = dps). To force fixed-point format always, set min_fixed = -inf, max_fixed = +inf. To force floating-point format, set min_fixed >= max_fixed. >>> from mpmath import * >>> nstr([+pi, ldexp(1,-500)]) '[3.14159, 3.05494e-151]' >>> nprint([+pi, ldexp(1,-500)]) [3.14159, 3.05494e-151] >>> nstr(mpf("5e-10"), 5) '5.0e-10' >>> nstr(mpf("5e-10"), 5, strip_zeros=False) '5.0000e-10' >>> nstr(mpf("5e-10"), 5, strip_zeros=False, min_fixed=-11) '0.00000000050000' >>> nstr(mpf(0), 5, show_zero_exponent=True) '0.0e+0' """ if isinstance(x, list): return "[%s]" % (", ".join(ctx.nstr(c, n, **kwargs) for c in x)) if isinstance(x, tuple): return "(%s)" % (", ".join(ctx.nstr(c, n, **kwargs) for c in x)) if hasattr(x, '_mpf_'): return to_str(x._mpf_, n, **kwargs) if hasattr(x, '_mpc_'): return "(" + mpc_to_str(x._mpc_, n, **kwargs) + ")" if isinstance(x, basestring): return repr(x) if isinstance(x, ctx.matrix): return x.__nstr__(n, **kwargs) return str(x) def _convert_fallback(ctx, x, strings): if strings and isinstance(x, basestring): if 'j' in x.lower(): x = x.lower().replace(' ', '') match = get_complex.match(x) re = match.group('re') if not re: re = 0 im = match.group('im').rstrip('j') return ctx.mpc(ctx.convert(re), ctx.convert(im)) if hasattr(x, "_mpi_"): a, b = x._mpi_ if a == b: return ctx.make_mpf(a) else: raise ValueError("can only create mpf from zero-width interval") raise TypeError("cannot create mpf from " + repr(x)) def mpmathify(ctx, *args, **kwargs): return ctx.convert(*args, **kwargs) def _parse_prec(ctx, kwargs): if kwargs: if kwargs.get('exact'): return 0, 'f' prec, rounding = ctx._prec_rounding if 'rounding' in kwargs: rounding = kwargs['rounding'] if 'prec' in kwargs: prec = kwargs['prec'] if prec == ctx.inf: return 0, 'f' else: prec = int(prec) elif 'dps' in kwargs: dps = kwargs['dps'] if dps == ctx.inf: return 0, 'f' prec = dps_to_prec(dps) return prec, rounding return ctx._prec_rounding _exact_overflow_msg = "the exact result does not fit in memory" _hypsum_msg = """hypsum() failed to converge to the requested %i bits of accuracy using a working precision of %i bits. Try with a higher maxprec, maxterms, or set zeroprec.""" def hypsum(ctx, p, q, flags, coeffs, z, accurate_small=True, **kwargs): if hasattr(z, "_mpf_"): key = p, q, flags, 'R' v = z._mpf_ elif hasattr(z, "_mpc_"): key = p, q, flags, 'C' v = z._mpc_ if key not in ctx.hyp_summators: ctx.hyp_summators[key] = libmp.make_hyp_summator(key)[1] summator = ctx.hyp_summators[key] prec = ctx.prec maxprec = kwargs.get('maxprec', ctx._default_hyper_maxprec(prec)) extraprec = 50 epsshift = 25 # Jumps in magnitude occur when parameters are close to negative # integers. We must ensure that these terms are included in # the sum and added accurately magnitude_check = {} max_total_jump = 0 for i, c in enumerate(coeffs): if flags[i] == 'Z': if i >= p and c <= 0: ok = False for ii, cc in enumerate(coeffs[:p]): # Note: c <= cc or c < cc, depending on convention if flags[ii] == 'Z' and cc <= 0 and c <= cc: ok = True if not ok: raise ZeroDivisionError("pole in hypergeometric series") continue n, d = ctx.nint_distance(c) n = -int(n) d = -d if i >= p and n >= 0 and d > 4: if n in magnitude_check: magnitude_check[n] += d else: magnitude_check[n] = d extraprec = max(extraprec, d - prec + 60) max_total_jump += abs(d) while 1: if extraprec > maxprec: raise ValueError(ctx._hypsum_msg % (prec, prec+extraprec)) wp = prec + extraprec if magnitude_check: mag_dict = dict((n,None) for n in magnitude_check) else: mag_dict = {} zv, have_complex, magnitude = summator(coeffs, v, prec, wp, \ epsshift, mag_dict, **kwargs) cancel = -magnitude jumps_resolved = True if extraprec < max_total_jump: for n in mag_dict.values(): if (n is None) or (n < prec): jumps_resolved = False break accurate = (cancel < extraprec-25-5 or not accurate_small) if jumps_resolved: if accurate: break # zero? zeroprec = kwargs.get('zeroprec') if zeroprec is not None: if cancel > zeroprec: if have_complex: return ctx.mpc(0) else: return ctx.zero # Some near-singularities were not included, so increase # precision and repeat until they are extraprec *= 2 # Possible workaround for bad roundoff in fixed-point arithmetic epsshift += 5 extraprec += 5 if type(zv) is tuple: if have_complex: return ctx.make_mpc(zv) else: return ctx.make_mpf(zv) else: return zv def ldexp(ctx, x, n): r""" Computes `x 2^n` efficiently. No rounding is performed. The argument `x` must be a real floating-point number (or possible to convert into one) and `n` must be a Python ``int``. >>> from mpmath import * >>> mp.dps = 15; mp.pretty = False >>> ldexp(1, 10) mpf('1024.0') >>> ldexp(1, -3) mpf('0.125') """ x = ctx.convert(x) return ctx.make_mpf(libmp.mpf_shift(x._mpf_, n)) def frexp(ctx, x): r""" Given a real number `x`, returns `(y, n)` with `y \in [0.5, 1)`, `n` a Python integer, and such that `x = y 2^n`. No rounding is performed. >>> from mpmath import * >>> mp.dps = 15; mp.pretty = False >>> frexp(7.5) (mpf('0.9375'), 3) """ x = ctx.convert(x) y, n = libmp.mpf_frexp(x._mpf_) return ctx.make_mpf(y), n def fneg(ctx, x, **kwargs): """ Negates the number *x*, giving a floating-point result, optionally using a custom precision and rounding mode. See the documentation of :func:`~mpmath.fadd` for a detailed description of how to specify precision and rounding. **Examples** An mpmath number is returned:: >>> from mpmath import * >>> mp.dps = 15; mp.pretty = False >>> fneg(2.5) mpf('-2.5') >>> fneg(-5+2j) mpc(real='5.0', imag='-2.0') Precise control over rounding is possible:: >>> x = fadd(2, 1e-100, exact=True) >>> fneg(x) mpf('-2.0') >>> fneg(x, rounding='f') mpf('-2.0000000000000004') Negating with and without roundoff:: >>> n = 200000000000000000000001 >>> print(int(-mpf(n))) -200000000000000016777216 >>> print(int(fneg(n))) -200000000000000016777216 >>> print(int(fneg(n, prec=log(n,2)+1))) -200000000000000000000001 >>> print(int(fneg(n, dps=log(n,10)+1))) -200000000000000000000001 >>> print(int(fneg(n, prec=inf))) -200000000000000000000001 >>> print(int(fneg(n, dps=inf))) -200000000000000000000001 >>> print(int(fneg(n, exact=True))) -200000000000000000000001 """ prec, rounding = ctx._parse_prec(kwargs) x = ctx.convert(x) if hasattr(x, '_mpf_'): return ctx.make_mpf(mpf_neg(x._mpf_, prec, rounding)) if hasattr(x, '_mpc_'): return ctx.make_mpc(mpc_neg(x._mpc_, prec, rounding)) raise ValueError("Arguments need to be mpf or mpc compatible numbers") def fadd(ctx, x, y, **kwargs): """ Adds the numbers *x* and *y*, giving a floating-point result, optionally using a custom precision and rounding mode. The default precision is the working precision of the context. You can specify a custom precision in bits by passing the *prec* keyword argument, or by providing an equivalent decimal precision with the *dps* keyword argument. If the precision is set to ``+inf``, or if the flag *exact=True* is passed, an exact addition with no rounding is performed. When the precision is finite, the optional *rounding* keyword argument specifies the direction of rounding. Valid options are ``'n'`` for nearest (default), ``'f'`` for floor, ``'c'`` for ceiling, ``'d'`` for down, ``'u'`` for up. **Examples** Using :func:`~mpmath.fadd` with precision and rounding control:: >>> from mpmath import * >>> mp.dps = 15; mp.pretty = False >>> fadd(2, 1e-20) mpf('2.0') >>> fadd(2, 1e-20, rounding='u') mpf('2.0000000000000004') >>> nprint(fadd(2, 1e-20, prec=100), 25) 2.00000000000000000001 >>> nprint(fadd(2, 1e-20, dps=15), 25) 2.0 >>> nprint(fadd(2, 1e-20, dps=25), 25) 2.00000000000000000001 >>> nprint(fadd(2, 1e-20, exact=True), 25) 2.00000000000000000001 Exact addition avoids cancellation errors, enforcing familiar laws of numbers such as `x+y-x = y`, which don't hold in floating-point arithmetic with finite precision:: >>> x, y = mpf(2), mpf('1e-1000') >>> print(x + y - x) 0.0 >>> print(fadd(x, y, prec=inf) - x) 1.0e-1000 >>> print(fadd(x, y, exact=True) - x) 1.0e-1000 Exact addition can be inefficient and may be impossible to perform with large magnitude differences:: >>> fadd(1, '1e-100000000000000000000', prec=inf) Traceback (most recent call last): ... OverflowError: the exact result does not fit in memory """ prec, rounding = ctx._parse_prec(kwargs) x = ctx.convert(x) y = ctx.convert(y) try: if hasattr(x, '_mpf_'): if hasattr(y, '_mpf_'): return ctx.make_mpf(mpf_add(x._mpf_, y._mpf_, prec, rounding)) if hasattr(y, '_mpc_'): return ctx.make_mpc(mpc_add_mpf(y._mpc_, x._mpf_, prec, rounding)) if hasattr(x, '_mpc_'): if hasattr(y, '_mpf_'): return ctx.make_mpc(mpc_add_mpf(x._mpc_, y._mpf_, prec, rounding)) if hasattr(y, '_mpc_'): return ctx.make_mpc(mpc_add(x._mpc_, y._mpc_, prec, rounding)) except (ValueError, OverflowError): raise OverflowError(ctx._exact_overflow_msg) raise ValueError("Arguments need to be mpf or mpc compatible numbers") def fsub(ctx, x, y, **kwargs): """ Subtracts the numbers *x* and *y*, giving a floating-point result, optionally using a custom precision and rounding mode. See the documentation of :func:`~mpmath.fadd` for a detailed description of how to specify precision and rounding. **Examples** Using :func:`~mpmath.fsub` with precision and rounding control:: >>> from mpmath import * >>> mp.dps = 15; mp.pretty = False >>> fsub(2, 1e-20) mpf('2.0') >>> fsub(2, 1e-20, rounding='d') mpf('1.9999999999999998') >>> nprint(fsub(2, 1e-20, prec=100), 25) 1.99999999999999999999 >>> nprint(fsub(2, 1e-20, dps=15), 25) 2.0 >>> nprint(fsub(2, 1e-20, dps=25), 25) 1.99999999999999999999 >>> nprint(fsub(2, 1e-20, exact=True), 25) 1.99999999999999999999 Exact subtraction avoids cancellation errors, enforcing familiar laws of numbers such as `x-y+y = x`, which don't hold in floating-point arithmetic with finite precision:: >>> x, y = mpf(2), mpf('1e1000') >>> print(x - y + y) 0.0 >>> print(fsub(x, y, prec=inf) + y) 2.0 >>> print(fsub(x, y, exact=True) + y) 2.0 Exact addition can be inefficient and may be impossible to perform with large magnitude differences:: >>> fsub(1, '1e-100000000000000000000', prec=inf) Traceback (most recent call last): ... OverflowError: the exact result does not fit in memory """ prec, rounding = ctx._parse_prec(kwargs) x = ctx.convert(x) y = ctx.convert(y) try: if hasattr(x, '_mpf_'): if hasattr(y, '_mpf_'): return ctx.make_mpf(mpf_sub(x._mpf_, y._mpf_, prec, rounding)) if hasattr(y, '_mpc_'): return ctx.make_mpc(mpc_sub((x._mpf_, fzero), y._mpc_, prec, rounding)) if hasattr(x, '_mpc_'): if hasattr(y, '_mpf_'): return ctx.make_mpc(mpc_sub_mpf(x._mpc_, y._mpf_, prec, rounding)) if hasattr(y, '_mpc_'): return ctx.make_mpc(mpc_sub(x._mpc_, y._mpc_, prec, rounding)) except (ValueError, OverflowError): raise OverflowError(ctx._exact_overflow_msg) raise ValueError("Arguments need to be mpf or mpc compatible numbers") def fmul(ctx, x, y, **kwargs): """ Multiplies the numbers *x* and *y*, giving a floating-point result, optionally using a custom precision and rounding mode. See the documentation of :func:`~mpmath.fadd` for a detailed description of how to specify precision and rounding. **Examples** The result is an mpmath number:: >>> from mpmath import * >>> mp.dps = 15; mp.pretty = False >>> fmul(2, 5.0) mpf('10.0') >>> fmul(0.5j, 0.5) mpc(real='0.0', imag='0.25') Avoiding roundoff:: >>> x, y = 10**10+1, 10**15+1 >>> print(x*y) 10000000001000010000000001 >>> print(mpf(x) * mpf(y)) 1.0000000001e+25 >>> print(int(mpf(x) * mpf(y))) 10000000001000011026399232 >>> print(int(fmul(x, y))) 10000000001000011026399232 >>> print(int(fmul(x, y, dps=25))) 10000000001000010000000001 >>> print(int(fmul(x, y, exact=True))) 10000000001000010000000001 Exact multiplication with complex numbers can be inefficient and may be impossible to perform with large magnitude differences between real and imaginary parts:: >>> x = 1+2j >>> y = mpc(2, '1e-100000000000000000000') >>> fmul(x, y) mpc(real='2.0', imag='4.0') >>> fmul(x, y, rounding='u') mpc(real='2.0', imag='4.0000000000000009') >>> fmul(x, y, exact=True) Traceback (most recent call last): ... OverflowError: the exact result does not fit in memory """ prec, rounding = ctx._parse_prec(kwargs) x = ctx.convert(x) y = ctx.convert(y) try: if hasattr(x, '_mpf_'): if hasattr(y, '_mpf_'): return ctx.make_mpf(mpf_mul(x._mpf_, y._mpf_, prec, rounding)) if hasattr(y, '_mpc_'): return ctx.make_mpc(mpc_mul_mpf(y._mpc_, x._mpf_, prec, rounding)) if hasattr(x, '_mpc_'): if hasattr(y, '_mpf_'): return ctx.make_mpc(mpc_mul_mpf(x._mpc_, y._mpf_, prec, rounding)) if hasattr(y, '_mpc_'): return ctx.make_mpc(mpc_mul(x._mpc_, y._mpc_, prec, rounding)) except (ValueError, OverflowError): raise OverflowError(ctx._exact_overflow_msg) raise ValueError("Arguments need to be mpf or mpc compatible numbers") def fdiv(ctx, x, y, **kwargs): """ Divides the numbers *x* and *y*, giving a floating-point result, optionally using a custom precision and rounding mode. See the documentation of :func:`~mpmath.fadd` for a detailed description of how to specify precision and rounding. **Examples** The result is an mpmath number:: >>> from mpmath import * >>> mp.dps = 15; mp.pretty = False >>> fdiv(3, 2) mpf('1.5') >>> fdiv(2, 3) mpf('0.66666666666666663') >>> fdiv(2+4j, 0.5) mpc(real='4.0', imag='8.0') The rounding direction and precision can be controlled:: >>> fdiv(2, 3, dps=3) # Should be accurate to at least 3 digits mpf('0.6666259765625') >>> fdiv(2, 3, rounding='d') mpf('0.66666666666666663') >>> fdiv(2, 3, prec=60) mpf('0.66666666666666667') >>> fdiv(2, 3, rounding='u') mpf('0.66666666666666674') Checking the error of a division by performing it at higher precision:: >>> fdiv(2, 3) - fdiv(2, 3, prec=100) mpf('-3.7007434154172148e-17') Unlike :func:`~mpmath.fadd`, :func:`~mpmath.fmul`, etc., exact division is not allowed since the quotient of two floating-point numbers generally does not have an exact floating-point representation. (In the future this might be changed to allow the case where the division is actually exact.) >>> fdiv(2, 3, exact=True) Traceback (most recent call last): ... ValueError: division is not an exact operation """ prec, rounding = ctx._parse_prec(kwargs) if not prec: raise ValueError("division is not an exact operation") x = ctx.convert(x) y = ctx.convert(y) if hasattr(x, '_mpf_'): if hasattr(y, '_mpf_'): return ctx.make_mpf(mpf_div(x._mpf_, y._mpf_, prec, rounding)) if hasattr(y, '_mpc_'): return ctx.make_mpc(mpc_div((x._mpf_, fzero), y._mpc_, prec, rounding)) if hasattr(x, '_mpc_'): if hasattr(y, '_mpf_'): return ctx.make_mpc(mpc_div_mpf(x._mpc_, y._mpf_, prec, rounding)) if hasattr(y, '_mpc_'): return ctx.make_mpc(mpc_div(x._mpc_, y._mpc_, prec, rounding)) raise ValueError("Arguments need to be mpf or mpc compatible numbers") def nint_distance(ctx, x): r""" Return `(n,d)` where `n` is the nearest integer to `x` and `d` is an estimate of `\log_2(|x-n|)`. If `d < 0`, `-d` gives the precision (measured in bits) lost to cancellation when computing `x-n`. >>> from mpmath import * >>> n, d = nint_distance(5) >>> print(n); print(d) 5 -inf >>> n, d = nint_distance(mpf(5)) >>> print(n); print(d) 5 -inf >>> n, d = nint_distance(mpf(5.00000001)) >>> print(n); print(d) 5 -26 >>> n, d = nint_distance(mpf(4.99999999)) >>> print(n); print(d) 5 -26 >>> n, d = nint_distance(mpc(5,10)) >>> print(n); print(d) 5 4 >>> n, d = nint_distance(mpc(5,0.000001)) >>> print(n); print(d) 5 -19 """ typx = type(x) if typx in int_types: return int(x), ctx.ninf elif typx is rational.mpq: p, q = x._mpq_ n, r = divmod(p, q) if 2*r >= q: n += 1 elif not r: return n, ctx.ninf # log(p/q-n) = log((p-nq)/q) = log(p-nq) - log(q) d = bitcount(abs(p-n*q)) - bitcount(q) return n, d if hasattr(x, "_mpf_"): re = x._mpf_ im_dist = ctx.ninf elif hasattr(x, "_mpc_"): re, im = x._mpc_ isign, iman, iexp, ibc = im if iman: im_dist = iexp + ibc elif im == fzero: im_dist = ctx.ninf else: raise ValueError("requires a finite number") else: x = ctx.convert(x) if hasattr(x, "_mpf_") or hasattr(x, "_mpc_"): return ctx.nint_distance(x) else: raise TypeError("requires an mpf/mpc") sign, man, exp, bc = re mag = exp+bc # |x| < 0.5 if mag < 0: n = 0 re_dist = mag elif man: # exact integer if exp >= 0: n = man << exp re_dist = ctx.ninf # exact half-integer elif exp == -1: n = (man>>1)+1 re_dist = 0 else: d = (-exp-1) t = man >> d if t & 1: t += 1 man = (t<>1 # int(t)>>1 re_dist = exp+bitcount(man) if sign: n = -n elif re == fzero: re_dist = ctx.ninf n = 0 else: raise ValueError("requires a finite number") return n, max(re_dist, im_dist) def fprod(ctx, factors): r""" Calculates a product containing a finite number of factors (for infinite products, see :func:`~mpmath.nprod`). The factors will be converted to mpmath numbers. >>> from mpmath import * >>> mp.dps = 15; mp.pretty = False >>> fprod([1, 2, 0.5, 7]) mpf('7.0') """ orig = ctx.prec try: v = ctx.one for p in factors: v *= p finally: ctx.prec = orig return +v def rand(ctx): """ Returns an ``mpf`` with value chosen randomly from `[0, 1)`. The number of randomly generated bits in the mantissa is equal to the working precision. """ return ctx.make_mpf(mpf_rand(ctx._prec)) def fraction(ctx, p, q): """ Given Python integers `(p, q)`, returns a lazy ``mpf`` representing the fraction `p/q`. The value is updated with the precision. >>> from mpmath import * >>> mp.dps = 15 >>> a = fraction(1,100) >>> b = mpf(1)/100 >>> print(a); print(b) 0.01 0.01 >>> mp.dps = 30 >>> print(a); print(b) # a will be accurate 0.01 0.0100000000000000002081668171172 >>> mp.dps = 15 """ return ctx.constant(lambda prec, rnd: from_rational(p, q, prec, rnd), '%s/%s' % (p, q)) def absmin(ctx, x): return abs(ctx.convert(x)) def absmax(ctx, x): return abs(ctx.convert(x)) def _as_points(ctx, x): # XXX: remove this? if hasattr(x, '_mpi_'): a, b = x._mpi_ return [ctx.make_mpf(a), ctx.make_mpf(b)] return x ''' def _zetasum(ctx, s, a, b): """ Computes sum of k^(-s) for k = a, a+1, ..., b with a, b both small integers. """ a = int(a) b = int(b) s = ctx.convert(s) prec, rounding = ctx._prec_rounding if hasattr(s, '_mpf_'): v = ctx.make_mpf(libmp.mpf_zetasum(s._mpf_, a, b, prec)) elif hasattr(s, '_mpc_'): v = ctx.make_mpc(libmp.mpc_zetasum(s._mpc_, a, b, prec)) return v ''' def _zetasum_fast(ctx, s, a, n, derivatives=[0], reflect=False): if not (ctx.isint(a) and hasattr(s, "_mpc_")): raise NotImplementedError a = int(a) prec = ctx._prec xs, ys = libmp.mpc_zetasum(s._mpc_, a, n, derivatives, reflect, prec) xs = [ctx.make_mpc(x) for x in xs] ys = [ctx.make_mpc(y) for y in ys] return xs, ys class PrecisionManager: def __init__(self, ctx, precfun, dpsfun, normalize_output=False): self.ctx = ctx self.precfun = precfun self.dpsfun = dpsfun self.normalize_output = normalize_output def __call__(self, f): @functools.wraps(f) def g(*args, **kwargs): orig = self.ctx.prec try: if self.precfun: self.ctx.prec = self.precfun(self.ctx.prec) else: self.ctx.dps = self.dpsfun(self.ctx.dps) if self.normalize_output: v = f(*args, **kwargs) if type(v) is tuple: return tuple([+a for a in v]) return +v else: return f(*args, **kwargs) finally: self.ctx.prec = orig return g def __enter__(self): self.origp = self.ctx.prec if self.precfun: self.ctx.prec = self.precfun(self.ctx.prec) else: self.ctx.dps = self.dpsfun(self.ctx.dps) def __exit__(self, exc_type, exc_val, exc_tb): self.ctx.prec = self.origp return False if __name__ == '__main__': import doctest doctest.testmod()