ctx_mp_python.py 37 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147
  1. #from ctx_base import StandardBaseContext
  2. from .libmp.backend import basestring, exec_
  3. from .libmp import (MPZ, MPZ_ZERO, MPZ_ONE, int_types, repr_dps,
  4. round_floor, round_ceiling, dps_to_prec, round_nearest, prec_to_dps,
  5. ComplexResult, to_pickable, from_pickable, normalize,
  6. from_int, from_float, from_npfloat, from_Decimal, from_str, to_int, to_float, to_str,
  7. from_rational, from_man_exp,
  8. fone, fzero, finf, fninf, fnan,
  9. mpf_abs, mpf_pos, mpf_neg, mpf_add, mpf_sub, mpf_mul, mpf_mul_int,
  10. mpf_div, mpf_rdiv_int, mpf_pow_int, mpf_mod,
  11. mpf_eq, mpf_cmp, mpf_lt, mpf_gt, mpf_le, mpf_ge,
  12. mpf_hash, mpf_rand,
  13. mpf_sum,
  14. bitcount, to_fixed,
  15. mpc_to_str,
  16. mpc_to_complex, mpc_hash, mpc_pos, mpc_is_nonzero, mpc_neg, mpc_conjugate,
  17. mpc_abs, mpc_add, mpc_add_mpf, mpc_sub, mpc_sub_mpf, mpc_mul, mpc_mul_mpf,
  18. mpc_mul_int, mpc_div, mpc_div_mpf, mpc_pow, mpc_pow_mpf, mpc_pow_int,
  19. mpc_mpf_div,
  20. mpf_pow,
  21. mpf_pi, mpf_degree, mpf_e, mpf_phi, mpf_ln2, mpf_ln10,
  22. mpf_euler, mpf_catalan, mpf_apery, mpf_khinchin,
  23. mpf_glaisher, mpf_twinprime, mpf_mertens,
  24. int_types)
  25. from . import rational
  26. from . import function_docs
  27. new = object.__new__
  28. class mpnumeric(object):
  29. """Base class for mpf and mpc."""
  30. __slots__ = []
  31. def __new__(cls, val):
  32. raise NotImplementedError
  33. class _mpf(mpnumeric):
  34. """
  35. An mpf instance holds a real-valued floating-point number. mpf:s
  36. work analogously to Python floats, but support arbitrary-precision
  37. arithmetic.
  38. """
  39. __slots__ = ['_mpf_']
  40. def __new__(cls, val=fzero, **kwargs):
  41. """A new mpf can be created from a Python float, an int, a
  42. or a decimal string representing a number in floating-point
  43. format."""
  44. prec, rounding = cls.context._prec_rounding
  45. if kwargs:
  46. prec = kwargs.get('prec', prec)
  47. if 'dps' in kwargs:
  48. prec = dps_to_prec(kwargs['dps'])
  49. rounding = kwargs.get('rounding', rounding)
  50. if type(val) is cls:
  51. sign, man, exp, bc = val._mpf_
  52. if (not man) and exp:
  53. return val
  54. v = new(cls)
  55. v._mpf_ = normalize(sign, man, exp, bc, prec, rounding)
  56. return v
  57. elif type(val) is tuple:
  58. if len(val) == 2:
  59. v = new(cls)
  60. v._mpf_ = from_man_exp(val[0], val[1], prec, rounding)
  61. return v
  62. if len(val) == 4:
  63. sign, man, exp, bc = val
  64. v = new(cls)
  65. v._mpf_ = normalize(sign, MPZ(man), exp, bc, prec, rounding)
  66. return v
  67. raise ValueError
  68. else:
  69. v = new(cls)
  70. v._mpf_ = mpf_pos(cls.mpf_convert_arg(val, prec, rounding), prec, rounding)
  71. return v
  72. @classmethod
  73. def mpf_convert_arg(cls, x, prec, rounding):
  74. if isinstance(x, int_types): return from_int(x)
  75. if isinstance(x, float): return from_float(x)
  76. if isinstance(x, basestring): return from_str(x, prec, rounding)
  77. if isinstance(x, cls.context.constant): return x.func(prec, rounding)
  78. if hasattr(x, '_mpf_'): return x._mpf_
  79. if hasattr(x, '_mpmath_'):
  80. t = cls.context.convert(x._mpmath_(prec, rounding))
  81. if hasattr(t, '_mpf_'):
  82. return t._mpf_
  83. if hasattr(x, '_mpi_'):
  84. a, b = x._mpi_
  85. if a == b:
  86. return a
  87. raise ValueError("can only create mpf from zero-width interval")
  88. raise TypeError("cannot create mpf from " + repr(x))
  89. @classmethod
  90. def mpf_convert_rhs(cls, x):
  91. if isinstance(x, int_types): return from_int(x)
  92. if isinstance(x, float): return from_float(x)
  93. if isinstance(x, complex_types): return cls.context.mpc(x)
  94. if isinstance(x, rational.mpq):
  95. p, q = x._mpq_
  96. return from_rational(p, q, cls.context.prec)
  97. if hasattr(x, '_mpf_'): return x._mpf_
  98. if hasattr(x, '_mpmath_'):
  99. t = cls.context.convert(x._mpmath_(*cls.context._prec_rounding))
  100. if hasattr(t, '_mpf_'):
  101. return t._mpf_
  102. return t
  103. return NotImplemented
  104. @classmethod
  105. def mpf_convert_lhs(cls, x):
  106. x = cls.mpf_convert_rhs(x)
  107. if type(x) is tuple:
  108. return cls.context.make_mpf(x)
  109. return x
  110. man_exp = property(lambda self: self._mpf_[1:3])
  111. man = property(lambda self: self._mpf_[1])
  112. exp = property(lambda self: self._mpf_[2])
  113. bc = property(lambda self: self._mpf_[3])
  114. real = property(lambda self: self)
  115. imag = property(lambda self: self.context.zero)
  116. conjugate = lambda self: self
  117. def __getstate__(self): return to_pickable(self._mpf_)
  118. def __setstate__(self, val): self._mpf_ = from_pickable(val)
  119. def __repr__(s):
  120. if s.context.pretty:
  121. return str(s)
  122. return "mpf('%s')" % to_str(s._mpf_, s.context._repr_digits)
  123. def __str__(s): return to_str(s._mpf_, s.context._str_digits)
  124. def __hash__(s): return mpf_hash(s._mpf_)
  125. def __int__(s): return int(to_int(s._mpf_))
  126. def __long__(s): return long(to_int(s._mpf_))
  127. def __float__(s): return to_float(s._mpf_, rnd=s.context._prec_rounding[1])
  128. def __complex__(s): return complex(float(s))
  129. def __nonzero__(s): return s._mpf_ != fzero
  130. __bool__ = __nonzero__
  131. def __abs__(s):
  132. cls, new, (prec, rounding) = s._ctxdata
  133. v = new(cls)
  134. v._mpf_ = mpf_abs(s._mpf_, prec, rounding)
  135. return v
  136. def __pos__(s):
  137. cls, new, (prec, rounding) = s._ctxdata
  138. v = new(cls)
  139. v._mpf_ = mpf_pos(s._mpf_, prec, rounding)
  140. return v
  141. def __neg__(s):
  142. cls, new, (prec, rounding) = s._ctxdata
  143. v = new(cls)
  144. v._mpf_ = mpf_neg(s._mpf_, prec, rounding)
  145. return v
  146. def _cmp(s, t, func):
  147. if hasattr(t, '_mpf_'):
  148. t = t._mpf_
  149. else:
  150. t = s.mpf_convert_rhs(t)
  151. if t is NotImplemented:
  152. return t
  153. return func(s._mpf_, t)
  154. def __cmp__(s, t): return s._cmp(t, mpf_cmp)
  155. def __lt__(s, t): return s._cmp(t, mpf_lt)
  156. def __gt__(s, t): return s._cmp(t, mpf_gt)
  157. def __le__(s, t): return s._cmp(t, mpf_le)
  158. def __ge__(s, t): return s._cmp(t, mpf_ge)
  159. def __ne__(s, t):
  160. v = s.__eq__(t)
  161. if v is NotImplemented:
  162. return v
  163. return not v
  164. def __rsub__(s, t):
  165. cls, new, (prec, rounding) = s._ctxdata
  166. if type(t) in int_types:
  167. v = new(cls)
  168. v._mpf_ = mpf_sub(from_int(t), s._mpf_, prec, rounding)
  169. return v
  170. t = s.mpf_convert_lhs(t)
  171. if t is NotImplemented:
  172. return t
  173. return t - s
  174. def __rdiv__(s, t):
  175. cls, new, (prec, rounding) = s._ctxdata
  176. if isinstance(t, int_types):
  177. v = new(cls)
  178. v._mpf_ = mpf_rdiv_int(t, s._mpf_, prec, rounding)
  179. return v
  180. t = s.mpf_convert_lhs(t)
  181. if t is NotImplemented:
  182. return t
  183. return t / s
  184. def __rpow__(s, t):
  185. t = s.mpf_convert_lhs(t)
  186. if t is NotImplemented:
  187. return t
  188. return t ** s
  189. def __rmod__(s, t):
  190. t = s.mpf_convert_lhs(t)
  191. if t is NotImplemented:
  192. return t
  193. return t % s
  194. def sqrt(s):
  195. return s.context.sqrt(s)
  196. def ae(s, t, rel_eps=None, abs_eps=None):
  197. return s.context.almosteq(s, t, rel_eps, abs_eps)
  198. def to_fixed(self, prec):
  199. return to_fixed(self._mpf_, prec)
  200. def __round__(self, *args):
  201. return round(float(self), *args)
  202. mpf_binary_op = """
  203. def %NAME%(self, other):
  204. mpf, new, (prec, rounding) = self._ctxdata
  205. sval = self._mpf_
  206. if hasattr(other, '_mpf_'):
  207. tval = other._mpf_
  208. %WITH_MPF%
  209. ttype = type(other)
  210. if ttype in int_types:
  211. %WITH_INT%
  212. elif ttype is float:
  213. tval = from_float(other)
  214. %WITH_MPF%
  215. elif hasattr(other, '_mpc_'):
  216. tval = other._mpc_
  217. mpc = type(other)
  218. %WITH_MPC%
  219. elif ttype is complex:
  220. tval = from_float(other.real), from_float(other.imag)
  221. mpc = self.context.mpc
  222. %WITH_MPC%
  223. if isinstance(other, mpnumeric):
  224. return NotImplemented
  225. try:
  226. other = mpf.context.convert(other, strings=False)
  227. except TypeError:
  228. return NotImplemented
  229. return self.%NAME%(other)
  230. """
  231. return_mpf = "; obj = new(mpf); obj._mpf_ = val; return obj"
  232. return_mpc = "; obj = new(mpc); obj._mpc_ = val; return obj"
  233. mpf_pow_same = """
  234. try:
  235. val = mpf_pow(sval, tval, prec, rounding) %s
  236. except ComplexResult:
  237. if mpf.context.trap_complex:
  238. raise
  239. mpc = mpf.context.mpc
  240. val = mpc_pow((sval, fzero), (tval, fzero), prec, rounding) %s
  241. """ % (return_mpf, return_mpc)
  242. def binary_op(name, with_mpf='', with_int='', with_mpc=''):
  243. code = mpf_binary_op
  244. code = code.replace("%WITH_INT%", with_int)
  245. code = code.replace("%WITH_MPC%", with_mpc)
  246. code = code.replace("%WITH_MPF%", with_mpf)
  247. code = code.replace("%NAME%", name)
  248. np = {}
  249. exec_(code, globals(), np)
  250. return np[name]
  251. _mpf.__eq__ = binary_op('__eq__',
  252. 'return mpf_eq(sval, tval)',
  253. 'return mpf_eq(sval, from_int(other))',
  254. 'return (tval[1] == fzero) and mpf_eq(tval[0], sval)')
  255. _mpf.__add__ = binary_op('__add__',
  256. 'val = mpf_add(sval, tval, prec, rounding)' + return_mpf,
  257. 'val = mpf_add(sval, from_int(other), prec, rounding)' + return_mpf,
  258. 'val = mpc_add_mpf(tval, sval, prec, rounding)' + return_mpc)
  259. _mpf.__sub__ = binary_op('__sub__',
  260. 'val = mpf_sub(sval, tval, prec, rounding)' + return_mpf,
  261. 'val = mpf_sub(sval, from_int(other), prec, rounding)' + return_mpf,
  262. 'val = mpc_sub((sval, fzero), tval, prec, rounding)' + return_mpc)
  263. _mpf.__mul__ = binary_op('__mul__',
  264. 'val = mpf_mul(sval, tval, prec, rounding)' + return_mpf,
  265. 'val = mpf_mul_int(sval, other, prec, rounding)' + return_mpf,
  266. 'val = mpc_mul_mpf(tval, sval, prec, rounding)' + return_mpc)
  267. _mpf.__div__ = binary_op('__div__',
  268. 'val = mpf_div(sval, tval, prec, rounding)' + return_mpf,
  269. 'val = mpf_div(sval, from_int(other), prec, rounding)' + return_mpf,
  270. 'val = mpc_mpf_div(sval, tval, prec, rounding)' + return_mpc)
  271. _mpf.__mod__ = binary_op('__mod__',
  272. 'val = mpf_mod(sval, tval, prec, rounding)' + return_mpf,
  273. 'val = mpf_mod(sval, from_int(other), prec, rounding)' + return_mpf,
  274. 'raise NotImplementedError("complex modulo")')
  275. _mpf.__pow__ = binary_op('__pow__',
  276. mpf_pow_same,
  277. 'val = mpf_pow_int(sval, other, prec, rounding)' + return_mpf,
  278. 'val = mpc_pow((sval, fzero), tval, prec, rounding)' + return_mpc)
  279. _mpf.__radd__ = _mpf.__add__
  280. _mpf.__rmul__ = _mpf.__mul__
  281. _mpf.__truediv__ = _mpf.__div__
  282. _mpf.__rtruediv__ = _mpf.__rdiv__
  283. class _constant(_mpf):
  284. """Represents a mathematical constant with dynamic precision.
  285. When printed or used in an arithmetic operation, a constant
  286. is converted to a regular mpf at the working precision. A
  287. regular mpf can also be obtained using the operation +x."""
  288. def __new__(cls, func, name, docname=''):
  289. a = object.__new__(cls)
  290. a.name = name
  291. a.func = func
  292. a.__doc__ = getattr(function_docs, docname, '')
  293. return a
  294. def __call__(self, prec=None, dps=None, rounding=None):
  295. prec2, rounding2 = self.context._prec_rounding
  296. if not prec: prec = prec2
  297. if not rounding: rounding = rounding2
  298. if dps: prec = dps_to_prec(dps)
  299. return self.context.make_mpf(self.func(prec, rounding))
  300. @property
  301. def _mpf_(self):
  302. prec, rounding = self.context._prec_rounding
  303. return self.func(prec, rounding)
  304. def __repr__(self):
  305. return "<%s: %s~>" % (self.name, self.context.nstr(self(dps=15)))
  306. class _mpc(mpnumeric):
  307. """
  308. An mpc represents a complex number using a pair of mpf:s (one
  309. for the real part and another for the imaginary part.) The mpc
  310. class behaves fairly similarly to Python's complex type.
  311. """
  312. __slots__ = ['_mpc_']
  313. def __new__(cls, real=0, imag=0):
  314. s = object.__new__(cls)
  315. if isinstance(real, complex_types):
  316. real, imag = real.real, real.imag
  317. elif hasattr(real, '_mpc_'):
  318. s._mpc_ = real._mpc_
  319. return s
  320. real = cls.context.mpf(real)
  321. imag = cls.context.mpf(imag)
  322. s._mpc_ = (real._mpf_, imag._mpf_)
  323. return s
  324. real = property(lambda self: self.context.make_mpf(self._mpc_[0]))
  325. imag = property(lambda self: self.context.make_mpf(self._mpc_[1]))
  326. def __getstate__(self):
  327. return to_pickable(self._mpc_[0]), to_pickable(self._mpc_[1])
  328. def __setstate__(self, val):
  329. self._mpc_ = from_pickable(val[0]), from_pickable(val[1])
  330. def __repr__(s):
  331. if s.context.pretty:
  332. return str(s)
  333. r = repr(s.real)[4:-1]
  334. i = repr(s.imag)[4:-1]
  335. return "%s(real=%s, imag=%s)" % (type(s).__name__, r, i)
  336. def __str__(s):
  337. return "(%s)" % mpc_to_str(s._mpc_, s.context._str_digits)
  338. def __complex__(s):
  339. return mpc_to_complex(s._mpc_, rnd=s.context._prec_rounding[1])
  340. def __pos__(s):
  341. cls, new, (prec, rounding) = s._ctxdata
  342. v = new(cls)
  343. v._mpc_ = mpc_pos(s._mpc_, prec, rounding)
  344. return v
  345. def __abs__(s):
  346. prec, rounding = s.context._prec_rounding
  347. v = new(s.context.mpf)
  348. v._mpf_ = mpc_abs(s._mpc_, prec, rounding)
  349. return v
  350. def __neg__(s):
  351. cls, new, (prec, rounding) = s._ctxdata
  352. v = new(cls)
  353. v._mpc_ = mpc_neg(s._mpc_, prec, rounding)
  354. return v
  355. def conjugate(s):
  356. cls, new, (prec, rounding) = s._ctxdata
  357. v = new(cls)
  358. v._mpc_ = mpc_conjugate(s._mpc_, prec, rounding)
  359. return v
  360. def __nonzero__(s):
  361. return mpc_is_nonzero(s._mpc_)
  362. __bool__ = __nonzero__
  363. def __hash__(s):
  364. return mpc_hash(s._mpc_)
  365. @classmethod
  366. def mpc_convert_lhs(cls, x):
  367. try:
  368. y = cls.context.convert(x)
  369. return y
  370. except TypeError:
  371. return NotImplemented
  372. def __eq__(s, t):
  373. if not hasattr(t, '_mpc_'):
  374. if isinstance(t, str):
  375. return False
  376. t = s.mpc_convert_lhs(t)
  377. if t is NotImplemented:
  378. return t
  379. return s.real == t.real and s.imag == t.imag
  380. def __ne__(s, t):
  381. b = s.__eq__(t)
  382. if b is NotImplemented:
  383. return b
  384. return not b
  385. def _compare(*args):
  386. raise TypeError("no ordering relation is defined for complex numbers")
  387. __gt__ = _compare
  388. __le__ = _compare
  389. __gt__ = _compare
  390. __ge__ = _compare
  391. def __add__(s, t):
  392. cls, new, (prec, rounding) = s._ctxdata
  393. if not hasattr(t, '_mpc_'):
  394. t = s.mpc_convert_lhs(t)
  395. if t is NotImplemented:
  396. return t
  397. if hasattr(t, '_mpf_'):
  398. v = new(cls)
  399. v._mpc_ = mpc_add_mpf(s._mpc_, t._mpf_, prec, rounding)
  400. return v
  401. v = new(cls)
  402. v._mpc_ = mpc_add(s._mpc_, t._mpc_, prec, rounding)
  403. return v
  404. def __sub__(s, t):
  405. cls, new, (prec, rounding) = s._ctxdata
  406. if not hasattr(t, '_mpc_'):
  407. t = s.mpc_convert_lhs(t)
  408. if t is NotImplemented:
  409. return t
  410. if hasattr(t, '_mpf_'):
  411. v = new(cls)
  412. v._mpc_ = mpc_sub_mpf(s._mpc_, t._mpf_, prec, rounding)
  413. return v
  414. v = new(cls)
  415. v._mpc_ = mpc_sub(s._mpc_, t._mpc_, prec, rounding)
  416. return v
  417. def __mul__(s, t):
  418. cls, new, (prec, rounding) = s._ctxdata
  419. if not hasattr(t, '_mpc_'):
  420. if isinstance(t, int_types):
  421. v = new(cls)
  422. v._mpc_ = mpc_mul_int(s._mpc_, t, prec, rounding)
  423. return v
  424. t = s.mpc_convert_lhs(t)
  425. if t is NotImplemented:
  426. return t
  427. if hasattr(t, '_mpf_'):
  428. v = new(cls)
  429. v._mpc_ = mpc_mul_mpf(s._mpc_, t._mpf_, prec, rounding)
  430. return v
  431. t = s.mpc_convert_lhs(t)
  432. v = new(cls)
  433. v._mpc_ = mpc_mul(s._mpc_, t._mpc_, prec, rounding)
  434. return v
  435. def __div__(s, t):
  436. cls, new, (prec, rounding) = s._ctxdata
  437. if not hasattr(t, '_mpc_'):
  438. t = s.mpc_convert_lhs(t)
  439. if t is NotImplemented:
  440. return t
  441. if hasattr(t, '_mpf_'):
  442. v = new(cls)
  443. v._mpc_ = mpc_div_mpf(s._mpc_, t._mpf_, prec, rounding)
  444. return v
  445. v = new(cls)
  446. v._mpc_ = mpc_div(s._mpc_, t._mpc_, prec, rounding)
  447. return v
  448. def __pow__(s, t):
  449. cls, new, (prec, rounding) = s._ctxdata
  450. if isinstance(t, int_types):
  451. v = new(cls)
  452. v._mpc_ = mpc_pow_int(s._mpc_, t, prec, rounding)
  453. return v
  454. t = s.mpc_convert_lhs(t)
  455. if t is NotImplemented:
  456. return t
  457. v = new(cls)
  458. if hasattr(t, '_mpf_'):
  459. v._mpc_ = mpc_pow_mpf(s._mpc_, t._mpf_, prec, rounding)
  460. else:
  461. v._mpc_ = mpc_pow(s._mpc_, t._mpc_, prec, rounding)
  462. return v
  463. __radd__ = __add__
  464. def __rsub__(s, t):
  465. t = s.mpc_convert_lhs(t)
  466. if t is NotImplemented:
  467. return t
  468. return t - s
  469. def __rmul__(s, t):
  470. cls, new, (prec, rounding) = s._ctxdata
  471. if isinstance(t, int_types):
  472. v = new(cls)
  473. v._mpc_ = mpc_mul_int(s._mpc_, t, prec, rounding)
  474. return v
  475. t = s.mpc_convert_lhs(t)
  476. if t is NotImplemented:
  477. return t
  478. return t * s
  479. def __rdiv__(s, t):
  480. t = s.mpc_convert_lhs(t)
  481. if t is NotImplemented:
  482. return t
  483. return t / s
  484. def __rpow__(s, t):
  485. t = s.mpc_convert_lhs(t)
  486. if t is NotImplemented:
  487. return t
  488. return t ** s
  489. __truediv__ = __div__
  490. __rtruediv__ = __rdiv__
  491. def ae(s, t, rel_eps=None, abs_eps=None):
  492. return s.context.almosteq(s, t, rel_eps, abs_eps)
  493. complex_types = (complex, _mpc)
  494. class PythonMPContext(object):
  495. def __init__(ctx):
  496. ctx._prec_rounding = [53, round_nearest]
  497. ctx.mpf = type('mpf', (_mpf,), {})
  498. ctx.mpc = type('mpc', (_mpc,), {})
  499. ctx.mpf._ctxdata = [ctx.mpf, new, ctx._prec_rounding]
  500. ctx.mpc._ctxdata = [ctx.mpc, new, ctx._prec_rounding]
  501. ctx.mpf.context = ctx
  502. ctx.mpc.context = ctx
  503. ctx.constant = type('constant', (_constant,), {})
  504. ctx.constant._ctxdata = [ctx.mpf, new, ctx._prec_rounding]
  505. ctx.constant.context = ctx
  506. def make_mpf(ctx, v):
  507. a = new(ctx.mpf)
  508. a._mpf_ = v
  509. return a
  510. def make_mpc(ctx, v):
  511. a = new(ctx.mpc)
  512. a._mpc_ = v
  513. return a
  514. def default(ctx):
  515. ctx._prec = ctx._prec_rounding[0] = 53
  516. ctx._dps = 15
  517. ctx.trap_complex = False
  518. def _set_prec(ctx, n):
  519. ctx._prec = ctx._prec_rounding[0] = max(1, int(n))
  520. ctx._dps = prec_to_dps(n)
  521. def _set_dps(ctx, n):
  522. ctx._prec = ctx._prec_rounding[0] = dps_to_prec(n)
  523. ctx._dps = max(1, int(n))
  524. prec = property(lambda ctx: ctx._prec, _set_prec)
  525. dps = property(lambda ctx: ctx._dps, _set_dps)
  526. def convert(ctx, x, strings=True):
  527. """
  528. Converts *x* to an ``mpf`` or ``mpc``. If *x* is of type ``mpf``,
  529. ``mpc``, ``int``, ``float``, ``complex``, the conversion
  530. will be performed losslessly.
  531. If *x* is a string, the result will be rounded to the present
  532. working precision. Strings representing fractions or complex
  533. numbers are permitted.
  534. >>> from mpmath import *
  535. >>> mp.dps = 15; mp.pretty = False
  536. >>> mpmathify(3.5)
  537. mpf('3.5')
  538. >>> mpmathify('2.1')
  539. mpf('2.1000000000000001')
  540. >>> mpmathify('3/4')
  541. mpf('0.75')
  542. >>> mpmathify('2+3j')
  543. mpc(real='2.0', imag='3.0')
  544. """
  545. if type(x) in ctx.types: return x
  546. if isinstance(x, int_types): return ctx.make_mpf(from_int(x))
  547. if isinstance(x, float): return ctx.make_mpf(from_float(x))
  548. if isinstance(x, complex):
  549. return ctx.make_mpc((from_float(x.real), from_float(x.imag)))
  550. if type(x).__module__ == 'numpy': return ctx.npconvert(x)
  551. if isinstance(x, numbers.Rational): # e.g. Fraction
  552. try: x = rational.mpq(int(x.numerator), int(x.denominator))
  553. except: pass
  554. prec, rounding = ctx._prec_rounding
  555. if isinstance(x, rational.mpq):
  556. p, q = x._mpq_
  557. return ctx.make_mpf(from_rational(p, q, prec))
  558. if strings and isinstance(x, basestring):
  559. try:
  560. _mpf_ = from_str(x, prec, rounding)
  561. return ctx.make_mpf(_mpf_)
  562. except ValueError:
  563. pass
  564. if hasattr(x, '_mpf_'): return ctx.make_mpf(x._mpf_)
  565. if hasattr(x, '_mpc_'): return ctx.make_mpc(x._mpc_)
  566. if hasattr(x, '_mpmath_'):
  567. return ctx.convert(x._mpmath_(prec, rounding))
  568. if type(x).__module__ == 'decimal':
  569. try: return ctx.make_mpf(from_Decimal(x, prec, rounding))
  570. except: pass
  571. return ctx._convert_fallback(x, strings)
  572. def npconvert(ctx, x):
  573. """
  574. Converts *x* to an ``mpf`` or ``mpc``. *x* should be a numpy
  575. scalar.
  576. """
  577. import numpy as np
  578. if isinstance(x, np.integer): return ctx.make_mpf(from_int(int(x)))
  579. if isinstance(x, np.floating): return ctx.make_mpf(from_npfloat(x))
  580. if isinstance(x, np.complexfloating):
  581. return ctx.make_mpc((from_npfloat(x.real), from_npfloat(x.imag)))
  582. raise TypeError("cannot create mpf from " + repr(x))
  583. def isnan(ctx, x):
  584. """
  585. Return *True* if *x* is a NaN (not-a-number), or for a complex
  586. number, whether either the real or complex part is NaN;
  587. otherwise return *False*::
  588. >>> from mpmath import *
  589. >>> isnan(3.14)
  590. False
  591. >>> isnan(nan)
  592. True
  593. >>> isnan(mpc(3.14,2.72))
  594. False
  595. >>> isnan(mpc(3.14,nan))
  596. True
  597. """
  598. if hasattr(x, "_mpf_"):
  599. return x._mpf_ == fnan
  600. if hasattr(x, "_mpc_"):
  601. return fnan in x._mpc_
  602. if isinstance(x, int_types) or isinstance(x, rational.mpq):
  603. return False
  604. x = ctx.convert(x)
  605. if hasattr(x, '_mpf_') or hasattr(x, '_mpc_'):
  606. return ctx.isnan(x)
  607. raise TypeError("isnan() needs a number as input")
  608. def isinf(ctx, x):
  609. """
  610. Return *True* if the absolute value of *x* is infinite;
  611. otherwise return *False*::
  612. >>> from mpmath import *
  613. >>> isinf(inf)
  614. True
  615. >>> isinf(-inf)
  616. True
  617. >>> isinf(3)
  618. False
  619. >>> isinf(3+4j)
  620. False
  621. >>> isinf(mpc(3,inf))
  622. True
  623. >>> isinf(mpc(inf,3))
  624. True
  625. """
  626. if hasattr(x, "_mpf_"):
  627. return x._mpf_ in (finf, fninf)
  628. if hasattr(x, "_mpc_"):
  629. re, im = x._mpc_
  630. return re in (finf, fninf) or im in (finf, fninf)
  631. if isinstance(x, int_types) or isinstance(x, rational.mpq):
  632. return False
  633. x = ctx.convert(x)
  634. if hasattr(x, '_mpf_') or hasattr(x, '_mpc_'):
  635. return ctx.isinf(x)
  636. raise TypeError("isinf() needs a number as input")
  637. def isnormal(ctx, x):
  638. """
  639. Determine whether *x* is "normal" in the sense of floating-point
  640. representation; that is, return *False* if *x* is zero, an
  641. infinity or NaN; otherwise return *True*. By extension, a
  642. complex number *x* is considered "normal" if its magnitude is
  643. normal::
  644. >>> from mpmath import *
  645. >>> isnormal(3)
  646. True
  647. >>> isnormal(0)
  648. False
  649. >>> isnormal(inf); isnormal(-inf); isnormal(nan)
  650. False
  651. False
  652. False
  653. >>> isnormal(0+0j)
  654. False
  655. >>> isnormal(0+3j)
  656. True
  657. >>> isnormal(mpc(2,nan))
  658. False
  659. """
  660. if hasattr(x, "_mpf_"):
  661. return bool(x._mpf_[1])
  662. if hasattr(x, "_mpc_"):
  663. re, im = x._mpc_
  664. re_normal = bool(re[1])
  665. im_normal = bool(im[1])
  666. if re == fzero: return im_normal
  667. if im == fzero: return re_normal
  668. return re_normal and im_normal
  669. if isinstance(x, int_types) or isinstance(x, rational.mpq):
  670. return bool(x)
  671. x = ctx.convert(x)
  672. if hasattr(x, '_mpf_') or hasattr(x, '_mpc_'):
  673. return ctx.isnormal(x)
  674. raise TypeError("isnormal() needs a number as input")
  675. def isint(ctx, x, gaussian=False):
  676. """
  677. Return *True* if *x* is integer-valued; otherwise return
  678. *False*::
  679. >>> from mpmath import *
  680. >>> isint(3)
  681. True
  682. >>> isint(mpf(3))
  683. True
  684. >>> isint(3.2)
  685. False
  686. >>> isint(inf)
  687. False
  688. Optionally, Gaussian integers can be checked for::
  689. >>> isint(3+0j)
  690. True
  691. >>> isint(3+2j)
  692. False
  693. >>> isint(3+2j, gaussian=True)
  694. True
  695. """
  696. if isinstance(x, int_types):
  697. return True
  698. if hasattr(x, "_mpf_"):
  699. sign, man, exp, bc = xval = x._mpf_
  700. return bool((man and exp >= 0) or xval == fzero)
  701. if hasattr(x, "_mpc_"):
  702. re, im = x._mpc_
  703. rsign, rman, rexp, rbc = re
  704. isign, iman, iexp, ibc = im
  705. re_isint = (rman and rexp >= 0) or re == fzero
  706. if gaussian:
  707. im_isint = (iman and iexp >= 0) or im == fzero
  708. return re_isint and im_isint
  709. return re_isint and im == fzero
  710. if isinstance(x, rational.mpq):
  711. p, q = x._mpq_
  712. return p % q == 0
  713. x = ctx.convert(x)
  714. if hasattr(x, '_mpf_') or hasattr(x, '_mpc_'):
  715. return ctx.isint(x, gaussian)
  716. raise TypeError("isint() needs a number as input")
  717. def fsum(ctx, terms, absolute=False, squared=False):
  718. """
  719. Calculates a sum containing a finite number of terms (for infinite
  720. series, see :func:`~mpmath.nsum`). The terms will be converted to
  721. mpmath numbers. For len(terms) > 2, this function is generally
  722. faster and produces more accurate results than the builtin
  723. Python function :func:`sum`.
  724. >>> from mpmath import *
  725. >>> mp.dps = 15; mp.pretty = False
  726. >>> fsum([1, 2, 0.5, 7])
  727. mpf('10.5')
  728. With squared=True each term is squared, and with absolute=True
  729. the absolute value of each term is used.
  730. """
  731. prec, rnd = ctx._prec_rounding
  732. real = []
  733. imag = []
  734. for term in terms:
  735. reval = imval = 0
  736. if hasattr(term, "_mpf_"):
  737. reval = term._mpf_
  738. elif hasattr(term, "_mpc_"):
  739. reval, imval = term._mpc_
  740. else:
  741. term = ctx.convert(term)
  742. if hasattr(term, "_mpf_"):
  743. reval = term._mpf_
  744. elif hasattr(term, "_mpc_"):
  745. reval, imval = term._mpc_
  746. else:
  747. raise NotImplementedError
  748. if imval:
  749. if squared:
  750. if absolute:
  751. real.append(mpf_mul(reval,reval))
  752. real.append(mpf_mul(imval,imval))
  753. else:
  754. reval, imval = mpc_pow_int((reval,imval),2,prec+10)
  755. real.append(reval)
  756. imag.append(imval)
  757. elif absolute:
  758. real.append(mpc_abs((reval,imval), prec))
  759. else:
  760. real.append(reval)
  761. imag.append(imval)
  762. else:
  763. if squared:
  764. reval = mpf_mul(reval, reval)
  765. elif absolute:
  766. reval = mpf_abs(reval)
  767. real.append(reval)
  768. s = mpf_sum(real, prec, rnd, absolute)
  769. if imag:
  770. s = ctx.make_mpc((s, mpf_sum(imag, prec, rnd)))
  771. else:
  772. s = ctx.make_mpf(s)
  773. return s
  774. def fdot(ctx, A, B=None, conjugate=False):
  775. r"""
  776. Computes the dot product of the iterables `A` and `B`,
  777. .. math ::
  778. \sum_{k=0} A_k B_k.
  779. Alternatively, :func:`~mpmath.fdot` accepts a single iterable of pairs.
  780. In other words, ``fdot(A,B)`` and ``fdot(zip(A,B))`` are equivalent.
  781. The elements are automatically converted to mpmath numbers.
  782. With ``conjugate=True``, the elements in the second vector
  783. will be conjugated:
  784. .. math ::
  785. \sum_{k=0} A_k \overline{B_k}
  786. **Examples**
  787. >>> from mpmath import *
  788. >>> mp.dps = 15; mp.pretty = False
  789. >>> A = [2, 1.5, 3]
  790. >>> B = [1, -1, 2]
  791. >>> fdot(A, B)
  792. mpf('6.5')
  793. >>> list(zip(A, B))
  794. [(2, 1), (1.5, -1), (3, 2)]
  795. >>> fdot(_)
  796. mpf('6.5')
  797. >>> A = [2, 1.5, 3j]
  798. >>> B = [1+j, 3, -1-j]
  799. >>> fdot(A, B)
  800. mpc(real='9.5', imag='-1.0')
  801. >>> fdot(A, B, conjugate=True)
  802. mpc(real='3.5', imag='-5.0')
  803. """
  804. if B is not None:
  805. A = zip(A, B)
  806. prec, rnd = ctx._prec_rounding
  807. real = []
  808. imag = []
  809. hasattr_ = hasattr
  810. types = (ctx.mpf, ctx.mpc)
  811. for a, b in A:
  812. if type(a) not in types: a = ctx.convert(a)
  813. if type(b) not in types: b = ctx.convert(b)
  814. a_real = hasattr_(a, "_mpf_")
  815. b_real = hasattr_(b, "_mpf_")
  816. if a_real and b_real:
  817. real.append(mpf_mul(a._mpf_, b._mpf_))
  818. continue
  819. a_complex = hasattr_(a, "_mpc_")
  820. b_complex = hasattr_(b, "_mpc_")
  821. if a_real and b_complex:
  822. aval = a._mpf_
  823. bre, bim = b._mpc_
  824. if conjugate:
  825. bim = mpf_neg(bim)
  826. real.append(mpf_mul(aval, bre))
  827. imag.append(mpf_mul(aval, bim))
  828. elif b_real and a_complex:
  829. are, aim = a._mpc_
  830. bval = b._mpf_
  831. real.append(mpf_mul(are, bval))
  832. imag.append(mpf_mul(aim, bval))
  833. elif a_complex and b_complex:
  834. #re, im = mpc_mul(a._mpc_, b._mpc_, prec+20)
  835. are, aim = a._mpc_
  836. bre, bim = b._mpc_
  837. if conjugate:
  838. bim = mpf_neg(bim)
  839. real.append(mpf_mul(are, bre))
  840. real.append(mpf_neg(mpf_mul(aim, bim)))
  841. imag.append(mpf_mul(are, bim))
  842. imag.append(mpf_mul(aim, bre))
  843. else:
  844. raise NotImplementedError
  845. s = mpf_sum(real, prec, rnd)
  846. if imag:
  847. s = ctx.make_mpc((s, mpf_sum(imag, prec, rnd)))
  848. else:
  849. s = ctx.make_mpf(s)
  850. return s
  851. def _wrap_libmp_function(ctx, mpf_f, mpc_f=None, mpi_f=None, doc="<no doc>"):
  852. """
  853. Given a low-level mpf_ function, and optionally similar functions
  854. for mpc_ and mpi_, defines the function as a context method.
  855. It is assumed that the return type is the same as that of
  856. the input; the exception is that propagation from mpf to mpc is possible
  857. by raising ComplexResult.
  858. """
  859. def f(x, **kwargs):
  860. if type(x) not in ctx.types:
  861. x = ctx.convert(x)
  862. prec, rounding = ctx._prec_rounding
  863. if kwargs:
  864. prec = kwargs.get('prec', prec)
  865. if 'dps' in kwargs:
  866. prec = dps_to_prec(kwargs['dps'])
  867. rounding = kwargs.get('rounding', rounding)
  868. if hasattr(x, '_mpf_'):
  869. try:
  870. return ctx.make_mpf(mpf_f(x._mpf_, prec, rounding))
  871. except ComplexResult:
  872. # Handle propagation to complex
  873. if ctx.trap_complex:
  874. raise
  875. return ctx.make_mpc(mpc_f((x._mpf_, fzero), prec, rounding))
  876. elif hasattr(x, '_mpc_'):
  877. return ctx.make_mpc(mpc_f(x._mpc_, prec, rounding))
  878. raise NotImplementedError("%s of a %s" % (name, type(x)))
  879. name = mpf_f.__name__[4:]
  880. f.__doc__ = function_docs.__dict__.get(name, "Computes the %s of x" % doc)
  881. return f
  882. # Called by SpecialFunctions.__init__()
  883. @classmethod
  884. def _wrap_specfun(cls, name, f, wrap):
  885. if wrap:
  886. def f_wrapped(ctx, *args, **kwargs):
  887. convert = ctx.convert
  888. args = [convert(a) for a in args]
  889. prec = ctx.prec
  890. try:
  891. ctx.prec += 10
  892. retval = f(ctx, *args, **kwargs)
  893. finally:
  894. ctx.prec = prec
  895. return +retval
  896. else:
  897. f_wrapped = f
  898. f_wrapped.__doc__ = function_docs.__dict__.get(name, f.__doc__)
  899. setattr(cls, name, f_wrapped)
  900. def _convert_param(ctx, x):
  901. if hasattr(x, "_mpc_"):
  902. v, im = x._mpc_
  903. if im != fzero:
  904. return x, 'C'
  905. elif hasattr(x, "_mpf_"):
  906. v = x._mpf_
  907. else:
  908. if type(x) in int_types:
  909. return int(x), 'Z'
  910. p = None
  911. if isinstance(x, tuple):
  912. p, q = x
  913. elif hasattr(x, '_mpq_'):
  914. p, q = x._mpq_
  915. elif isinstance(x, basestring) and '/' in x:
  916. p, q = x.split('/')
  917. p = int(p)
  918. q = int(q)
  919. if p is not None:
  920. if not p % q:
  921. return p // q, 'Z'
  922. return ctx.mpq(p,q), 'Q'
  923. x = ctx.convert(x)
  924. if hasattr(x, "_mpc_"):
  925. v, im = x._mpc_
  926. if im != fzero:
  927. return x, 'C'
  928. elif hasattr(x, "_mpf_"):
  929. v = x._mpf_
  930. else:
  931. return x, 'U'
  932. sign, man, exp, bc = v
  933. if man:
  934. if exp >= -4:
  935. if sign:
  936. man = -man
  937. if exp >= 0:
  938. return int(man) << exp, 'Z'
  939. if exp >= -4:
  940. p, q = int(man), (1<<(-exp))
  941. return ctx.mpq(p,q), 'Q'
  942. x = ctx.make_mpf(v)
  943. return x, 'R'
  944. elif not exp:
  945. return 0, 'Z'
  946. else:
  947. return x, 'U'
  948. def _mpf_mag(ctx, x):
  949. sign, man, exp, bc = x
  950. if man:
  951. return exp+bc
  952. if x == fzero:
  953. return ctx.ninf
  954. if x == finf or x == fninf:
  955. return ctx.inf
  956. return ctx.nan
  957. def mag(ctx, x):
  958. """
  959. Quick logarithmic magnitude estimate of a number. Returns an
  960. integer or infinity `m` such that `|x| <= 2^m`. It is not
  961. guaranteed that `m` is an optimal bound, but it will never
  962. be too large by more than 2 (and probably not more than 1).
  963. **Examples**
  964. >>> from mpmath import *
  965. >>> mp.pretty = True
  966. >>> mag(10), mag(10.0), mag(mpf(10)), int(ceil(log(10,2)))
  967. (4, 4, 4, 4)
  968. >>> mag(10j), mag(10+10j)
  969. (4, 5)
  970. >>> mag(0.01), int(ceil(log(0.01,2)))
  971. (-6, -6)
  972. >>> mag(0), mag(inf), mag(-inf), mag(nan)
  973. (-inf, +inf, +inf, nan)
  974. """
  975. if hasattr(x, "_mpf_"):
  976. return ctx._mpf_mag(x._mpf_)
  977. elif hasattr(x, "_mpc_"):
  978. r, i = x._mpc_
  979. if r == fzero:
  980. return ctx._mpf_mag(i)
  981. if i == fzero:
  982. return ctx._mpf_mag(r)
  983. return 1+max(ctx._mpf_mag(r), ctx._mpf_mag(i))
  984. elif isinstance(x, int_types):
  985. if x:
  986. return bitcount(abs(x))
  987. return ctx.ninf
  988. elif isinstance(x, rational.mpq):
  989. p, q = x._mpq_
  990. if p:
  991. return 1 + bitcount(abs(p)) - bitcount(q)
  992. return ctx.ninf
  993. else:
  994. x = ctx.convert(x)
  995. if hasattr(x, "_mpf_") or hasattr(x, "_mpc_"):
  996. return ctx.mag(x)
  997. else:
  998. raise TypeError("requires an mpf/mpc")
  999. # Register with "numbers" ABC
  1000. # We do not subclass, hence we do not use the @abstractmethod checks. While
  1001. # this is less invasive it may turn out that we do not actually support
  1002. # parts of the expected interfaces. See
  1003. # http://docs.python.org/2/library/numbers.html for list of abstract
  1004. # methods.
  1005. try:
  1006. import numbers
  1007. numbers.Complex.register(_mpc)
  1008. numbers.Real.register(_mpf)
  1009. except ImportError:
  1010. pass