Coverage for /usr/lib/python3/dist-packages/sympy/polys/domains/mpelements.py: 52%

1"""Real and complex elements. """ 

2 

3 

4from sympy.polys.domains.domainelement import DomainElement 

5from sympy.utilities import public 

6 

7from mpmath.ctx_mp_python import PythonMPContext, _mpf, _mpc, _constant 

8from mpmath.libmp import (MPZ_ONE, fzero, fone, finf, fninf, fnan, 

9 round_nearest, mpf_mul, repr_dps, int_types, 

10 from_int, from_float, from_str, to_rational) 

11from mpmath.rational import mpq 

12 

13 

14@public 

15class RealElement(_mpf, DomainElement): 

16 """An element of a real domain. """ 

17 

18 __slots__ = ('__mpf__',) 

19 

20 def _set_mpf(self, val): 

21 self.__mpf__ = val 

22 

23 _mpf_ = property(lambda self: self.__mpf__, _set_mpf) 

24 

25 def parent(self): 

26 return self.context._parent 

27 

28@public 

29class ComplexElement(_mpc, DomainElement): 

30 """An element of a complex domain. """ 

31 

32 __slots__ = ('__mpc__',) 

33 

34 def _set_mpc(self, val): 

35 self.__mpc__ = val 

36 

37 _mpc_ = property(lambda self: self.__mpc__, _set_mpc) 

38 

39 def parent(self): 

40 return self.context._parent 

41 

42new = object.__new__ 

43 

44@public 

45class MPContext(PythonMPContext): 

46 

47 def __init__(ctx, prec=53, dps=None, tol=None, real=False): 

48 ctx._prec_rounding = [prec, round_nearest] 

49 

50 if dps is None: 

51 ctx._set_prec(prec) 

52 else: 

53 ctx._set_dps(dps) 

54 

55 ctx.mpf = RealElement 

56 ctx.mpc = ComplexElement 

57 ctx.mpf._ctxdata = [ctx.mpf, new, ctx._prec_rounding] 

58 ctx.mpc._ctxdata = [ctx.mpc, new, ctx._prec_rounding] 

59 

60 if real: 

61 ctx.mpf.context = ctx 

62 else: 

63 ctx.mpc.context = ctx 

64 

65 ctx.constant = _constant 

66 ctx.constant._ctxdata = [ctx.mpf, new, ctx._prec_rounding] 

67 ctx.constant.context = ctx 

68 

69 ctx.types = [ctx.mpf, ctx.mpc, ctx.constant] 

70 ctx.trap_complex = True 

71 ctx.pretty = True 

72 

73 if tol is None: 

74 ctx.tol = ctx._make_tol() 

75 elif tol is False: 

76 ctx.tol = fzero 

77 else: 

78 ctx.tol = ctx._convert_tol(tol) 

79 

80 ctx.tolerance = ctx.make_mpf(ctx.tol) 

81 

82 if not ctx.tolerance: 

83 ctx.max_denom = 1000000 

84 else: 

85 ctx.max_denom = int(1/ctx.tolerance) 

86 

87 ctx.zero = ctx.make_mpf(fzero) 

88 ctx.one = ctx.make_mpf(fone) 

89 ctx.j = ctx.make_mpc((fzero, fone)) 

90 ctx.inf = ctx.make_mpf(finf) 

91 ctx.ninf = ctx.make_mpf(fninf) 

92 ctx.nan = ctx.make_mpf(fnan) 

93 

94 def _make_tol(ctx): 

95 hundred = (0, 25, 2, 5) 

96 eps = (0, MPZ_ONE, 1-ctx.prec, 1) 

97 return mpf_mul(hundred, eps) 

98 

99 def make_tol(ctx): 

100 return ctx.make_mpf(ctx._make_tol()) 

101 

102 def _convert_tol(ctx, tol): 

103 if isinstance(tol, int_types): 

104 return from_int(tol) 

105 if isinstance(tol, float): 

106 return from_float(tol) 

107 if hasattr(tol, "_mpf_"): 

108 return tol._mpf_ 

109 prec, rounding = ctx._prec_rounding 

110 if isinstance(tol, str): 

111 return from_str(tol, prec, rounding) 

112 raise ValueError("expected a real number, got %s" % tol) 

113 

114 def _convert_fallback(ctx, x, strings): 

115 raise TypeError("cannot create mpf from " + repr(x)) 

116 

117 @property 

118 def _repr_digits(ctx): 

119 return repr_dps(ctx._prec) 

120 

121 @property 

122 def _str_digits(ctx): 

123 return ctx._dps 

124 

125 def to_rational(ctx, s, limit=True): 

126 p, q = to_rational(s._mpf_) 

127 

128 if not limit or q <= ctx.max_denom: 

129 return p, q 

130 

131 p0, q0, p1, q1 = 0, 1, 1, 0 

132 n, d = p, q 

133 

134 while True: 

135 a = n//d 

136 q2 = q0 + a*q1 

137 if q2 > ctx.max_denom: 

138 break 

139 p0, q0, p1, q1 = p1, q1, p0 + a*p1, q2 

140 n, d = d, n - a*d 

141 

142 k = (ctx.max_denom - q0)//q1 

143 

144 number = mpq(p, q) 

145 bound1 = mpq(p0 + k*p1, q0 + k*q1) 

146 bound2 = mpq(p1, q1) 

147 

148 if not bound2 or not bound1: 

149 return p, q 

150 elif abs(bound2 - number) <= abs(bound1 - number): 

151 return bound2._mpq_ 

152 else: 

153 return bound1._mpq_ 

154 

155 def almosteq(ctx, s, t, rel_eps=None, abs_eps=None): 

156 t = ctx.convert(t) 

157 if abs_eps is None and rel_eps is None: 

158 rel_eps = abs_eps = ctx.tolerance or ctx.make_tol() 

159 if abs_eps is None: 

160 abs_eps = ctx.convert(rel_eps) 

161 elif rel_eps is None: 

162 rel_eps = ctx.convert(abs_eps) 

163 diff = abs(s-t) 

164 if diff <= abs_eps: 

165 return True 

166 abss = abs(s) 

167 abst = abs(t) 

168 if abss < abst: 

169 err = diff/abst 

170 else: 

171 err = diff/abss 

172 return err <= rel_eps