mpelements.py 4.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181
  1. #
  2. # This module is deprecated and should not be used any more. The actual
  3. # implementation of RR and CC now uses mpmath's mpf and mpc types directly.
  4. #
  5. """Real and complex elements. """
  6. from sympy.external.gmpy import MPQ
  7. from sympy.polys.domains.domainelement import DomainElement
  8. from sympy.utilities import public
  9. from mpmath.ctx_mp_python import PythonMPContext, _mpf, _mpc, _constant
  10. from mpmath.libmp import (MPZ_ONE, fzero, fone, finf, fninf, fnan,
  11. round_nearest, mpf_mul, repr_dps, int_types,
  12. from_int, from_float, from_str, to_rational)
  13. @public
  14. class RealElement(_mpf, DomainElement):
  15. """An element of a real domain. """
  16. __slots__ = ('__mpf__',)
  17. def _set_mpf(self, val):
  18. self.__mpf__ = val
  19. _mpf_ = property(lambda self: self.__mpf__, _set_mpf)
  20. def parent(self):
  21. return self.context._parent
  22. @public
  23. class ComplexElement(_mpc, DomainElement):
  24. """An element of a complex domain. """
  25. __slots__ = ('__mpc__',)
  26. def _set_mpc(self, val):
  27. self.__mpc__ = val
  28. _mpc_ = property(lambda self: self.__mpc__, _set_mpc)
  29. def parent(self):
  30. return self.context._parent
  31. new = object.__new__
  32. @public
  33. class MPContext(PythonMPContext):
  34. def __init__(ctx, prec=53, dps=None, tol=None, real=False):
  35. ctx._prec_rounding = [prec, round_nearest]
  36. if dps is None:
  37. ctx._set_prec(prec)
  38. else:
  39. ctx._set_dps(dps)
  40. ctx.mpf = RealElement
  41. ctx.mpc = ComplexElement
  42. ctx.mpf._ctxdata = [ctx.mpf, new, ctx._prec_rounding]
  43. ctx.mpc._ctxdata = [ctx.mpc, new, ctx._prec_rounding]
  44. if real:
  45. ctx.mpf.context = ctx
  46. else:
  47. ctx.mpc.context = ctx
  48. ctx.constant = _constant
  49. ctx.constant._ctxdata = [ctx.mpf, new, ctx._prec_rounding]
  50. ctx.constant.context = ctx
  51. ctx.types = [ctx.mpf, ctx.mpc, ctx.constant]
  52. ctx.trap_complex = True
  53. ctx.pretty = True
  54. if tol is None:
  55. ctx.tol = ctx._make_tol()
  56. elif tol is False:
  57. ctx.tol = fzero
  58. else:
  59. ctx.tol = ctx._convert_tol(tol)
  60. ctx.tolerance = ctx.make_mpf(ctx.tol)
  61. if not ctx.tolerance:
  62. ctx.max_denom = 1000000
  63. else:
  64. ctx.max_denom = int(1/ctx.tolerance)
  65. ctx.zero = ctx.make_mpf(fzero)
  66. ctx.one = ctx.make_mpf(fone)
  67. ctx.j = ctx.make_mpc((fzero, fone))
  68. ctx.inf = ctx.make_mpf(finf)
  69. ctx.ninf = ctx.make_mpf(fninf)
  70. ctx.nan = ctx.make_mpf(fnan)
  71. def _make_tol(ctx):
  72. hundred = (0, 25, 2, 5)
  73. eps = (0, MPZ_ONE, 1-ctx.prec, 1)
  74. return mpf_mul(hundred, eps)
  75. def make_tol(ctx):
  76. return ctx.make_mpf(ctx._make_tol())
  77. def _convert_tol(ctx, tol):
  78. if isinstance(tol, int_types):
  79. return from_int(tol)
  80. if isinstance(tol, float):
  81. return from_float(tol)
  82. if hasattr(tol, "_mpf_"):
  83. return tol._mpf_
  84. prec, rounding = ctx._prec_rounding
  85. if isinstance(tol, str):
  86. return from_str(tol, prec, rounding)
  87. raise ValueError("expected a real number, got %s" % tol)
  88. def _convert_fallback(ctx, x, strings):
  89. raise TypeError("cannot create mpf from " + repr(x))
  90. @property
  91. def _repr_digits(ctx):
  92. return repr_dps(ctx._prec)
  93. @property
  94. def _str_digits(ctx):
  95. return ctx._dps
  96. def to_rational(ctx, s, limit=True):
  97. p, q = to_rational(s._mpf_)
  98. # Needed for GROUND_TYPES=flint if gmpy2 is installed because mpmath's
  99. # to_rational() function returns a gmpy2.mpz instance and if MPQ is
  100. # flint.fmpq then MPQ(p, q) will fail.
  101. p = int(p)
  102. if not limit or q <= ctx.max_denom:
  103. return p, q
  104. p0, q0, p1, q1 = 0, 1, 1, 0
  105. n, d = p, q
  106. while True:
  107. a = n//d
  108. q2 = q0 + a*q1
  109. if q2 > ctx.max_denom:
  110. break
  111. p0, q0, p1, q1 = p1, q1, p0 + a*p1, q2
  112. n, d = d, n - a*d
  113. k = (ctx.max_denom - q0)//q1
  114. number = MPQ(p, q)
  115. bound1 = MPQ(p0 + k*p1, q0 + k*q1)
  116. bound2 = MPQ(p1, q1)
  117. if not bound2 or not bound1:
  118. return p, q
  119. elif abs(bound2 - number) <= abs(bound1 - number):
  120. return bound2.numerator, bound2.denominator
  121. else:
  122. return bound1.numerator, bound1.denominator
  123. def almosteq(ctx, s, t, rel_eps=None, abs_eps=None):
  124. t = ctx.convert(t)
  125. if abs_eps is None and rel_eps is None:
  126. rel_eps = abs_eps = ctx.tolerance or ctx.make_tol()
  127. if abs_eps is None:
  128. abs_eps = ctx.convert(rel_eps)
  129. elif rel_eps is None:
  130. rel_eps = ctx.convert(abs_eps)
  131. diff = abs(s-t)
  132. if diff <= abs_eps:
  133. return True
  134. abss = abs(s)
  135. abst = abs(t)
  136. if abss < abst:
  137. err = diff/abst
  138. else:
  139. err = diff/abss
  140. return err <= rel_eps