/
ratsimp.py
213 lines (166 loc) · 7.14 KB
/
ratsimp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
from itertools import combinations_with_replacement
from ..core import Add, Dummy, Rational, symbols
from ..polys import (ComputationFailedError, Poly, cancel,
parallel_poly_from_expr, reduced)
from ..polys.monomials import Monomial
def ratsimp(expr):
"""
Put an expression over a common denominator, cancel and reduce.
Examples
========
>>> ratsimp(1/x + 1/y)
(x + y)/(x*y)
"""
f, g = cancel(expr).as_numer_denom()
try:
Q, r = reduced(f, [g], field=True, expand=False)
except ComputationFailedError:
return f/g
return Add(*Q) + cancel(r/g)
def ratsimpmodprime(expr, G, *gens, **args):
"""
Simplifies a rational expression ``expr`` modulo the prime ideal
generated by ``G``. ``G`` should be a Gröbner basis of the
ideal.
>>> eq = (x + y**5 + y)/(x - y)
>>> ratsimpmodprime(eq, [x*y**5 - x - y], x, y, order='lex')
(x**2 + x*y + x + y)/(x**2 - x*y)
If ``polynomial`` is False, the algorithm computes a rational
simplification which minimizes the sum of the total degrees of
the numerator and the denominator.
If ``polynomial`` is True, this function just brings numerator and
denominator into a canonical form. This is much faster, but has
potentially worse results.
References
==========
M. Monagan, R. Pearce, Rational Simplification Modulo a Polynomial
Ideal,
http://citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.163.6984
(specifically, the second algorithm)
"""
from ..matrices import Matrix, zeros
from ..solvers.solvers import minsolve_linear_system
quick = args.pop('quick', True)
polynomial = args.pop('polynomial', False)
# usual preparation of polynomials:
num, denom = cancel(expr).as_numer_denom()
polys, opt = parallel_poly_from_expr([num, denom] + G, *gens, **args)
domain = opt.domain
opt.domain = domain.field
# compute only once
leading_monomials = [g.LM(opt.order) for g in polys[2:]]
tested = set()
def staircase(n):
"""
Compute all monomials with degree less than ``n`` that are
not divisible by any element of ``leading_monomials``.
"""
if n == 0:
return [1]
S = []
for mi in combinations_with_replacement(range(len(opt.gens)), n):
m = [0]*len(opt.gens)
for i in mi:
m[i] += 1
if all(not lmg.divides(m) for lmg in leading_monomials):
S.append(m)
return [Monomial(s).as_expr(*opt.gens) for s in S] + staircase(n - 1)
def _ratsimpmodprime(a, b, allsol, N=0, D=0):
r"""
Computes a rational simplification of ``a/b`` which minimizes
the sum of the total degrees of the numerator and the denominator.
The algorithm proceeds by looking at ``a * d - b * c`` modulo
the ideal generated by ``G`` for some ``c`` and ``d`` with degree
less than ``a`` and ``b`` respectively.
The coefficients of ``c`` and ``d`` are indeterminates and thus
the coefficients of the normalform of ``a * d - b * c`` are
linear polynomials in these indeterminates.
If these linear polynomials, considered as system of
equations, have a nontrivial solution, then `\frac{a}{b}
\equiv \frac{c}{d}` modulo the ideal generated by ``G``. So,
by construction, the degree of ``c`` and ``d`` is less than
the degree of ``a`` and ``b``, so a simpler representation
has been found.
After a simpler representation has been found, the algorithm
tries to reduce the degree of the numerator and denominator
and returns the result afterwards.
As an extension, if quick=False, we look at all possible degrees such
that the total degree is less than *or equal to* the best current
solution. We retain a list of all solutions of minimal degree, and try
to find the best one at the end.
"""
c, d = a, b
steps = 0
maxdeg = a.total_degree() + b.total_degree()
if quick:
bound = maxdeg - 1
else:
bound = maxdeg
while N + D <= bound:
if (N, D) in tested:
break
tested.add((N, D))
M1 = staircase(N)
M2 = staircase(D)
Cs = symbols(f'c:{len(M1):d}', cls=Dummy)
Ds = symbols(f'd:{len(M2):d}', cls=Dummy)
ng = Cs + Ds
c_hat = Poly(
sum(Cs[i] * M1[i] for i in range(len(M1))), *(opt.gens + ng))
d_hat = Poly(
sum(Ds[i] * M2[i] for i in range(len(M2))), *(opt.gens + ng))
r = reduced(a * d_hat - b * c_hat, G, *(opt.gens + ng),
order=opt.order, polys=True)[1]
S = Poly(r, *opt.gens).coeffs()
Sv = Cs + Ds
Sm = zeros(len(S), len(Sv) + 1)
Sm[:, :-1] = Matrix([[a.diff(b) for b in Sv] for a in S])
sol = minsolve_linear_system(Sm, *Sv, quick=True)
if sol and not all(s == 0 for s in sol.values()):
c = c_hat.subs(sol)
d = d_hat.subs(sol)
# The "free" variables occuring before as parameters
# might still be in the substituted c, d, so set them
# to the value chosen before:
c = c.subs(dict(zip(Cs + Ds, [1] * (len(Cs) + len(Ds)))))
d = d.subs(dict(zip(Cs + Ds, [1] * (len(Cs) + len(Ds)))))
c = Poly(c, *opt.gens)
d = Poly(d, *opt.gens)
if d == 0:
raise ValueError('Ideal not prime?')
allsol.append((c_hat, d_hat, S, Cs + Ds))
if N + D != maxdeg:
allsol = [allsol[-1]]
break
steps += 1
N += 1
D += 1
if steps > 0:
c, d, allsol = _ratsimpmodprime(c, d, allsol, N, D - steps)
c, d, allsol = _ratsimpmodprime(c, d, allsol, N - steps, D)
return c, d, allsol
# preprocessing. this improves performance a bit when deg(num)
# and deg(denom) are large:
num = reduced(num, G, *opt.gens, order=opt.order)[1]
denom = reduced(denom, G, *opt.gens, order=opt.order)[1]
if polynomial:
return (num/denom).cancel()
c, d, allsol = _ratsimpmodprime(
Poly(num, *opt.gens, domain=opt.domain), Poly(denom, *opt.gens, domain=opt.domain), [])
if not quick and allsol:
# Looking for best minimal solution.
newsol = []
for c_hat, d_hat, S, ng in allsol:
Sm = zeros(len(S), len(ng) + 1)
Sm[:, :-1] = Matrix([[a.diff(b) for b in ng] for a in S])
sol = minsolve_linear_system(Sm, *ng)
newsol.append((c_hat.subs(sol), d_hat.subs(sol)))
c, d = min(newsol, key=lambda x: len(x[0].terms()) + len(x[1].terms()))
if not domain.is_Field:
cn, c = c.clear_denoms(convert=True)
dn, d = d.clear_denoms(convert=True)
else:
cn, dn = 1, 1
cf, c, d = cancel((c, d), *opt.gens, order=opt.order) # canonicalize signs
r = cf*Rational(cn, dn)
return (c*r.denominator)/(d*r.numerator)