Coverage for /usr/lib/python3/dist-packages/mpmath/ctx_mp.py: 27%
628 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2This module defines the mpf, mpc classes, and standard functions for
3operating with them.
4"""
5__docformat__ = 'plaintext'
7import functools
9import re
11from .ctx_base import StandardBaseContext
13from .libmp.backend import basestring, BACKEND
15from . import libmp
17from .libmp import (MPZ, MPZ_ZERO, MPZ_ONE, int_types, repr_dps,
18 round_floor, round_ceiling, dps_to_prec, round_nearest, prec_to_dps,
19 ComplexResult, to_pickable, from_pickable, normalize,
20 from_int, from_float, from_str, to_int, to_float, to_str,
21 from_rational, from_man_exp,
22 fone, fzero, finf, fninf, fnan,
23 mpf_abs, mpf_pos, mpf_neg, mpf_add, mpf_sub, mpf_mul, mpf_mul_int,
24 mpf_div, mpf_rdiv_int, mpf_pow_int, mpf_mod,
25 mpf_eq, mpf_cmp, mpf_lt, mpf_gt, mpf_le, mpf_ge,
26 mpf_hash, mpf_rand,
27 mpf_sum,
28 bitcount, to_fixed,
29 mpc_to_str,
30 mpc_to_complex, mpc_hash, mpc_pos, mpc_is_nonzero, mpc_neg, mpc_conjugate,
31 mpc_abs, mpc_add, mpc_add_mpf, mpc_sub, mpc_sub_mpf, mpc_mul, mpc_mul_mpf,
32 mpc_mul_int, mpc_div, mpc_div_mpf, mpc_pow, mpc_pow_mpf, mpc_pow_int,
33 mpc_mpf_div,
34 mpf_pow,
35 mpf_pi, mpf_degree, mpf_e, mpf_phi, mpf_ln2, mpf_ln10,
36 mpf_euler, mpf_catalan, mpf_apery, mpf_khinchin,
37 mpf_glaisher, mpf_twinprime, mpf_mertens,
38 int_types)
40from . import function_docs
41from . import rational
43new = object.__new__
45get_complex = re.compile(r'^\(?(?P<re>[\+\-]?\d*(\.\d*)?(e[\+\-]?\d+)?)??'
46 r'(?P<im>[\+\-]?\d*(\.\d*)?(e[\+\-]?\d+)?j)?\)?$')
48if BACKEND == 'sage':
49 from sage.libs.mpmath.ext_main import Context as BaseMPContext
50 # pickle hack
51 import sage.libs.mpmath.ext_main as _mpf_module
52else:
53 from .ctx_mp_python import PythonMPContext as BaseMPContext
54 from . import ctx_mp_python as _mpf_module
56from .ctx_mp_python import _mpf, _mpc, mpnumeric
58class MPContext(BaseMPContext, StandardBaseContext):
59 """
60 Context for multiprecision arithmetic with a global precision.
61 """
63 def __init__(ctx):
64 BaseMPContext.__init__(ctx)
65 ctx.trap_complex = False
66 ctx.pretty = False
67 ctx.types = [ctx.mpf, ctx.mpc, ctx.constant]
68 ctx._mpq = rational.mpq
69 ctx.default()
70 StandardBaseContext.__init__(ctx)
72 ctx.mpq = rational.mpq
73 ctx.init_builtins()
75 ctx.hyp_summators = {}
77 ctx._init_aliases()
79 # XXX: automate
80 try:
81 ctx.bernoulli.im_func.func_doc = function_docs.bernoulli
82 ctx.primepi.im_func.func_doc = function_docs.primepi
83 ctx.psi.im_func.func_doc = function_docs.psi
84 ctx.atan2.im_func.func_doc = function_docs.atan2
85 except AttributeError:
86 # python 3
87 ctx.bernoulli.__func__.func_doc = function_docs.bernoulli
88 ctx.primepi.__func__.func_doc = function_docs.primepi
89 ctx.psi.__func__.func_doc = function_docs.psi
90 ctx.atan2.__func__.func_doc = function_docs.atan2
92 ctx.digamma.func_doc = function_docs.digamma
93 ctx.cospi.func_doc = function_docs.cospi
94 ctx.sinpi.func_doc = function_docs.sinpi
96 def init_builtins(ctx):
98 mpf = ctx.mpf
99 mpc = ctx.mpc
101 # Exact constants
102 ctx.one = ctx.make_mpf(fone)
103 ctx.zero = ctx.make_mpf(fzero)
104 ctx.j = ctx.make_mpc((fzero,fone))
105 ctx.inf = ctx.make_mpf(finf)
106 ctx.ninf = ctx.make_mpf(fninf)
107 ctx.nan = ctx.make_mpf(fnan)
109 eps = ctx.constant(lambda prec, rnd: (0, MPZ_ONE, 1-prec, 1),
110 "epsilon of working precision", "eps")
111 ctx.eps = eps
113 # Approximate constants
114 ctx.pi = ctx.constant(mpf_pi, "pi", "pi")
115 ctx.ln2 = ctx.constant(mpf_ln2, "ln(2)", "ln2")
116 ctx.ln10 = ctx.constant(mpf_ln10, "ln(10)", "ln10")
117 ctx.phi = ctx.constant(mpf_phi, "Golden ratio phi", "phi")
118 ctx.e = ctx.constant(mpf_e, "e = exp(1)", "e")
119 ctx.euler = ctx.constant(mpf_euler, "Euler's constant", "euler")
120 ctx.catalan = ctx.constant(mpf_catalan, "Catalan's constant", "catalan")
121 ctx.khinchin = ctx.constant(mpf_khinchin, "Khinchin's constant", "khinchin")
122 ctx.glaisher = ctx.constant(mpf_glaisher, "Glaisher's constant", "glaisher")
123 ctx.apery = ctx.constant(mpf_apery, "Apery's constant", "apery")
124 ctx.degree = ctx.constant(mpf_degree, "1 deg = pi / 180", "degree")
125 ctx.twinprime = ctx.constant(mpf_twinprime, "Twin prime constant", "twinprime")
126 ctx.mertens = ctx.constant(mpf_mertens, "Mertens' constant", "mertens")
128 # Standard functions
129 ctx.sqrt = ctx._wrap_libmp_function(libmp.mpf_sqrt, libmp.mpc_sqrt)
130 ctx.cbrt = ctx._wrap_libmp_function(libmp.mpf_cbrt, libmp.mpc_cbrt)
131 ctx.ln = ctx._wrap_libmp_function(libmp.mpf_log, libmp.mpc_log)
132 ctx.atan = ctx._wrap_libmp_function(libmp.mpf_atan, libmp.mpc_atan)
133 ctx.exp = ctx._wrap_libmp_function(libmp.mpf_exp, libmp.mpc_exp)
134 ctx.expj = ctx._wrap_libmp_function(libmp.mpf_expj, libmp.mpc_expj)
135 ctx.expjpi = ctx._wrap_libmp_function(libmp.mpf_expjpi, libmp.mpc_expjpi)
136 ctx.sin = ctx._wrap_libmp_function(libmp.mpf_sin, libmp.mpc_sin)
137 ctx.cos = ctx._wrap_libmp_function(libmp.mpf_cos, libmp.mpc_cos)
138 ctx.tan = ctx._wrap_libmp_function(libmp.mpf_tan, libmp.mpc_tan)
139 ctx.sinh = ctx._wrap_libmp_function(libmp.mpf_sinh, libmp.mpc_sinh)
140 ctx.cosh = ctx._wrap_libmp_function(libmp.mpf_cosh, libmp.mpc_cosh)
141 ctx.tanh = ctx._wrap_libmp_function(libmp.mpf_tanh, libmp.mpc_tanh)
142 ctx.asin = ctx._wrap_libmp_function(libmp.mpf_asin, libmp.mpc_asin)
143 ctx.acos = ctx._wrap_libmp_function(libmp.mpf_acos, libmp.mpc_acos)
144 ctx.atan = ctx._wrap_libmp_function(libmp.mpf_atan, libmp.mpc_atan)
145 ctx.asinh = ctx._wrap_libmp_function(libmp.mpf_asinh, libmp.mpc_asinh)
146 ctx.acosh = ctx._wrap_libmp_function(libmp.mpf_acosh, libmp.mpc_acosh)
147 ctx.atanh = ctx._wrap_libmp_function(libmp.mpf_atanh, libmp.mpc_atanh)
148 ctx.sinpi = ctx._wrap_libmp_function(libmp.mpf_sin_pi, libmp.mpc_sin_pi)
149 ctx.cospi = ctx._wrap_libmp_function(libmp.mpf_cos_pi, libmp.mpc_cos_pi)
150 ctx.floor = ctx._wrap_libmp_function(libmp.mpf_floor, libmp.mpc_floor)
151 ctx.ceil = ctx._wrap_libmp_function(libmp.mpf_ceil, libmp.mpc_ceil)
152 ctx.nint = ctx._wrap_libmp_function(libmp.mpf_nint, libmp.mpc_nint)
153 ctx.frac = ctx._wrap_libmp_function(libmp.mpf_frac, libmp.mpc_frac)
154 ctx.fib = ctx.fibonacci = ctx._wrap_libmp_function(libmp.mpf_fibonacci, libmp.mpc_fibonacci)
156 ctx.gamma = ctx._wrap_libmp_function(libmp.mpf_gamma, libmp.mpc_gamma)
157 ctx.rgamma = ctx._wrap_libmp_function(libmp.mpf_rgamma, libmp.mpc_rgamma)
158 ctx.loggamma = ctx._wrap_libmp_function(libmp.mpf_loggamma, libmp.mpc_loggamma)
159 ctx.fac = ctx.factorial = ctx._wrap_libmp_function(libmp.mpf_factorial, libmp.mpc_factorial)
161 ctx.digamma = ctx._wrap_libmp_function(libmp.mpf_psi0, libmp.mpc_psi0)
162 ctx.harmonic = ctx._wrap_libmp_function(libmp.mpf_harmonic, libmp.mpc_harmonic)
163 ctx.ei = ctx._wrap_libmp_function(libmp.mpf_ei, libmp.mpc_ei)
164 ctx.e1 = ctx._wrap_libmp_function(libmp.mpf_e1, libmp.mpc_e1)
165 ctx._ci = ctx._wrap_libmp_function(libmp.mpf_ci, libmp.mpc_ci)
166 ctx._si = ctx._wrap_libmp_function(libmp.mpf_si, libmp.mpc_si)
167 ctx.ellipk = ctx._wrap_libmp_function(libmp.mpf_ellipk, libmp.mpc_ellipk)
168 ctx._ellipe = ctx._wrap_libmp_function(libmp.mpf_ellipe, libmp.mpc_ellipe)
169 ctx.agm1 = ctx._wrap_libmp_function(libmp.mpf_agm1, libmp.mpc_agm1)
170 ctx._erf = ctx._wrap_libmp_function(libmp.mpf_erf, None)
171 ctx._erfc = ctx._wrap_libmp_function(libmp.mpf_erfc, None)
172 ctx._zeta = ctx._wrap_libmp_function(libmp.mpf_zeta, libmp.mpc_zeta)
173 ctx._altzeta = ctx._wrap_libmp_function(libmp.mpf_altzeta, libmp.mpc_altzeta)
175 # Faster versions
176 ctx.sqrt = getattr(ctx, "_sage_sqrt", ctx.sqrt)
177 ctx.exp = getattr(ctx, "_sage_exp", ctx.exp)
178 ctx.ln = getattr(ctx, "_sage_ln", ctx.ln)
179 ctx.cos = getattr(ctx, "_sage_cos", ctx.cos)
180 ctx.sin = getattr(ctx, "_sage_sin", ctx.sin)
182 def to_fixed(ctx, x, prec):
183 return x.to_fixed(prec)
185 def hypot(ctx, x, y):
186 r"""
187 Computes the Euclidean norm of the vector `(x, y)`, equal
188 to `\sqrt{x^2 + y^2}`. Both `x` and `y` must be real."""
189 x = ctx.convert(x)
190 y = ctx.convert(y)
191 return ctx.make_mpf(libmp.mpf_hypot(x._mpf_, y._mpf_, *ctx._prec_rounding))
193 def _gamma_upper_int(ctx, n, z):
194 n = int(ctx._re(n))
195 if n == 0:
196 return ctx.e1(z)
197 if not hasattr(z, '_mpf_'):
198 raise NotImplementedError
199 prec, rounding = ctx._prec_rounding
200 real, imag = libmp.mpf_expint(n, z._mpf_, prec, rounding, gamma=True)
201 if imag is None:
202 return ctx.make_mpf(real)
203 else:
204 return ctx.make_mpc((real, imag))
206 def _expint_int(ctx, n, z):
207 n = int(n)
208 if n == 1:
209 return ctx.e1(z)
210 if not hasattr(z, '_mpf_'):
211 raise NotImplementedError
212 prec, rounding = ctx._prec_rounding
213 real, imag = libmp.mpf_expint(n, z._mpf_, prec, rounding)
214 if imag is None:
215 return ctx.make_mpf(real)
216 else:
217 return ctx.make_mpc((real, imag))
219 def _nthroot(ctx, x, n):
220 if hasattr(x, '_mpf_'):
221 try:
222 return ctx.make_mpf(libmp.mpf_nthroot(x._mpf_, n, *ctx._prec_rounding))
223 except ComplexResult:
224 if ctx.trap_complex:
225 raise
226 x = (x._mpf_, libmp.fzero)
227 else:
228 x = x._mpc_
229 return ctx.make_mpc(libmp.mpc_nthroot(x, n, *ctx._prec_rounding))
231 def _besselj(ctx, n, z):
232 prec, rounding = ctx._prec_rounding
233 if hasattr(z, '_mpf_'):
234 return ctx.make_mpf(libmp.mpf_besseljn(n, z._mpf_, prec, rounding))
235 elif hasattr(z, '_mpc_'):
236 return ctx.make_mpc(libmp.mpc_besseljn(n, z._mpc_, prec, rounding))
238 def _agm(ctx, a, b=1):
239 prec, rounding = ctx._prec_rounding
240 if hasattr(a, '_mpf_') and hasattr(b, '_mpf_'):
241 try:
242 v = libmp.mpf_agm(a._mpf_, b._mpf_, prec, rounding)
243 return ctx.make_mpf(v)
244 except ComplexResult:
245 pass
246 if hasattr(a, '_mpf_'): a = (a._mpf_, libmp.fzero)
247 else: a = a._mpc_
248 if hasattr(b, '_mpf_'): b = (b._mpf_, libmp.fzero)
249 else: b = b._mpc_
250 return ctx.make_mpc(libmp.mpc_agm(a, b, prec, rounding))
252 def bernoulli(ctx, n):
253 return ctx.make_mpf(libmp.mpf_bernoulli(int(n), *ctx._prec_rounding))
255 def _zeta_int(ctx, n):
256 return ctx.make_mpf(libmp.mpf_zeta_int(int(n), *ctx._prec_rounding))
258 def atan2(ctx, y, x):
259 x = ctx.convert(x)
260 y = ctx.convert(y)
261 return ctx.make_mpf(libmp.mpf_atan2(y._mpf_, x._mpf_, *ctx._prec_rounding))
263 def psi(ctx, m, z):
264 z = ctx.convert(z)
265 m = int(m)
266 if ctx._is_real_type(z):
267 return ctx.make_mpf(libmp.mpf_psi(m, z._mpf_, *ctx._prec_rounding))
268 else:
269 return ctx.make_mpc(libmp.mpc_psi(m, z._mpc_, *ctx._prec_rounding))
271 def cos_sin(ctx, x, **kwargs):
272 if type(x) not in ctx.types:
273 x = ctx.convert(x)
274 prec, rounding = ctx._parse_prec(kwargs)
275 if hasattr(x, '_mpf_'):
276 c, s = libmp.mpf_cos_sin(x._mpf_, prec, rounding)
277 return ctx.make_mpf(c), ctx.make_mpf(s)
278 elif hasattr(x, '_mpc_'):
279 c, s = libmp.mpc_cos_sin(x._mpc_, prec, rounding)
280 return ctx.make_mpc(c), ctx.make_mpc(s)
281 else:
282 return ctx.cos(x, **kwargs), ctx.sin(x, **kwargs)
284 def cospi_sinpi(ctx, x, **kwargs):
285 if type(x) not in ctx.types:
286 x = ctx.convert(x)
287 prec, rounding = ctx._parse_prec(kwargs)
288 if hasattr(x, '_mpf_'):
289 c, s = libmp.mpf_cos_sin_pi(x._mpf_, prec, rounding)
290 return ctx.make_mpf(c), ctx.make_mpf(s)
291 elif hasattr(x, '_mpc_'):
292 c, s = libmp.mpc_cos_sin_pi(x._mpc_, prec, rounding)
293 return ctx.make_mpc(c), ctx.make_mpc(s)
294 else:
295 return ctx.cos(x, **kwargs), ctx.sin(x, **kwargs)
297 def clone(ctx):
298 """
299 Create a copy of the context, with the same working precision.
300 """
301 a = ctx.__class__()
302 a.prec = ctx.prec
303 return a
305 # Several helper methods
306 # TODO: add more of these, make consistent, write docstrings, ...
308 def _is_real_type(ctx, x):
309 if hasattr(x, '_mpc_') or type(x) is complex:
310 return False
311 return True
313 def _is_complex_type(ctx, x):
314 if hasattr(x, '_mpc_') or type(x) is complex:
315 return True
316 return False
318 def isnan(ctx, x):
319 """
320 Return *True* if *x* is a NaN (not-a-number), or for a complex
321 number, whether either the real or complex part is NaN;
322 otherwise return *False*::
324 >>> from mpmath import *
325 >>> isnan(3.14)
326 False
327 >>> isnan(nan)
328 True
329 >>> isnan(mpc(3.14,2.72))
330 False
331 >>> isnan(mpc(3.14,nan))
332 True
334 """
335 if hasattr(x, "_mpf_"):
336 return x._mpf_ == fnan
337 if hasattr(x, "_mpc_"):
338 return fnan in x._mpc_
339 if isinstance(x, int_types) or isinstance(x, rational.mpq):
340 return False
341 x = ctx.convert(x)
342 if hasattr(x, '_mpf_') or hasattr(x, '_mpc_'):
343 return ctx.isnan(x)
344 raise TypeError("isnan() needs a number as input")
346 def isfinite(ctx, x):
347 """
348 Return *True* if *x* is a finite number, i.e. neither
349 an infinity or a NaN.
351 >>> from mpmath import *
352 >>> isfinite(inf)
353 False
354 >>> isfinite(-inf)
355 False
356 >>> isfinite(3)
357 True
358 >>> isfinite(nan)
359 False
360 >>> isfinite(3+4j)
361 True
362 >>> isfinite(mpc(3,inf))
363 False
364 >>> isfinite(mpc(nan,3))
365 False
367 """
368 if ctx.isinf(x) or ctx.isnan(x):
369 return False
370 return True
372 def isnpint(ctx, x):
373 """
374 Determine if *x* is a nonpositive integer.
375 """
376 if not x:
377 return True
378 if hasattr(x, '_mpf_'):
379 sign, man, exp, bc = x._mpf_
380 return sign and exp >= 0
381 if hasattr(x, '_mpc_'):
382 return not x.imag and ctx.isnpint(x.real)
383 if type(x) in int_types:
384 return x <= 0
385 if isinstance(x, ctx.mpq):
386 p, q = x._mpq_
387 if not p:
388 return True
389 return q == 1 and p <= 0
390 return ctx.isnpint(ctx.convert(x))
392 def __str__(ctx):
393 lines = ["Mpmath settings:",
394 (" mp.prec = %s" % ctx.prec).ljust(30) + "[default: 53]",
395 (" mp.dps = %s" % ctx.dps).ljust(30) + "[default: 15]",
396 (" mp.trap_complex = %s" % ctx.trap_complex).ljust(30) + "[default: False]",
397 ]
398 return "\n".join(lines)
400 @property
401 def _repr_digits(ctx):
402 return repr_dps(ctx._prec)
404 @property
405 def _str_digits(ctx):
406 return ctx._dps
408 def extraprec(ctx, n, normalize_output=False):
409 """
410 The block
412 with extraprec(n):
413 <code>
415 increases the precision n bits, executes <code>, and then
416 restores the precision.
418 extraprec(n)(f) returns a decorated version of the function f
419 that increases the working precision by n bits before execution,
420 and restores the parent precision afterwards. With
421 normalize_output=True, it rounds the return value to the parent
422 precision.
423 """
424 return PrecisionManager(ctx, lambda p: p + n, None, normalize_output)
426 def extradps(ctx, n, normalize_output=False):
427 """
428 This function is analogous to extraprec (see documentation)
429 but changes the decimal precision instead of the number of bits.
430 """
431 return PrecisionManager(ctx, None, lambda d: d + n, normalize_output)
433 def workprec(ctx, n, normalize_output=False):
434 """
435 The block
437 with workprec(n):
438 <code>
440 sets the precision to n bits, executes <code>, and then restores
441 the precision.
443 workprec(n)(f) returns a decorated version of the function f
444 that sets the precision to n bits before execution,
445 and restores the precision afterwards. With normalize_output=True,
446 it rounds the return value to the parent precision.
447 """
448 return PrecisionManager(ctx, lambda p: n, None, normalize_output)
450 def workdps(ctx, n, normalize_output=False):
451 """
452 This function is analogous to workprec (see documentation)
453 but changes the decimal precision instead of the number of bits.
454 """
455 return PrecisionManager(ctx, None, lambda d: n, normalize_output)
457 def autoprec(ctx, f, maxprec=None, catch=(), verbose=False):
458 r"""
459 Return a wrapped copy of *f* that repeatedly evaluates *f*
460 with increasing precision until the result converges to the
461 full precision used at the point of the call.
463 This heuristically protects against rounding errors, at the cost of
464 roughly a 2x slowdown compared to manually setting the optimal
465 precision. This method can, however, easily be fooled if the results
466 from *f* depend "discontinuously" on the precision, for instance
467 if catastrophic cancellation can occur. Therefore, :func:`~mpmath.autoprec`
468 should be used judiciously.
470 **Examples**
472 Many functions are sensitive to perturbations of the input arguments.
473 If the arguments are decimal numbers, they may have to be converted
474 to binary at a much higher precision. If the amount of required
475 extra precision is unknown, :func:`~mpmath.autoprec` is convenient::
477 >>> from mpmath import *
478 >>> mp.dps = 15
479 >>> mp.pretty = True
480 >>> besselj(5, 125 * 10**28) # Exact input
481 -8.03284785591801e-17
482 >>> besselj(5, '1.25e30') # Bad
483 7.12954868316652e-16
484 >>> autoprec(besselj)(5, '1.25e30') # Good
485 -8.03284785591801e-17
487 The following fails to converge because `\sin(\pi) = 0` whereas all
488 finite-precision approximations of `\pi` give nonzero values::
490 >>> autoprec(sin)(pi) # doctest: +IGNORE_EXCEPTION_DETAIL
491 Traceback (most recent call last):
492 ...
493 NoConvergence: autoprec: prec increased to 2910 without convergence
495 As the following example shows, :func:`~mpmath.autoprec` can protect against
496 cancellation, but is fooled by too severe cancellation::
498 >>> x = 1e-10
499 >>> exp(x)-1; expm1(x); autoprec(lambda t: exp(t)-1)(x)
500 1.00000008274037e-10
501 1.00000000005e-10
502 1.00000000005e-10
503 >>> x = 1e-50
504 >>> exp(x)-1; expm1(x); autoprec(lambda t: exp(t)-1)(x)
505 0.0
506 1.0e-50
507 0.0
509 With *catch*, an exception or list of exceptions to intercept
510 may be specified. The raised exception is interpreted
511 as signaling insufficient precision. This permits, for example,
512 evaluating a function where a too low precision results in a
513 division by zero::
515 >>> f = lambda x: 1/(exp(x)-1)
516 >>> f(1e-30)
517 Traceback (most recent call last):
518 ...
519 ZeroDivisionError
520 >>> autoprec(f, catch=ZeroDivisionError)(1e-30)
521 1.0e+30
524 """
525 def f_autoprec_wrapped(*args, **kwargs):
526 prec = ctx.prec
527 if maxprec is None:
528 maxprec2 = ctx._default_hyper_maxprec(prec)
529 else:
530 maxprec2 = maxprec
531 try:
532 ctx.prec = prec + 10
533 try:
534 v1 = f(*args, **kwargs)
535 except catch:
536 v1 = ctx.nan
537 prec2 = prec + 20
538 while 1:
539 ctx.prec = prec2
540 try:
541 v2 = f(*args, **kwargs)
542 except catch:
543 v2 = ctx.nan
544 if v1 == v2:
545 break
546 err = ctx.mag(v2-v1) - ctx.mag(v2)
547 if err < (-prec):
548 break
549 if verbose:
550 print("autoprec: target=%s, prec=%s, accuracy=%s" \
551 % (prec, prec2, -err))
552 v1 = v2
553 if prec2 >= maxprec2:
554 raise ctx.NoConvergence(\
555 "autoprec: prec increased to %i without convergence"\
556 % prec2)
557 prec2 += int(prec2*2)
558 prec2 = min(prec2, maxprec2)
559 finally:
560 ctx.prec = prec
561 return +v2
562 return f_autoprec_wrapped
564 def nstr(ctx, x, n=6, **kwargs):
565 """
566 Convert an ``mpf`` or ``mpc`` to a decimal string literal with *n*
567 significant digits. The small default value for *n* is chosen to
568 make this function useful for printing collections of numbers
569 (lists, matrices, etc).
571 If *x* is a list or tuple, :func:`~mpmath.nstr` is applied recursively
572 to each element. For unrecognized classes, :func:`~mpmath.nstr`
573 simply returns ``str(x)``.
575 The companion function :func:`~mpmath.nprint` prints the result
576 instead of returning it.
578 The keyword arguments *strip_zeros*, *min_fixed*, *max_fixed*
579 and *show_zero_exponent* are forwarded to :func:`~mpmath.libmp.to_str`.
581 The number will be printed in fixed-point format if the position
582 of the leading digit is strictly between min_fixed
583 (default = min(-dps/3,-5)) and max_fixed (default = dps).
585 To force fixed-point format always, set min_fixed = -inf,
586 max_fixed = +inf. To force floating-point format, set
587 min_fixed >= max_fixed.
589 >>> from mpmath import *
590 >>> nstr([+pi, ldexp(1,-500)])
591 '[3.14159, 3.05494e-151]'
592 >>> nprint([+pi, ldexp(1,-500)])
593 [3.14159, 3.05494e-151]
594 >>> nstr(mpf("5e-10"), 5)
595 '5.0e-10'
596 >>> nstr(mpf("5e-10"), 5, strip_zeros=False)
597 '5.0000e-10'
598 >>> nstr(mpf("5e-10"), 5, strip_zeros=False, min_fixed=-11)
599 '0.00000000050000'
600 >>> nstr(mpf(0), 5, show_zero_exponent=True)
601 '0.0e+0'
603 """
604 if isinstance(x, list):
605 return "[%s]" % (", ".join(ctx.nstr(c, n, **kwargs) for c in x))
606 if isinstance(x, tuple):
607 return "(%s)" % (", ".join(ctx.nstr(c, n, **kwargs) for c in x))
608 if hasattr(x, '_mpf_'):
609 return to_str(x._mpf_, n, **kwargs)
610 if hasattr(x, '_mpc_'):
611 return "(" + mpc_to_str(x._mpc_, n, **kwargs) + ")"
612 if isinstance(x, basestring):
613 return repr(x)
614 if isinstance(x, ctx.matrix):
615 return x.__nstr__(n, **kwargs)
616 return str(x)
618 def _convert_fallback(ctx, x, strings):
619 if strings and isinstance(x, basestring):
620 if 'j' in x.lower():
621 x = x.lower().replace(' ', '')
622 match = get_complex.match(x)
623 re = match.group('re')
624 if not re:
625 re = 0
626 im = match.group('im').rstrip('j')
627 return ctx.mpc(ctx.convert(re), ctx.convert(im))
628 if hasattr(x, "_mpi_"):
629 a, b = x._mpi_
630 if a == b:
631 return ctx.make_mpf(a)
632 else:
633 raise ValueError("can only create mpf from zero-width interval")
634 raise TypeError("cannot create mpf from " + repr(x))
636 def mpmathify(ctx, *args, **kwargs):
637 return ctx.convert(*args, **kwargs)
639 def _parse_prec(ctx, kwargs):
640 if kwargs:
641 if kwargs.get('exact'):
642 return 0, 'f'
643 prec, rounding = ctx._prec_rounding
644 if 'rounding' in kwargs:
645 rounding = kwargs['rounding']
646 if 'prec' in kwargs:
647 prec = kwargs['prec']
648 if prec == ctx.inf:
649 return 0, 'f'
650 else:
651 prec = int(prec)
652 elif 'dps' in kwargs:
653 dps = kwargs['dps']
654 if dps == ctx.inf:
655 return 0, 'f'
656 prec = dps_to_prec(dps)
657 return prec, rounding
658 return ctx._prec_rounding
660 _exact_overflow_msg = "the exact result does not fit in memory"
662 _hypsum_msg = """hypsum() failed to converge to the requested %i bits of accuracy
663using a working precision of %i bits. Try with a higher maxprec,
664maxterms, or set zeroprec."""
666 def hypsum(ctx, p, q, flags, coeffs, z, accurate_small=True, **kwargs):
667 if hasattr(z, "_mpf_"):
668 key = p, q, flags, 'R'
669 v = z._mpf_
670 elif hasattr(z, "_mpc_"):
671 key = p, q, flags, 'C'
672 v = z._mpc_
673 if key not in ctx.hyp_summators:
674 ctx.hyp_summators[key] = libmp.make_hyp_summator(key)[1]
675 summator = ctx.hyp_summators[key]
676 prec = ctx.prec
677 maxprec = kwargs.get('maxprec', ctx._default_hyper_maxprec(prec))
678 extraprec = 50
679 epsshift = 25
680 # Jumps in magnitude occur when parameters are close to negative
681 # integers. We must ensure that these terms are included in
682 # the sum and added accurately
683 magnitude_check = {}
684 max_total_jump = 0
685 for i, c in enumerate(coeffs):
686 if flags[i] == 'Z':
687 if i >= p and c <= 0:
688 ok = False
689 for ii, cc in enumerate(coeffs[:p]):
690 # Note: c <= cc or c < cc, depending on convention
691 if flags[ii] == 'Z' and cc <= 0 and c <= cc:
692 ok = True
693 if not ok:
694 raise ZeroDivisionError("pole in hypergeometric series")
695 continue
696 n, d = ctx.nint_distance(c)
697 n = -int(n)
698 d = -d
699 if i >= p and n >= 0 and d > 4:
700 if n in magnitude_check:
701 magnitude_check[n] += d
702 else:
703 magnitude_check[n] = d
704 extraprec = max(extraprec, d - prec + 60)
705 max_total_jump += abs(d)
706 while 1:
707 if extraprec > maxprec:
708 raise ValueError(ctx._hypsum_msg % (prec, prec+extraprec))
709 wp = prec + extraprec
710 if magnitude_check:
711 mag_dict = dict((n,None) for n in magnitude_check)
712 else:
713 mag_dict = {}
714 zv, have_complex, magnitude = summator(coeffs, v, prec, wp, \
715 epsshift, mag_dict, **kwargs)
716 cancel = -magnitude
717 jumps_resolved = True
718 if extraprec < max_total_jump:
719 for n in mag_dict.values():
720 if (n is None) or (n < prec):
721 jumps_resolved = False
722 break
723 accurate = (cancel < extraprec-25-5 or not accurate_small)
724 if jumps_resolved:
725 if accurate:
726 break
727 # zero?
728 zeroprec = kwargs.get('zeroprec')
729 if zeroprec is not None:
730 if cancel > zeroprec:
731 if have_complex:
732 return ctx.mpc(0)
733 else:
734 return ctx.zero
736 # Some near-singularities were not included, so increase
737 # precision and repeat until they are
738 extraprec *= 2
739 # Possible workaround for bad roundoff in fixed-point arithmetic
740 epsshift += 5
741 extraprec += 5
743 if type(zv) is tuple:
744 if have_complex:
745 return ctx.make_mpc(zv)
746 else:
747 return ctx.make_mpf(zv)
748 else:
749 return zv
751 def ldexp(ctx, x, n):
752 r"""
753 Computes `x 2^n` efficiently. No rounding is performed.
754 The argument `x` must be a real floating-point number (or
755 possible to convert into one) and `n` must be a Python ``int``.
757 >>> from mpmath import *
758 >>> mp.dps = 15; mp.pretty = False
759 >>> ldexp(1, 10)
760 mpf('1024.0')
761 >>> ldexp(1, -3)
762 mpf('0.125')
764 """
765 x = ctx.convert(x)
766 return ctx.make_mpf(libmp.mpf_shift(x._mpf_, n))
768 def frexp(ctx, x):
769 r"""
770 Given a real number `x`, returns `(y, n)` with `y \in [0.5, 1)`,
771 `n` a Python integer, and such that `x = y 2^n`. No rounding is
772 performed.
774 >>> from mpmath import *
775 >>> mp.dps = 15; mp.pretty = False
776 >>> frexp(7.5)
777 (mpf('0.9375'), 3)
779 """
780 x = ctx.convert(x)
781 y, n = libmp.mpf_frexp(x._mpf_)
782 return ctx.make_mpf(y), n
784 def fneg(ctx, x, **kwargs):
785 """
786 Negates the number *x*, giving a floating-point result, optionally
787 using a custom precision and rounding mode.
789 See the documentation of :func:`~mpmath.fadd` for a detailed description
790 of how to specify precision and rounding.
792 **Examples**
794 An mpmath number is returned::
796 >>> from mpmath import *
797 >>> mp.dps = 15; mp.pretty = False
798 >>> fneg(2.5)
799 mpf('-2.5')
800 >>> fneg(-5+2j)
801 mpc(real='5.0', imag='-2.0')
803 Precise control over rounding is possible::
805 >>> x = fadd(2, 1e-100, exact=True)
806 >>> fneg(x)
807 mpf('-2.0')
808 >>> fneg(x, rounding='f')
809 mpf('-2.0000000000000004')
811 Negating with and without roundoff::
813 >>> n = 200000000000000000000001
814 >>> print(int(-mpf(n)))
815 -200000000000000016777216
816 >>> print(int(fneg(n)))
817 -200000000000000016777216
818 >>> print(int(fneg(n, prec=log(n,2)+1)))
819 -200000000000000000000001
820 >>> print(int(fneg(n, dps=log(n,10)+1)))
821 -200000000000000000000001
822 >>> print(int(fneg(n, prec=inf)))
823 -200000000000000000000001
824 >>> print(int(fneg(n, dps=inf)))
825 -200000000000000000000001
826 >>> print(int(fneg(n, exact=True)))
827 -200000000000000000000001
829 """
830 prec, rounding = ctx._parse_prec(kwargs)
831 x = ctx.convert(x)
832 if hasattr(x, '_mpf_'):
833 return ctx.make_mpf(mpf_neg(x._mpf_, prec, rounding))
834 if hasattr(x, '_mpc_'):
835 return ctx.make_mpc(mpc_neg(x._mpc_, prec, rounding))
836 raise ValueError("Arguments need to be mpf or mpc compatible numbers")
838 def fadd(ctx, x, y, **kwargs):
839 """
840 Adds the numbers *x* and *y*, giving a floating-point result,
841 optionally using a custom precision and rounding mode.
843 The default precision is the working precision of the context.
844 You can specify a custom precision in bits by passing the *prec* keyword
845 argument, or by providing an equivalent decimal precision with the *dps*
846 keyword argument. If the precision is set to ``+inf``, or if the flag
847 *exact=True* is passed, an exact addition with no rounding is performed.
849 When the precision is finite, the optional *rounding* keyword argument
850 specifies the direction of rounding. Valid options are ``'n'`` for
851 nearest (default), ``'f'`` for floor, ``'c'`` for ceiling, ``'d'``
852 for down, ``'u'`` for up.
854 **Examples**
856 Using :func:`~mpmath.fadd` with precision and rounding control::
858 >>> from mpmath import *
859 >>> mp.dps = 15; mp.pretty = False
860 >>> fadd(2, 1e-20)
861 mpf('2.0')
862 >>> fadd(2, 1e-20, rounding='u')
863 mpf('2.0000000000000004')
864 >>> nprint(fadd(2, 1e-20, prec=100), 25)
865 2.00000000000000000001
866 >>> nprint(fadd(2, 1e-20, dps=15), 25)
867 2.0
868 >>> nprint(fadd(2, 1e-20, dps=25), 25)
869 2.00000000000000000001
870 >>> nprint(fadd(2, 1e-20, exact=True), 25)
871 2.00000000000000000001
873 Exact addition avoids cancellation errors, enforcing familiar laws
874 of numbers such as `x+y-x = y`, which don't hold in floating-point
875 arithmetic with finite precision::
877 >>> x, y = mpf(2), mpf('1e-1000')
878 >>> print(x + y - x)
879 0.0
880 >>> print(fadd(x, y, prec=inf) - x)
881 1.0e-1000
882 >>> print(fadd(x, y, exact=True) - x)
883 1.0e-1000
885 Exact addition can be inefficient and may be impossible to perform
886 with large magnitude differences::
888 >>> fadd(1, '1e-100000000000000000000', prec=inf)
889 Traceback (most recent call last):
890 ...
891 OverflowError: the exact result does not fit in memory
893 """
894 prec, rounding = ctx._parse_prec(kwargs)
895 x = ctx.convert(x)
896 y = ctx.convert(y)
897 try:
898 if hasattr(x, '_mpf_'):
899 if hasattr(y, '_mpf_'):
900 return ctx.make_mpf(mpf_add(x._mpf_, y._mpf_, prec, rounding))
901 if hasattr(y, '_mpc_'):
902 return ctx.make_mpc(mpc_add_mpf(y._mpc_, x._mpf_, prec, rounding))
903 if hasattr(x, '_mpc_'):
904 if hasattr(y, '_mpf_'):
905 return ctx.make_mpc(mpc_add_mpf(x._mpc_, y._mpf_, prec, rounding))
906 if hasattr(y, '_mpc_'):
907 return ctx.make_mpc(mpc_add(x._mpc_, y._mpc_, prec, rounding))
908 except (ValueError, OverflowError):
909 raise OverflowError(ctx._exact_overflow_msg)
910 raise ValueError("Arguments need to be mpf or mpc compatible numbers")
912 def fsub(ctx, x, y, **kwargs):
913 """
914 Subtracts the numbers *x* and *y*, giving a floating-point result,
915 optionally using a custom precision and rounding mode.
917 See the documentation of :func:`~mpmath.fadd` for a detailed description
918 of how to specify precision and rounding.
920 **Examples**
922 Using :func:`~mpmath.fsub` with precision and rounding control::
924 >>> from mpmath import *
925 >>> mp.dps = 15; mp.pretty = False
926 >>> fsub(2, 1e-20)
927 mpf('2.0')
928 >>> fsub(2, 1e-20, rounding='d')
929 mpf('1.9999999999999998')
930 >>> nprint(fsub(2, 1e-20, prec=100), 25)
931 1.99999999999999999999
932 >>> nprint(fsub(2, 1e-20, dps=15), 25)
933 2.0
934 >>> nprint(fsub(2, 1e-20, dps=25), 25)
935 1.99999999999999999999
936 >>> nprint(fsub(2, 1e-20, exact=True), 25)
937 1.99999999999999999999
939 Exact subtraction avoids cancellation errors, enforcing familiar laws
940 of numbers such as `x-y+y = x`, which don't hold in floating-point
941 arithmetic with finite precision::
943 >>> x, y = mpf(2), mpf('1e1000')
944 >>> print(x - y + y)
945 0.0
946 >>> print(fsub(x, y, prec=inf) + y)
947 2.0
948 >>> print(fsub(x, y, exact=True) + y)
949 2.0
951 Exact addition can be inefficient and may be impossible to perform
952 with large magnitude differences::
954 >>> fsub(1, '1e-100000000000000000000', prec=inf)
955 Traceback (most recent call last):
956 ...
957 OverflowError: the exact result does not fit in memory
959 """
960 prec, rounding = ctx._parse_prec(kwargs)
961 x = ctx.convert(x)
962 y = ctx.convert(y)
963 try:
964 if hasattr(x, '_mpf_'):
965 if hasattr(y, '_mpf_'):
966 return ctx.make_mpf(mpf_sub(x._mpf_, y._mpf_, prec, rounding))
967 if hasattr(y, '_mpc_'):
968 return ctx.make_mpc(mpc_sub((x._mpf_, fzero), y._mpc_, prec, rounding))
969 if hasattr(x, '_mpc_'):
970 if hasattr(y, '_mpf_'):
971 return ctx.make_mpc(mpc_sub_mpf(x._mpc_, y._mpf_, prec, rounding))
972 if hasattr(y, '_mpc_'):
973 return ctx.make_mpc(mpc_sub(x._mpc_, y._mpc_, prec, rounding))
974 except (ValueError, OverflowError):
975 raise OverflowError(ctx._exact_overflow_msg)
976 raise ValueError("Arguments need to be mpf or mpc compatible numbers")
978 def fmul(ctx, x, y, **kwargs):
979 """
980 Multiplies the numbers *x* and *y*, giving a floating-point result,
981 optionally using a custom precision and rounding mode.
983 See the documentation of :func:`~mpmath.fadd` for a detailed description
984 of how to specify precision and rounding.
986 **Examples**
988 The result is an mpmath number::
990 >>> from mpmath import *
991 >>> mp.dps = 15; mp.pretty = False
992 >>> fmul(2, 5.0)
993 mpf('10.0')
994 >>> fmul(0.5j, 0.5)
995 mpc(real='0.0', imag='0.25')
997 Avoiding roundoff::
999 >>> x, y = 10**10+1, 10**15+1
1000 >>> print(x*y)
1001 10000000001000010000000001
1002 >>> print(mpf(x) * mpf(y))
1003 1.0000000001e+25
1004 >>> print(int(mpf(x) * mpf(y)))
1005 10000000001000011026399232
1006 >>> print(int(fmul(x, y)))
1007 10000000001000011026399232
1008 >>> print(int(fmul(x, y, dps=25)))
1009 10000000001000010000000001
1010 >>> print(int(fmul(x, y, exact=True)))
1011 10000000001000010000000001
1013 Exact multiplication with complex numbers can be inefficient and may
1014 be impossible to perform with large magnitude differences between
1015 real and imaginary parts::
1017 >>> x = 1+2j
1018 >>> y = mpc(2, '1e-100000000000000000000')
1019 >>> fmul(x, y)
1020 mpc(real='2.0', imag='4.0')
1021 >>> fmul(x, y, rounding='u')
1022 mpc(real='2.0', imag='4.0000000000000009')
1023 >>> fmul(x, y, exact=True)
1024 Traceback (most recent call last):
1025 ...
1026 OverflowError: the exact result does not fit in memory
1028 """
1029 prec, rounding = ctx._parse_prec(kwargs)
1030 x = ctx.convert(x)
1031 y = ctx.convert(y)
1032 try:
1033 if hasattr(x, '_mpf_'):
1034 if hasattr(y, '_mpf_'):
1035 return ctx.make_mpf(mpf_mul(x._mpf_, y._mpf_, prec, rounding))
1036 if hasattr(y, '_mpc_'):
1037 return ctx.make_mpc(mpc_mul_mpf(y._mpc_, x._mpf_, prec, rounding))
1038 if hasattr(x, '_mpc_'):
1039 if hasattr(y, '_mpf_'):
1040 return ctx.make_mpc(mpc_mul_mpf(x._mpc_, y._mpf_, prec, rounding))
1041 if hasattr(y, '_mpc_'):
1042 return ctx.make_mpc(mpc_mul(x._mpc_, y._mpc_, prec, rounding))
1043 except (ValueError, OverflowError):
1044 raise OverflowError(ctx._exact_overflow_msg)
1045 raise ValueError("Arguments need to be mpf or mpc compatible numbers")
1047 def fdiv(ctx, x, y, **kwargs):
1048 """
1049 Divides the numbers *x* and *y*, giving a floating-point result,
1050 optionally using a custom precision and rounding mode.
1052 See the documentation of :func:`~mpmath.fadd` for a detailed description
1053 of how to specify precision and rounding.
1055 **Examples**
1057 The result is an mpmath number::
1059 >>> from mpmath import *
1060 >>> mp.dps = 15; mp.pretty = False
1061 >>> fdiv(3, 2)
1062 mpf('1.5')
1063 >>> fdiv(2, 3)
1064 mpf('0.66666666666666663')
1065 >>> fdiv(2+4j, 0.5)
1066 mpc(real='4.0', imag='8.0')
1068 The rounding direction and precision can be controlled::
1070 >>> fdiv(2, 3, dps=3) # Should be accurate to at least 3 digits
1071 mpf('0.6666259765625')
1072 >>> fdiv(2, 3, rounding='d')
1073 mpf('0.66666666666666663')
1074 >>> fdiv(2, 3, prec=60)
1075 mpf('0.66666666666666667')
1076 >>> fdiv(2, 3, rounding='u')
1077 mpf('0.66666666666666674')
1079 Checking the error of a division by performing it at higher precision::
1081 >>> fdiv(2, 3) - fdiv(2, 3, prec=100)
1082 mpf('-3.7007434154172148e-17')
1084 Unlike :func:`~mpmath.fadd`, :func:`~mpmath.fmul`, etc., exact division is not
1085 allowed since the quotient of two floating-point numbers generally
1086 does not have an exact floating-point representation. (In the
1087 future this might be changed to allow the case where the division
1088 is actually exact.)
1090 >>> fdiv(2, 3, exact=True)
1091 Traceback (most recent call last):
1092 ...
1093 ValueError: division is not an exact operation
1095 """
1096 prec, rounding = ctx._parse_prec(kwargs)
1097 if not prec:
1098 raise ValueError("division is not an exact operation")
1099 x = ctx.convert(x)
1100 y = ctx.convert(y)
1101 if hasattr(x, '_mpf_'):
1102 if hasattr(y, '_mpf_'):
1103 return ctx.make_mpf(mpf_div(x._mpf_, y._mpf_, prec, rounding))
1104 if hasattr(y, '_mpc_'):
1105 return ctx.make_mpc(mpc_div((x._mpf_, fzero), y._mpc_, prec, rounding))
1106 if hasattr(x, '_mpc_'):
1107 if hasattr(y, '_mpf_'):
1108 return ctx.make_mpc(mpc_div_mpf(x._mpc_, y._mpf_, prec, rounding))
1109 if hasattr(y, '_mpc_'):
1110 return ctx.make_mpc(mpc_div(x._mpc_, y._mpc_, prec, rounding))
1111 raise ValueError("Arguments need to be mpf or mpc compatible numbers")
1113 def nint_distance(ctx, x):
1114 r"""
1115 Return `(n,d)` where `n` is the nearest integer to `x` and `d` is
1116 an estimate of `\log_2(|x-n|)`. If `d < 0`, `-d` gives the precision
1117 (measured in bits) lost to cancellation when computing `x-n`.
1119 >>> from mpmath import *
1120 >>> n, d = nint_distance(5)
1121 >>> print(n); print(d)
1122 5
1123 -inf
1124 >>> n, d = nint_distance(mpf(5))
1125 >>> print(n); print(d)
1126 5
1127 -inf
1128 >>> n, d = nint_distance(mpf(5.00000001))
1129 >>> print(n); print(d)
1130 5
1131 -26
1132 >>> n, d = nint_distance(mpf(4.99999999))
1133 >>> print(n); print(d)
1134 5
1135 -26
1136 >>> n, d = nint_distance(mpc(5,10))
1137 >>> print(n); print(d)
1138 5
1139 4
1140 >>> n, d = nint_distance(mpc(5,0.000001))
1141 >>> print(n); print(d)
1142 5
1143 -19
1145 """
1146 typx = type(x)
1147 if typx in int_types:
1148 return int(x), ctx.ninf
1149 elif typx is rational.mpq:
1150 p, q = x._mpq_
1151 n, r = divmod(p, q)
1152 if 2*r >= q:
1153 n += 1
1154 elif not r:
1155 return n, ctx.ninf
1156 # log(p/q-n) = log((p-nq)/q) = log(p-nq) - log(q)
1157 d = bitcount(abs(p-n*q)) - bitcount(q)
1158 return n, d
1159 if hasattr(x, "_mpf_"):
1160 re = x._mpf_
1161 im_dist = ctx.ninf
1162 elif hasattr(x, "_mpc_"):
1163 re, im = x._mpc_
1164 isign, iman, iexp, ibc = im
1165 if iman:
1166 im_dist = iexp + ibc
1167 elif im == fzero:
1168 im_dist = ctx.ninf
1169 else:
1170 raise ValueError("requires a finite number")
1171 else:
1172 x = ctx.convert(x)
1173 if hasattr(x, "_mpf_") or hasattr(x, "_mpc_"):
1174 return ctx.nint_distance(x)
1175 else:
1176 raise TypeError("requires an mpf/mpc")
1177 sign, man, exp, bc = re
1178 mag = exp+bc
1179 # |x| < 0.5
1180 if mag < 0:
1181 n = 0
1182 re_dist = mag
1183 elif man:
1184 # exact integer
1185 if exp >= 0:
1186 n = man << exp
1187 re_dist = ctx.ninf
1188 # exact half-integer
1189 elif exp == -1:
1190 n = (man>>1)+1
1191 re_dist = 0
1192 else:
1193 d = (-exp-1)
1194 t = man >> d
1195 if t & 1:
1196 t += 1
1197 man = (t<<d) - man
1198 else:
1199 man -= (t<<d)
1200 n = t>>1 # int(t)>>1
1201 re_dist = exp+bitcount(man)
1202 if sign:
1203 n = -n
1204 elif re == fzero:
1205 re_dist = ctx.ninf
1206 n = 0
1207 else:
1208 raise ValueError("requires a finite number")
1209 return n, max(re_dist, im_dist)
1211 def fprod(ctx, factors):
1212 r"""
1213 Calculates a product containing a finite number of factors (for
1214 infinite products, see :func:`~mpmath.nprod`). The factors will be
1215 converted to mpmath numbers.
1217 >>> from mpmath import *
1218 >>> mp.dps = 15; mp.pretty = False
1219 >>> fprod([1, 2, 0.5, 7])
1220 mpf('7.0')
1222 """
1223 orig = ctx.prec
1224 try:
1225 v = ctx.one
1226 for p in factors:
1227 v *= p
1228 finally:
1229 ctx.prec = orig
1230 return +v
1232 def rand(ctx):
1233 """
1234 Returns an ``mpf`` with value chosen randomly from `[0, 1)`.
1235 The number of randomly generated bits in the mantissa is equal
1236 to the working precision.
1237 """
1238 return ctx.make_mpf(mpf_rand(ctx._prec))
1240 def fraction(ctx, p, q):
1241 """
1242 Given Python integers `(p, q)`, returns a lazy ``mpf`` representing
1243 the fraction `p/q`. The value is updated with the precision.
1245 >>> from mpmath import *
1246 >>> mp.dps = 15
1247 >>> a = fraction(1,100)
1248 >>> b = mpf(1)/100
1249 >>> print(a); print(b)
1250 0.01
1251 0.01
1252 >>> mp.dps = 30
1253 >>> print(a); print(b) # a will be accurate
1254 0.01
1255 0.0100000000000000002081668171172
1256 >>> mp.dps = 15
1257 """
1258 return ctx.constant(lambda prec, rnd: from_rational(p, q, prec, rnd),
1259 '%s/%s' % (p, q))
1261 def absmin(ctx, x):
1262 return abs(ctx.convert(x))
1264 def absmax(ctx, x):
1265 return abs(ctx.convert(x))
1267 def _as_points(ctx, x):
1268 # XXX: remove this?
1269 if hasattr(x, '_mpi_'):
1270 a, b = x._mpi_
1271 return [ctx.make_mpf(a), ctx.make_mpf(b)]
1272 return x
1274 '''
1275 def _zetasum(ctx, s, a, b):
1276 """
1277 Computes sum of k^(-s) for k = a, a+1, ..., b with a, b both small
1278 integers.
1279 """
1280 a = int(a)
1281 b = int(b)
1282 s = ctx.convert(s)
1283 prec, rounding = ctx._prec_rounding
1284 if hasattr(s, '_mpf_'):
1285 v = ctx.make_mpf(libmp.mpf_zetasum(s._mpf_, a, b, prec))
1286 elif hasattr(s, '_mpc_'):
1287 v = ctx.make_mpc(libmp.mpc_zetasum(s._mpc_, a, b, prec))
1288 return v
1289 '''
1291 def _zetasum_fast(ctx, s, a, n, derivatives=[0], reflect=False):
1292 if not (ctx.isint(a) and hasattr(s, "_mpc_")):
1293 raise NotImplementedError
1294 a = int(a)
1295 prec = ctx._prec
1296 xs, ys = libmp.mpc_zetasum(s._mpc_, a, n, derivatives, reflect, prec)
1297 xs = [ctx.make_mpc(x) for x in xs]
1298 ys = [ctx.make_mpc(y) for y in ys]
1299 return xs, ys
1301class PrecisionManager:
1302 def __init__(self, ctx, precfun, dpsfun, normalize_output=False):
1303 self.ctx = ctx
1304 self.precfun = precfun
1305 self.dpsfun = dpsfun
1306 self.normalize_output = normalize_output
1307 def __call__(self, f):
1308 @functools.wraps(f)
1309 def g(*args, **kwargs):
1310 orig = self.ctx.prec
1311 try:
1312 if self.precfun:
1313 self.ctx.prec = self.precfun(self.ctx.prec)
1314 else:
1315 self.ctx.dps = self.dpsfun(self.ctx.dps)
1316 if self.normalize_output:
1317 v = f(*args, **kwargs)
1318 if type(v) is tuple:
1319 return tuple([+a for a in v])
1320 return +v
1321 else:
1322 return f(*args, **kwargs)
1323 finally:
1324 self.ctx.prec = orig
1325 return g
1326 def __enter__(self):
1327 self.origp = self.ctx.prec
1328 if self.precfun:
1329 self.ctx.prec = self.precfun(self.ctx.prec)
1330 else:
1331 self.ctx.dps = self.dpsfun(self.ctx.dps)
1332 def __exit__(self, exc_type, exc_val, exc_tb):
1333 self.ctx.prec = self.origp
1334 return False
1337if __name__ == '__main__':
1338 import doctest
1339 doctest.testmod()