ctx_mp_python.py 37 KB

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