123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253 |
- from .ctx_base import StandardBaseContext
- import math
- import cmath
- from . import math2
- from . import function_docs
- from .libmp import mpf_bernoulli, to_float, int_types
- from . import libmp
- class FPContext(StandardBaseContext):
- """
- Context for fast low-precision arithmetic (53-bit precision, giving at most
- about 15-digit accuracy), using Python's builtin float and complex.
- """
- def __init__(ctx):
- StandardBaseContext.__init__(ctx)
- # Override SpecialFunctions implementation
- ctx.loggamma = math2.loggamma
- ctx._bernoulli_cache = {}
- ctx.pretty = False
- ctx._init_aliases()
- _mpq = lambda cls, x: float(x[0])/x[1]
- NoConvergence = libmp.NoConvergence
- def _get_prec(ctx): return 53
- def _set_prec(ctx, p): return
- def _get_dps(ctx): return 15
- def _set_dps(ctx, p): return
- _fixed_precision = True
- prec = property(_get_prec, _set_prec)
- dps = property(_get_dps, _set_dps)
- zero = 0.0
- one = 1.0
- eps = math2.EPS
- inf = math2.INF
- ninf = math2.NINF
- nan = math2.NAN
- j = 1j
- # Called by SpecialFunctions.__init__()
- @classmethod
- def _wrap_specfun(cls, name, f, wrap):
- if wrap:
- def f_wrapped(ctx, *args, **kwargs):
- convert = ctx.convert
- args = [convert(a) for a in args]
- return f(ctx, *args, **kwargs)
- else:
- f_wrapped = f
- f_wrapped.__doc__ = function_docs.__dict__.get(name, f.__doc__)
- setattr(cls, name, f_wrapped)
- def bernoulli(ctx, n):
- cache = ctx._bernoulli_cache
- if n in cache:
- return cache[n]
- cache[n] = to_float(mpf_bernoulli(n, 53, 'n'), strict=True)
- return cache[n]
- pi = math2.pi
- e = math2.e
- euler = math2.euler
- sqrt2 = 1.4142135623730950488
- sqrt5 = 2.2360679774997896964
- phi = 1.6180339887498948482
- ln2 = 0.69314718055994530942
- ln10 = 2.302585092994045684
- euler = 0.57721566490153286061
- catalan = 0.91596559417721901505
- khinchin = 2.6854520010653064453
- apery = 1.2020569031595942854
- glaisher = 1.2824271291006226369
- absmin = absmax = abs
- def is_special(ctx, x):
- return x - x != 0.0
- def isnan(ctx, x):
- return x != x
- def isinf(ctx, x):
- return abs(x) == math2.INF
- def isnormal(ctx, x):
- if x:
- return x - x == 0.0
- return False
- def isnpint(ctx, x):
- if type(x) is complex:
- if x.imag:
- return False
- x = x.real
- return x <= 0.0 and round(x) == x
- mpf = float
- mpc = complex
- def convert(ctx, x):
- try:
- return float(x)
- except:
- return complex(x)
- power = staticmethod(math2.pow)
- sqrt = staticmethod(math2.sqrt)
- exp = staticmethod(math2.exp)
- ln = log = staticmethod(math2.log)
- cos = staticmethod(math2.cos)
- sin = staticmethod(math2.sin)
- tan = staticmethod(math2.tan)
- cos_sin = staticmethod(math2.cos_sin)
- acos = staticmethod(math2.acos)
- asin = staticmethod(math2.asin)
- atan = staticmethod(math2.atan)
- cosh = staticmethod(math2.cosh)
- sinh = staticmethod(math2.sinh)
- tanh = staticmethod(math2.tanh)
- gamma = staticmethod(math2.gamma)
- rgamma = staticmethod(math2.rgamma)
- fac = factorial = staticmethod(math2.factorial)
- floor = staticmethod(math2.floor)
- ceil = staticmethod(math2.ceil)
- cospi = staticmethod(math2.cospi)
- sinpi = staticmethod(math2.sinpi)
- cbrt = staticmethod(math2.cbrt)
- _nthroot = staticmethod(math2.nthroot)
- _ei = staticmethod(math2.ei)
- _e1 = staticmethod(math2.e1)
- _zeta = _zeta_int = staticmethod(math2.zeta)
- # XXX: math2
- def arg(ctx, z):
- z = complex(z)
- return math.atan2(z.imag, z.real)
- def expj(ctx, x):
- return ctx.exp(ctx.j*x)
- def expjpi(ctx, x):
- return ctx.exp(ctx.j*ctx.pi*x)
- ldexp = math.ldexp
- frexp = math.frexp
- def mag(ctx, z):
- if z:
- return ctx.frexp(abs(z))[1]
- return ctx.ninf
- def isint(ctx, z):
- if hasattr(z, "imag"): # float/int don't have .real/.imag in py2.5
- if z.imag:
- return False
- z = z.real
- try:
- return z == int(z)
- except:
- return False
- def nint_distance(ctx, z):
- if hasattr(z, "imag"): # float/int don't have .real/.imag in py2.5
- n = round(z.real)
- else:
- n = round(z)
- if n == z:
- return n, ctx.ninf
- return n, ctx.mag(abs(z-n))
- def _convert_param(ctx, z):
- if type(z) is tuple:
- p, q = z
- return ctx.mpf(p) / q, 'R'
- if hasattr(z, "imag"): # float/int don't have .real/.imag in py2.5
- intz = int(z.real)
- else:
- intz = int(z)
- if z == intz:
- return intz, 'Z'
- return z, 'R'
- def _is_real_type(ctx, z):
- return isinstance(z, float) or isinstance(z, int_types)
- def _is_complex_type(ctx, z):
- return isinstance(z, complex)
- def hypsum(ctx, p, q, types, coeffs, z, maxterms=6000, **kwargs):
- coeffs = list(coeffs)
- num = range(p)
- den = range(p,p+q)
- tol = ctx.eps
- s = t = 1.0
- k = 0
- while 1:
- for i in num: t *= (coeffs[i]+k)
- for i in den: t /= (coeffs[i]+k)
- k += 1; t /= k; t *= z; s += t
- if abs(t) < tol:
- return s
- if k > maxterms:
- raise ctx.NoConvergence
- def atan2(ctx, x, y):
- return math.atan2(x, y)
- def psi(ctx, m, z):
- m = int(m)
- if m == 0:
- return ctx.digamma(z)
- return (-1)**(m+1) * ctx.fac(m) * ctx.zeta(m+1, z)
- digamma = staticmethod(math2.digamma)
- def harmonic(ctx, x):
- x = ctx.convert(x)
- if x == 0 or x == 1:
- return x
- return ctx.digamma(x+1) + ctx.euler
- nstr = str
- def to_fixed(ctx, x, prec):
- return int(math.ldexp(x, prec))
- def rand(ctx):
- import random
- return random.random()
- _erf = staticmethod(math2.erf)
- _erfc = staticmethod(math2.erfc)
- def sum_accurately(ctx, terms, check_step=1):
- s = ctx.zero
- k = 0
- for term in terms():
- s += term
- if (not k % check_step) and term:
- if abs(term) <= 1e-18*abs(s):
- break
- k += 1
- return s
|