ctx_fp.py 6.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253
  1. from .ctx_base import StandardBaseContext
  2. import math
  3. import cmath
  4. from . import math2
  5. from . import function_docs
  6. from .libmp import mpf_bernoulli, to_float, int_types
  7. from . import libmp
  8. class FPContext(StandardBaseContext):
  9. """
  10. Context for fast low-precision arithmetic (53-bit precision, giving at most
  11. about 15-digit accuracy), using Python's builtin float and complex.
  12. """
  13. def __init__(ctx):
  14. StandardBaseContext.__init__(ctx)
  15. # Override SpecialFunctions implementation
  16. ctx.loggamma = math2.loggamma
  17. ctx._bernoulli_cache = {}
  18. ctx.pretty = False
  19. ctx._init_aliases()
  20. _mpq = lambda cls, x: float(x[0])/x[1]
  21. NoConvergence = libmp.NoConvergence
  22. def _get_prec(ctx): return 53
  23. def _set_prec(ctx, p): return
  24. def _get_dps(ctx): return 15
  25. def _set_dps(ctx, p): return
  26. _fixed_precision = True
  27. prec = property(_get_prec, _set_prec)
  28. dps = property(_get_dps, _set_dps)
  29. zero = 0.0
  30. one = 1.0
  31. eps = math2.EPS
  32. inf = math2.INF
  33. ninf = math2.NINF
  34. nan = math2.NAN
  35. j = 1j
  36. # Called by SpecialFunctions.__init__()
  37. @classmethod
  38. def _wrap_specfun(cls, name, f, wrap):
  39. if wrap:
  40. def f_wrapped(ctx, *args, **kwargs):
  41. convert = ctx.convert
  42. args = [convert(a) for a in args]
  43. return f(ctx, *args, **kwargs)
  44. else:
  45. f_wrapped = f
  46. f_wrapped.__doc__ = function_docs.__dict__.get(name, f.__doc__)
  47. setattr(cls, name, f_wrapped)
  48. def bernoulli(ctx, n):
  49. cache = ctx._bernoulli_cache
  50. if n in cache:
  51. return cache[n]
  52. cache[n] = to_float(mpf_bernoulli(n, 53, 'n'), strict=True)
  53. return cache[n]
  54. pi = math2.pi
  55. e = math2.e
  56. euler = math2.euler
  57. sqrt2 = 1.4142135623730950488
  58. sqrt5 = 2.2360679774997896964
  59. phi = 1.6180339887498948482
  60. ln2 = 0.69314718055994530942
  61. ln10 = 2.302585092994045684
  62. euler = 0.57721566490153286061
  63. catalan = 0.91596559417721901505
  64. khinchin = 2.6854520010653064453
  65. apery = 1.2020569031595942854
  66. glaisher = 1.2824271291006226369
  67. absmin = absmax = abs
  68. def is_special(ctx, x):
  69. return x - x != 0.0
  70. def isnan(ctx, x):
  71. return x != x
  72. def isinf(ctx, x):
  73. return abs(x) == math2.INF
  74. def isnormal(ctx, x):
  75. if x:
  76. return x - x == 0.0
  77. return False
  78. def isnpint(ctx, x):
  79. if type(x) is complex:
  80. if x.imag:
  81. return False
  82. x = x.real
  83. return x <= 0.0 and round(x) == x
  84. mpf = float
  85. mpc = complex
  86. def convert(ctx, x):
  87. try:
  88. return float(x)
  89. except:
  90. return complex(x)
  91. power = staticmethod(math2.pow)
  92. sqrt = staticmethod(math2.sqrt)
  93. exp = staticmethod(math2.exp)
  94. ln = log = staticmethod(math2.log)
  95. cos = staticmethod(math2.cos)
  96. sin = staticmethod(math2.sin)
  97. tan = staticmethod(math2.tan)
  98. cos_sin = staticmethod(math2.cos_sin)
  99. acos = staticmethod(math2.acos)
  100. asin = staticmethod(math2.asin)
  101. atan = staticmethod(math2.atan)
  102. cosh = staticmethod(math2.cosh)
  103. sinh = staticmethod(math2.sinh)
  104. tanh = staticmethod(math2.tanh)
  105. gamma = staticmethod(math2.gamma)
  106. rgamma = staticmethod(math2.rgamma)
  107. fac = factorial = staticmethod(math2.factorial)
  108. floor = staticmethod(math2.floor)
  109. ceil = staticmethod(math2.ceil)
  110. cospi = staticmethod(math2.cospi)
  111. sinpi = staticmethod(math2.sinpi)
  112. cbrt = staticmethod(math2.cbrt)
  113. _nthroot = staticmethod(math2.nthroot)
  114. _ei = staticmethod(math2.ei)
  115. _e1 = staticmethod(math2.e1)
  116. _zeta = _zeta_int = staticmethod(math2.zeta)
  117. # XXX: math2
  118. def arg(ctx, z):
  119. z = complex(z)
  120. return math.atan2(z.imag, z.real)
  121. def expj(ctx, x):
  122. return ctx.exp(ctx.j*x)
  123. def expjpi(ctx, x):
  124. return ctx.exp(ctx.j*ctx.pi*x)
  125. ldexp = math.ldexp
  126. frexp = math.frexp
  127. def mag(ctx, z):
  128. if z:
  129. return ctx.frexp(abs(z))[1]
  130. return ctx.ninf
  131. def isint(ctx, z):
  132. if hasattr(z, "imag"): # float/int don't have .real/.imag in py2.5
  133. if z.imag:
  134. return False
  135. z = z.real
  136. try:
  137. return z == int(z)
  138. except:
  139. return False
  140. def nint_distance(ctx, z):
  141. if hasattr(z, "imag"): # float/int don't have .real/.imag in py2.5
  142. n = round(z.real)
  143. else:
  144. n = round(z)
  145. if n == z:
  146. return n, ctx.ninf
  147. return n, ctx.mag(abs(z-n))
  148. def _convert_param(ctx, z):
  149. if type(z) is tuple:
  150. p, q = z
  151. return ctx.mpf(p) / q, 'R'
  152. if hasattr(z, "imag"): # float/int don't have .real/.imag in py2.5
  153. intz = int(z.real)
  154. else:
  155. intz = int(z)
  156. if z == intz:
  157. return intz, 'Z'
  158. return z, 'R'
  159. def _is_real_type(ctx, z):
  160. return isinstance(z, float) or isinstance(z, int_types)
  161. def _is_complex_type(ctx, z):
  162. return isinstance(z, complex)
  163. def hypsum(ctx, p, q, types, coeffs, z, maxterms=6000, **kwargs):
  164. coeffs = list(coeffs)
  165. num = range(p)
  166. den = range(p,p+q)
  167. tol = ctx.eps
  168. s = t = 1.0
  169. k = 0
  170. while 1:
  171. for i in num: t *= (coeffs[i]+k)
  172. for i in den: t /= (coeffs[i]+k)
  173. k += 1; t /= k; t *= z; s += t
  174. if abs(t) < tol:
  175. return s
  176. if k > maxterms:
  177. raise ctx.NoConvergence
  178. def atan2(ctx, x, y):
  179. return math.atan2(x, y)
  180. def psi(ctx, m, z):
  181. m = int(m)
  182. if m == 0:
  183. return ctx.digamma(z)
  184. return (-1)**(m+1) * ctx.fac(m) * ctx.zeta(m+1, z)
  185. digamma = staticmethod(math2.digamma)
  186. def harmonic(ctx, x):
  187. x = ctx.convert(x)
  188. if x == 0 or x == 1:
  189. return x
  190. return ctx.digamma(x+1) + ctx.euler
  191. nstr = str
  192. def to_fixed(ctx, x, prec):
  193. return int(math.ldexp(x, prec))
  194. def rand(ctx):
  195. import random
  196. return random.random()
  197. _erf = staticmethod(math2.erf)
  198. _erfc = staticmethod(math2.erfc)
  199. def sum_accurately(ctx, terms, check_step=1):
  200. s = ctx.zero
  201. k = 0
  202. for term in terms():
  203. s += term
  204. if (not k % check_step) and term:
  205. if abs(term) <= 1e-18*abs(s):
  206. break
  207. k += 1
  208. return s