254 lines
6.4 KiB
Python
254 lines
6.4 KiB
Python
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
|