Coverage for /usr/lib/python3/dist-packages/sympy/simplify/ratsimp.py: 9%

1from itertools import combinations_with_replacement 

2from sympy.core import symbols, Add, Dummy 

3from sympy.core.numbers import Rational 

4from sympy.polys import cancel, ComputationFailed, parallel_poly_from_expr, reduced, Poly 

5from sympy.polys.monomials import Monomial, monomial_div 

6from sympy.polys.polyerrors import DomainError, PolificationFailed 

7from sympy.utilities.misc import debug, debugf 

8 

9def ratsimp(expr): 

10 """ 

11 Put an expression over a common denominator, cancel and reduce. 

12 

13 Examples 

14 ======== 

15 

16 >>> from sympy import ratsimp 

17 >>> from sympy.abc import x, y 

18 >>> ratsimp(1/x + 1/y) 

19 (x + y)/(x*y) 

20 """ 

21 

22 f, g = cancel(expr).as_numer_denom() 

23 try: 

24 Q, r = reduced(f, [g], field=True, expand=False) 

25 except ComputationFailed: 

26 return f/g 

27 

28 return Add(*Q) + cancel(r/g) 

29 

30 

31def ratsimpmodprime(expr, G, *gens, quick=True, polynomial=False, **args): 

32 """ 

33 Simplifies a rational expression ``expr`` modulo the prime ideal 

34 generated by ``G``. ``G`` should be a Groebner basis of the 

35 ideal. 

36 

37 Examples 

38 ======== 

39 

40 >>> from sympy.simplify.ratsimp import ratsimpmodprime 

41 >>> from sympy.abc import x, y 

42 >>> eq = (x + y**5 + y)/(x - y) 

43 >>> ratsimpmodprime(eq, [x*y**5 - x - y], x, y, order='lex') 

44 (-x**2 - x*y - x - y)/(-x**2 + x*y) 

45 

46 If ``polynomial`` is ``False``, the algorithm computes a rational 

47 simplification which minimizes the sum of the total degrees of 

48 the numerator and the denominator. 

49 

50 If ``polynomial`` is ``True``, this function just brings numerator and 

51 denominator into a canonical form. This is much faster, but has 

52 potentially worse results. 

53 

54 References 

55 ========== 

56 

57 .. [1] M. Monagan, R. Pearce, Rational Simplification Modulo a Polynomial 

58 Ideal, https://dl.acm.org/doi/pdf/10.1145/1145768.1145809 

59 (specifically, the second algorithm) 

60 """ 

61 from sympy.solvers.solvers import solve 

62 

63 debug('ratsimpmodprime', expr) 

64 

65 # usual preparation of polynomials: 

66 

67 num, denom = cancel(expr).as_numer_denom() 

68 

69 try: 

70 polys, opt = parallel_poly_from_expr([num, denom] + G, *gens, **args) 

71 except PolificationFailed: 

72 return expr 

73 

74 domain = opt.domain 

75 

76 if domain.has_assoc_Field: 

77 opt.domain = domain.get_field() 

78 else: 

79 raise DomainError( 

80 "Cannot compute rational simplification over %s" % domain) 

81 

82 # compute only once 

83 leading_monomials = [g.LM(opt.order) for g in polys[2:]] 

84 tested = set() 

85 

86 def staircase(n): 

87 """ 

88 Compute all monomials with degree less than ``n`` that are 

89 not divisible by any element of ``leading_monomials``. 

90 """ 

91 if n == 0: 

92 return [1] 

93 S = [] 

94 for mi in combinations_with_replacement(range(len(opt.gens)), n): 

95 m = [0]*len(opt.gens) 

96 for i in mi: 

97 m[i] += 1 

98 if all(monomial_div(m, lmg) is None for lmg in 

99 leading_monomials): 

100 S.append(m) 

101 

102 return [Monomial(s).as_expr(*opt.gens) for s in S] + staircase(n - 1) 

103 

104 def _ratsimpmodprime(a, b, allsol, N=0, D=0): 

105 r""" 

106 Computes a rational simplification of ``a/b`` which minimizes 

107 the sum of the total degrees of the numerator and the denominator. 

108 

109 Explanation 

110 =========== 

111 

112 The algorithm proceeds by looking at ``a * d - b * c`` modulo 

113 the ideal generated by ``G`` for some ``c`` and ``d`` with degree 

114 less than ``a`` and ``b`` respectively. 

115 The coefficients of ``c`` and ``d`` are indeterminates and thus 

116 the coefficients of the normalform of ``a * d - b * c`` are 

117 linear polynomials in these indeterminates. 

118 If these linear polynomials, considered as system of 

119 equations, have a nontrivial solution, then `\frac{a}{b} 

120 \equiv \frac{c}{d}` modulo the ideal generated by ``G``. So, 

121 by construction, the degree of ``c`` and ``d`` is less than 

122 the degree of ``a`` and ``b``, so a simpler representation 

123 has been found. 

124 After a simpler representation has been found, the algorithm 

125 tries to reduce the degree of the numerator and denominator 

126 and returns the result afterwards. 

127 

128 As an extension, if quick=False, we look at all possible degrees such 

129 that the total degree is less than *or equal to* the best current 

130 solution. We retain a list of all solutions of minimal degree, and try 

131 to find the best one at the end. 

132 """ 

133 c, d = a, b 

134 steps = 0 

135 

136 maxdeg = a.total_degree() + b.total_degree() 

137 if quick: 

138 bound = maxdeg - 1 

139 else: 

140 bound = maxdeg 

141 while N + D <= bound: 

142 if (N, D) in tested: 

143 break 

144 tested.add((N, D)) 

145 

146 M1 = staircase(N) 

147 M2 = staircase(D) 

148 debugf('%s / %s: %s, %s', (N, D, M1, M2)) 

149 

150 Cs = symbols("c:%d" % len(M1), cls=Dummy) 

151 Ds = symbols("d:%d" % len(M2), cls=Dummy) 

152 ng = Cs + Ds 

153 

154 c_hat = Poly( 

155 sum([Cs[i] * M1[i] for i in range(len(M1))]), opt.gens + ng) 

156 d_hat = Poly( 

157 sum([Ds[i] * M2[i] for i in range(len(M2))]), opt.gens + ng) 

158 

159 r = reduced(a * d_hat - b * c_hat, G, opt.gens + ng, 

160 order=opt.order, polys=True)[1] 

161 

162 S = Poly(r, gens=opt.gens).coeffs() 

163 sol = solve(S, Cs + Ds, particular=True, quick=True) 

164 

165 if sol and not all(s == 0 for s in sol.values()): 

166 c = c_hat.subs(sol) 

167 d = d_hat.subs(sol) 

168 

169 # The "free" variables occurring before as parameters 

170 # might still be in the substituted c, d, so set them 

171 # to the value chosen before: 

172 c = c.subs(dict(list(zip(Cs + Ds, [1] * (len(Cs) + len(Ds)))))) 

173 d = d.subs(dict(list(zip(Cs + Ds, [1] * (len(Cs) + len(Ds)))))) 

174 

175 c = Poly(c, opt.gens) 

176 d = Poly(d, opt.gens) 

177 if d == 0: 

178 raise ValueError('Ideal not prime?') 

179 

180 allsol.append((c_hat, d_hat, S, Cs + Ds)) 

181 if N + D != maxdeg: 

182 allsol = [allsol[-1]] 

183 

184 break 

185 

186 steps += 1 

187 N += 1 

188 D += 1 

189 

190 if steps > 0: 

191 c, d, allsol = _ratsimpmodprime(c, d, allsol, N, D - steps) 

192 c, d, allsol = _ratsimpmodprime(c, d, allsol, N - steps, D) 

193 

194 return c, d, allsol 

195 

196 # preprocessing. this improves performance a bit when deg(num) 

197 # and deg(denom) are large: 

198 num = reduced(num, G, opt.gens, order=opt.order)[1] 

199 denom = reduced(denom, G, opt.gens, order=opt.order)[1] 

200 

201 if polynomial: 

202 return (num/denom).cancel() 

203 

204 c, d, allsol = _ratsimpmodprime( 

205 Poly(num, opt.gens, domain=opt.domain), Poly(denom, opt.gens, domain=opt.domain), []) 

206 if not quick and allsol: 

207 debugf('Looking for best minimal solution. Got: %s', len(allsol)) 

208 newsol = [] 

209 for c_hat, d_hat, S, ng in allsol: 

210 sol = solve(S, ng, particular=True, quick=False) 

211 # all values of sol should be numbers; if not, solve is broken 

212 newsol.append((c_hat.subs(sol), d_hat.subs(sol))) 

213 c, d = min(newsol, key=lambda x: len(x[0].terms()) + len(x[1].terms())) 

214 

215 if not domain.is_Field: 

216 cn, c = c.clear_denoms(convert=True) 

217 dn, d = d.clear_denoms(convert=True) 

218 r = Rational(cn, dn) 

219 else: 

220 r = Rational(1) 

221 

222 return (c*r.q)/(d*r.p)