Source code for diofant.simplify.ratsimp

from itertools import combinations_with_replacement

from ..core import Add, Dummy, Rational, symbols
from ..polys import (ComputationFailed, Poly, cancel, parallel_poly_from_expr,
from ..polys.monomials import Monomial
from ..utilities.misc import debug

[docs]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 ComputationFailed: 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, (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) debug('ratsimpmodprime', expr) # 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) debug('%s / %s: %s, %s' % (N, D, M1, M2)) Cs = symbols("c:%d" % len(M1), cls=Dummy) Ds = symbols("d:%d" % len(M2), 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, gens=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: debug('Looking for best minimal solution. Got: %s' % len(allsol)) 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)