# Ticket #12417: trac_12417-pfd.3.patch

File trac_12417-pfd.3.patch, 41.8 KB (added by araichev, 9 years ago)

Tested on Sage 5.0.1

• ## sage/categories/quotient_fields.py

# HG changeset patch
# User Alex Raichev <alex.raichev@gmail.com>
# Date 1340772096 -43200
# Node ID 8b83d9f8c0a1dcfeb427f74e3ec20c8b7d17aa3e
# Parent  01e0779db32182ef38b6bbeec105ba5f821a40a3
#12417: Refactored quotient_fields.py so that all partial fraction decompositions go through the REFD class

diff --git a/sage/categories/quotient_fields.py b/sage/categories/quotient_fields.py
 a r""" Quotient fields REFERENCES: .. [Raic2012]  Alexander Raichev. "Leinartas's partial fraction decomposition". _. """ #***************************************************************************** #  Copyright (C) 2008 Teresa Gomez-Diaz (CNRS) from sage.categories.category_singleton import Category_singleton from sage.misc.cachefunc import cached_method from sage.misc.abstract_method import abstract_method from functools import total_ordering class QuotientFields(Category_singleton): """ """ Greatest common divisor NOTE: NOTE:: In a field, the greatest common divisor is not very informative, as it is only determined up to a unit. But in """ Least common multiple NOTE: NOTE:: In a field, the least common multiple is not very informative, as it is only determined up to a unit. But in return (self.numerator().factor(*args, **kwds) / self.denominator().factor(*args, **kwds)) def partial_fraction_decomposition(self): """ Decomposes fraction field element into a whole part and a list of fraction field elements over prime power denominators. The sum will be equal to the original fraction. AUTHORS: - Robert Bradshaw (2007-05-31) EXAMPLES:: sage: S. = QQ[] sage: q = 1/(t+1) + 2/(t+2) + 3/(t-3); q (6*t^2 + 4*t - 6)/(t^3 - 7*t - 6) sage: whole, parts = q.partial_fraction_decomposition(); parts [3/(t - 3), 1/(t + 1), 2/(t + 2)] sage: sum(parts) == q True sage: q = 1/(t^3+1) + 2/(t^2+2) + 3/(t-3)^5 sage: whole, parts = q.partial_fraction_decomposition(); parts [1/3/(t + 1), 3/(t^5 - 15*t^4 + 90*t^3 - 270*t^2 + 405*t - 243), (-1/3*t + 2/3)/(t^2 - t + 1), 2/(t^2 + 2)] sage: sum(parts) == q True We do the best we can over inexact fields:: sage: R. = RealField(20)[] sage: q = 1/(x^2 + x + 2)^2 + 1/(x-1); q (x^4 + 2.0000*x^3 + 5.0000*x^2 + 5.0000*x + 3.0000)/(x^5 + x^4 + 3.0000*x^3 - x^2 - 4.0000) sage: whole, parts = q.partial_fraction_decomposition(); parts [1.0000/(x - 1.0000), 1.0000/(x^4 + 2.0000*x^3 + 5.0000*x^2 + 4.0000*x + 4.0000)] sage: sum(parts) (x^4 + 2.0000*x^3 + 5.0000*x^2 + 5.0000*x + 3.0000)/(x^5 + x^4 + 3.0000*x^3 - x^2 - 4.0000) TESTS: We test partial fraction for irreducible denominators:: sage: R. = ZZ[] sage: q = x^2/(x-1) sage: q.partial_fraction_decomposition() (x + 1, [1/(x - 1)]) sage: q = x^10/(x-1)^5 sage: whole, parts = q.partial_fraction_decomposition() sage: whole + sum(parts) == q True And also over finite fields (see trac #6052, #9945):: sage: R. = GF(2)[] sage: q = (x+1)/(x^3+x+1) sage: q.partial_fraction_decomposition() (0, [(x + 1)/(x^3 + x + 1)]) sage: R. = GF(11)[] sage: q = x + 1 + 1/(x+1) + x^2/(x^3 + 2*x + 9) sage: q.partial_fraction_decomposition() (x + 1, [1/(x + 1), x^2/(x^3 + 2*x + 9)]) And even the rationals:: sage: (26/15).partial_fraction_decomposition() (1, [1/3, 2/5]) """ from sage.misc.misc import prod denom = self.denominator() whole, numer = self.numerator().quo_rem(denom) factors = denom.factor() if factors.unit() != 1: numer *= ~factors.unit() if len(factors) == 1: return whole, [numer/r**e for r,e in factors] if not self.parent().is_exact(): # factors not grouped in this case all = {} for r in factors: all[r[0]] = 0 for r in factors: all[r[0]] += r[1] factors = all.items() factors.sort() # for doctest consistency factors = [r**e for r,e in factors] parts = [] for d in factors: # note that the product below is non-empty, since the case # of only one factor has been dealt with above n = numer * prod([r for r in factors if r != d]).inverse_mod(d) % d # we know the inverse exists as the two are relatively prime parts.append(n/d) return whole, parts def derivative(self, *args): r""" The derivative of this rational function, with respect to variables return self.__class__(R, num, den, coerce=False, reduce=False) def partial_fraction_decomposition(self, kind=None, factor=True): r""" Return a partial fraction decomposition of self of kind kind. INPUT: - kind - (Optional; default=None) The string 'nullstellensatz' or 'algebraic_dependence', indicating what kind of multivariate decomposition to return in case self is multivariate. The default multivariate decomposition is a Leinartas decomposition. - factor - (Optional; default=True) If True, then the output will be a REFDSum instance (which features factored denominators). Note that factoring the denominator of self happens anyway as part of the algorithm. OUTPUT: A sum of fraction field elements or, if factor=True, a REFDSum instance that is a partial fraction decomposition of self. EXAMPLES:: No variables:: sage: f = 64/45 sage: print f.partial_fraction_decomposition() [(1, []), (2, [(3, 2)]), (1, [(5, 1)])] sage: print f.partial_fraction_decomposition(factor=False) [1, 2/9, 1/5] sage: print sum(f.partial_fraction_decomposition(factor=False)) == f True One variable:: sage: R. = PolynomialRing(QQ) sage: f = 1 + 2/x + 2*x/(x^2 - 4*x + 8) sage: print f (x^3 + 16)/(x^3 - 4*x^2 + 8*x) sage: d1 = f.partial_fraction_decomposition() sage: print d1 [(1, []), (2, [(x, 1)]), (2*x, [(x^2 - 4*x + 8, 1)])] sage: d1.sum().quotient() == f True sage: print f.partial_fraction_decomposition(factor=False) [1, 2/x, 2*x/(x^2 - 4*x + 8)] Two variables:: sage: R.= PolynomialRing(QQ) sage: f = 1/x + 1/y + 1/(x + y) sage: print f (x^2 + 3*x*y + y^2)/(x^2*y + x*y^2) sage: print f.partial_fraction_decomposition(kind='nullstellensatz') [(0, []), (x^2 + 3*x*y + y^2, [(y, 1), (x, 1), (x + y, 1)])] sage: print f.partial_fraction_decomposition(kind='algebraic_dependence') [(0, []), (x^2 + 3*x*y + y^2, [(y, 2), (x, 1)]), (-x^2 - 3*x*y - y^2, [(y, 2), (x + y, 1)])] sage: print f.partial_fraction_decomposition() [(0, []), (x^2 + 3*x*y + y^2, [(y, 2), (x, 1)]), (-x^2 - 3*x*y - y^2, [(y, 2), (x + y, 1)])] Three variables over a finite field:: sage: R.= PolynomialRing(GF(2)) sage: f = x + 1/x + 1/y + 1/z + 1/(x*y + z) sage: print f (x^3*y^2*z + x^2*y*z^2 + x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2)/(x^2*y^2*z + x*y*z^2) sage: d1 = f.partial_fraction_decomposition() sage: print d1 [(x, []), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (x*y + z, 1)]), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (y, 1), (x, 1)])] sage: d2 = f.partial_fraction_decomposition() sage: print d2 [(x, []), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (x*y + z, 1)]), (x^2*y^2 + x^2*y*z + x*y^2*z + x*z^2 + y*z^2, [(z, 2), (y, 1), (x, 1)])] sage: print d2.sum().quotient() == f True NOTE:: See the docstrings of the REFD class methods leinartas_decomposition(), nullstellensatz_decomposition(), and algebraic_dependence_decomposition() for definitions of the different kinds of partial fraction decomposition options available. AUTHORS: - Robert Bradshaw (2007-05-31) - Alexander Raichev (2012-06-25) """ f = REFD(quotient=self) if f.ring().ngens() <= 1: decomp = f.univariate_decomposition() else: if kind == 'nullstellensatz': decomp = f.nullstellensatz_decomposition() elif kind == 'algebraic_dependence': decomp = f.algebraic_dependence_decomposition() else: decomp = f.leinartas_decomposition() if factor: return decomp else: return [r.quotient() for r in decomp] @total_ordering class REFD(object): r""" Represents a rational expression with factored denominator (REFD) p/(q_1^{e_1} \cdots q_m^{e_m}) by storing the parts p and [(q_1,e_1), \ldots, (q_m,e_m)]. Here p, q_1, \ldots, q_m are elements of a 0- or 1-variate factorial polynomial ring R, q_1, \ldots, q_m are distinct irreducible elements of R, and e_1, \ldots, e_m are positive integers. An element r of R is represented as [r, (,)]. AUTHORS: - Alexander Raichev (2012-06-25) """ def __init__(self, numerator=None, denominator_list=None, quotient=None, reduce_=True): r""" Create a REFD instance. INPUT: - numerator - (Optional; default=None) An element p of a 0- or 1-variate factorial polynomial ring R. - denominator_list - (Optional; default=None) A list of the form [(q_1,e_1), \ldots, (q_m,e_m)] where the q_1, \ldots, q_m are distinct irreducible elements of R and the e_i are positive integers. - quotient - (Optional; default=None) An element of a field of fractions of a factorial ring. - reduce_ - (Optional; default=True) If True, then represent p/(q_1^{e_1} \cdots q_m^{e_m}) in lowest terms. If False, then won't attempt to divide p by any of the q_i. OUTPUT: A REFD instance representing the rational expression p/(q_1^{e_1} \cdots q_m^{e_m}). To get a non-None output, one of numerator or quotient must be non-None. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: f = REFD(64, [(3, 2), (5, 1)]) sage: print f (64, [(3, 2), (5, 1)]) :: sage: R. = PolynomialRing(QQ) sage: qs = [x, 1], [y, 1], [x*y+1, 1] sage: f = REFD(x, qs) sage: print f (1, [(y, 1), (x*y + 1, 1)]) sage: ff = REFD(x, qs, reduce_=False) sage: print ff (x, [(y, 1), (x, 1), (x*y + 1, 1)]) :: sage: f = REFD(x + y, [(x + y, 1)]) sage: print f (1, []) :: sage: f = 64/45 sage: print REFD(quotient=f) (64, [(3, 2), (5, 1)]) :: sage: R. = PolynomialRing(QQ) sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) sage: print REFD(quotient=f) (5*x^7 - 5*x^6 + 5/3*x^5 - 5/3*x^4 + 2*x^3 - 2/3*x^2 + 1/3*x - 1/3, [(x - 1, 1), (x, 1), (x^2 + 1/3, 1)]) :: sage: R. = PolynomialRing(QQ) sage: f = 2*y/(5*(x^3 - 1)*(y + 1)) sage: print REFD(quotient=f) (2/5*y, [(y + 1, 1), (x - 1, 1), (x^2 + x + 1, 1)]) Singular throws a 'not implemented' error when trying to factor in a multivariate polynomial ring over an inexact field :: sage: R.= PolynomialRing(CC) sage: f = (x + 1)/(x*y*(x*y + 1)^2) sage: REFD(quotient=f) Traceback (most recent call last): ... TypeError: Singular error: ? not implemented ? error occurred in or before STDIN line 17: def sage9=factorize(sage8); """ from sage.rings.polynomial.polynomial_ring import is_PolynomialRing from sage.rings.polynomial.multi_polynomial import is_MPolynomial if quotient is not None: p = quotient.numerator() q = quotient.denominator() R = p.parent() self._ring = R if (is_PolynomialRing(R) or is_MPolynomial(R)) and \ (q in R.base_ring()): # R is a polynomial ring and q is in its base ring. # No factoring needed. self._numerator = quotient self._denominator_list = [] else: # Otherwise factor q. try: qs = q.factor() except NotImplementedError: # Singular's factor() needs proof=False. qs = q.factor(proof=False) if qs.unit() != 1: p = p/qs.unit() qs = sorted([tuple(t) for t in qs]) # Sort for consitency. self._numerator = p self._denominator_list = qs # Done. No reducing needed as Sage automatically reduces fractions. return if numerator is not None: R = numerator.parent() else: R = None self._numerator = numerator self._denominator_list = sorted([tuple(t) for t in denominator_list]) self._ring = R if (numerator is not None) and reduce_: # Reduce rational expression if possible by canceling # irreducible factors in denominator_list. p = numerator qs = denominator_list p_new = p qs_new = [] for (q, e) in qs: a = 0 while q.divides(p_new) and a < e: a += 1 p_new = R(p_new/q) if e - a > 0: qs_new.append((q,e-a)) self._numerator = p_new self._denominator_list = qs_new def numerator(self): r""" Return the numerator of self. """ return self._numerator def denominator_list(self): r""" Return the denominator list of self. """ return self._denominator_list def denominator(self): r""" Return the denominator of self. """ from sage.misc.misc import prod return prod([q**e for q,e in self.denominator_list()]) def ring(self): r""" Return the ring of self. """ return self._ring def list(self): r""" Convert self into a list. """ return (self.numerator(), self.denominator_list()) def quotient(self): r""" Convert self into a field of fractions element. """ return self.numerator()/self.denominator() def __str__(self): return str(self.list()) def __eq__(self, other): r""" Two REFD instances are equal iff they represent the same rational expression. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: R.= PolynomialRing(QQ) sage: qs = [x,1], [y,1], [x*y+1,1] sage: f = REFD(x, qs) sage: ff = REFD(x, qs, reduce_=False) sage: f == ff True sage: g = REFD(y, qs) sage: g == f False """ return self.quotient() == other.quotient() def __ne__(self, other): r""" EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: R.= PolynomialRing(QQ) sage: qs = [x,1], [y,1], [x*y+1,1] sage: f = REFD(x, qs) sage: ff = REFD(x, qs, reduce_=False) sage: f != ff False sage: g = REFD(y, qs) sage: g != f True """ return not (self == other) def __lt__(self, other): r""" REFD A is less than REFD B iff the denominator list of A is shorter than that of B, or the denominator list lengths are equal and the denominator of A is less than that of B in the ambient ring, or the denominator list lengths are equal and the denominators are equal and the the numerator of A is less than that of B in the ambient ring. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: R.= PolynomialRing(QQ) sage: qs = [x,1], [y,1], [x*y+1,1] sage: f = REFD(x, qs) sage: ff = REFD(x, qs, reduce_=False) sage: g = REFD(x + 1, qs) sage: h = REFD(x + 2, qs) sage: print f < ff False sage: print f < g True sage: print g < f False sage: print g < h True """ sn = self.numerator() on = other.numerator() sdl = self.denominator_list() odl = other.denominator_list() sd = self.denominator() od = other.denominator() if self == other: return False elif len(sdl) < len(odl) or\ (len(sdl) == len(odl) and sd < od) or\ (len(sdl) == len(odl) and sd == od and sn < on): return True else: return False def univariate_decomposition(self): r""" Return the usual 0- or 1-variate partial fraction decomposition of self as a REFDSum instance. Assume that self lies in the field of fractions of a 0- or 1-variate factorial polynomial ring. EXAMPLES:: No variables:: sage: from sage.categories.quotient_fields import REFD sage: f = 64/45 sage: decomp = REFD(quotient=f).univariate_decomposition() sage: print decomp [(1, []), (2, [(3, 2)]), (1, [(5, 1)])] sage: print decomp.sum().quotient() == f True One variable:: sage: R. = PolynomialRing(QQ) sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) sage: print f (15*x^7 - 15*x^6 + 5*x^5 - 5*x^4 + 6*x^3 - 2*x^2 + x - 1)/(3*x^4 - 3*x^3 + x^2 - x) sage: decomp = REFD(quotient=f).univariate_decomposition() sage: print decomp [(5*x^3, []), (1, [(x - 1, 1)]), (1, [(x, 1)]), (1/3, [(x^2 + 1/3, 1)])] sage: print decomp.sum().quotient() == f True One variable over a finite field:: sage: R. = PolynomialRing(GF(2)) sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) sage: print f (x^6 + x^4 + 1)/(x^3 + x) sage: decomp = REFD(quotient=f).univariate_decomposition() sage: print decomp [(x^3, []), (1, [(x, 1)]), (x, [(x + 1, 2)])] sage: print decomp.sum().quotient() == f True One variable over an inexact field:: sage: R. = PolynomialRing(CC) sage: f = 5*x^3 + 1/x + 1/(x-1) + 1/(3*x^2 + 1) sage: print f (15.0000000000000*x^7 - 15.0000000000000*x^6 + 5.00000000000000*x^5 - 5.00000000000000*x^4 + 6.00000000000000*x^3 - 2.00000000000000*x^2 + x - 1.00000000000000)/(3.00000000000000*x^4 - 3.00000000000000*x^3 + x^2 - x) sage: decomp = REFD(quotient=f).univariate_decomposition() sage: print decomp [(5.00000000000000*x^3, []), (1.00000000000000, [(x - 1.00000000000000, 1)]), (-0.288675134594813*I, [(x - 0.577350269189626*I, 1)]), (1.00000000000000, [(x, 1)]), (0.288675134594813*I, [(x + 0.577350269189626*I, 1)])] sage: print decomp.sum().quotient() == f # Rounding error coming False NOTE:: Let f = p/q be a rational expression where p and q lie in a 0- or 1-variate factorial polynomial ring R. Let q_1^{e_1} \cdots q_m^{e_m} be the unique factorization of q in R into irreducible factors. Then f can be written uniquely as (*)   p_0 + \sum_{i=1}^{m} p_i/ q_i^{e_i}, for some p_j \in R. I call (*) the *usual partial fraction decomposition* of f. AUTHORS: - Robert Bradshaw (2007-05-31) - Alexander Raichev (2012-06-25) """ from sage.misc.misc import prod p = self.numerator() q = self.denominator() whole, p = p.quo_rem(q) qs = self.denominator_list() decomp = [REFD(whole, [])] for (a, m) in qs: numer = p * prod([b**n for (b, n) in qs if b != a])\ .inverse_mod(a**m) % (a**m) # The inverse exists because the product and a**m # are relatively prime. decomp.append(REFD(numer, [(a, m)])) return REFDSum(decomp) def nullstellensatz_certificate(self): r""" Let [(q_1, e_1), \ldots, (q_m,e_m)] be the denominator list of self. Return a list of polynomials h_1, \ldots, h_m in self.ring() that satisfies h_1 q_1 + \cdots + h_m q_m = 1 if it exists. Otherwise return None. Only works for multivariate self. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: R. = PolynomialRing(QQ) sage: f = 1/(x^2 * (x*y + 1)) sage: print f 1/(x^3*y + x^2) sage: ff = REFD(quotient=f) sage: L = ff.nullstellensatz_certificate() sage: print L [y^2, -x*y + 1] sage: qs = ff.denominator_list() sage: sum([L[i]*qs[i][0]**qs[i][1] for i in range(len(qs))]) == 1 True :: sage: f = 1/(x*y) sage: L = REFD(quotient=f).nullstellensatz_certificate() sage: L is None True """ R = self.ring() qs = self.denominator_list() J = R.ideal([q**e for q, e in qs]) if R(1) in J: return R(1).lift(J) else: return None def nullstellensatz_decomposition(self): r""" Return a Nullstellensatz decomposition of self as a REFDSum instance. Recursive. Only works for multivariate self. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: R. = PolynomialRing(QQ) sage: f = 1/(x*(x*y + 1)) sage: decomp = REFD(quotient=f).nullstellensatz_decomposition() sage: print decomp [(0, []), (1, [(x, 1)]), (-y, [(x*y + 1, 1)])] sage: decomp.sum().quotient() == f True sage: for r in decomp: ...       L = r.nullstellensatz_certificate() ...       L is None True True True NOTE:: Let f = p/q be a rational expression where p and q lie in a d-variate polynomial ring K[X] for some field K and d \ge 1. Let q_1^{e_1} \cdots q_m^{e_m} be the unique factorization of q in K[X] into irreducible factors and let V_i be the algebraic variety \{x\in L^d: q_i(x) = 0\} of q_i over the algebraic closure L of K. By [Raic2012]_, f can be written as (*)   \sum p_A/\prod_{i \in A} q_i^{e_i}, where the p_A are in K[X] and the sum is taken over all subsets A \subseteq \{1, \ldots, m\} such that \cap_{i\in A} T_i \neq \emptyset. I call (*) a *Nullstellensatz decomposition* of f. Nullstellensatz decompositions are not unique. The algorithm used comes from [Raic2012]_. """ decomp = REFDSum() L = self.nullstellensatz_certificate() if L is None: # No decomposing possible. decomp = REFDSum([self]) else: # Decompose recursively. p = self.numerator() qs = self.denominator_list() m = len(qs) iteration1 = REFDSum([REFD(p*L[i],[(qs[j][0], qs[j][1]) for j in range(m) if j != i]) for i in range(m) if L[i] != 0]) # Now decompose each REFD of iteration1. for r in iteration1: decomp.extend(r.nullstellensatz_decomposition()) # Simplify and return result. return decomp.combine_like_terms().whole_and_parts() def algebraic_dependence_certificate(self): r""" Return the ideal J of annihilating polynomials for the set of polynomials [q**e for (q,e) in self.denominator_list()], which could be the zero ideal. The ideal J lies in a polynomial ring over the field self.ring().base_ring() that has m = len(self.denominator_list()) indeterminates. EXAMPLES: sage: from sage.categories.quotient_fields import REFD sage: R. = PolynomialRing(QQ) sage: f = 1/(x^2 * (x*y + 1) * y^3) sage: ff = REFD(quotient=f) sage: J = ff.algebraic_dependence_certificate() sage: print J Ideal (1 - 6*T2 + 15*T2^2 - 20*T2^3 + 15*T2^4 - T0^2*T1^3 - 6*T2^5 + T2^6) of Multivariate Polynomial Ring in T0, T1, T2 over Rational Field sage: g = J.gens()[0] sage: qs = ff.denominator_list() sage: g(*(q**e for q, e in qs)) == 0 True :: sage: f = 1/(x^3 * y^2) sage: J = REFD(quotient=f).algebraic_dependence_certificate() sage: print J Ideal (0) of Multivariate Polynomial Ring in T0, T1 over Rational Field """ from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing R = self.ring() qs = self.denominator_list() if not qs: return R.ideal()    # The zero ideal. m = len(qs) F = R.base_ring() Xs = list(R.gens()) d = len(Xs) # Expand R by 2*m new variables. S = 'S' while S in [str(x) for x in Xs]: S = S + 'S' Ss = [S + str(i) for i in range(m)] T = 'T' while T in [str(x) for x in Xs]: T = T + 'T' Ts = [T + str(i) for i in range(m)] Vs = [str(x) for x in Xs] + Ss + Ts RR = PolynomialRing(F, Vs) Xs = RR.gens()[:d] Ss = RR.gens()[d: d + m] Ts = RR.gens()[d + m: d + 2 * m] # Compute the appropriate elimination ideal. J = RR.ideal([ Ss[j] - RR(qs[j][0]) for j in range(m)] +\ [ Ss[j]**qs[j][1] - Ts[j] for j in range(m)]) J = J.elimination_ideal(Xs + Ss) # Coerce J into the polynomial ring in the indeteminates Ts[m:]. # I choose the negdeglex order because i find it useful in my work. RRR = PolynomialRing(F, [str(t) for t in Ts], order='negdeglex') return RRR.ideal(J) def algebraic_dependence_decomposition(self, whole_and_parts=True): r""" Return an algebraic dependence decomposition of self as a REFDSum instance. Recursive. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: R. = PolynomialRing(QQ) sage: f = 1/(x^2 * (x*y + 1) * y^3) sage: decomp = REFD(quotient=f).algebraic_dependence_decomposition() sage: print decomp [(0, []), (-x, [(x*y + 1, 1)]), (x^2*y^2 - x*y + 1, [(y, 3), (x, 2)])] sage: print decomp.sum().quotient() == f True sage: for r in decomp: ...       J = r.algebraic_dependence_certificate() ...       Z = J.ring().ideal()    # The zero ideal ...       J == Z True True True NOTE:: Let f = p/q be a rational expression where p and q lie in a d-variate polynomial ring K[X] for some field K. Let q_1^{e_1} \cdots q_m^{e_m} be the unique factorization of q in K[X] into irreducible factors and let V_i be the algebraic variety \{x\in L^d: q_i(x) = 0\} of q_i over the algebraic closure L of K. By [Raic2012]_, f can be written as (*)   \sum p_A/\prod_{i \in A} q_i^{b_i}, where the b_i are positive integers, the p_A are in K[X], and the sum is taken over all subsets A \subseteq \{1, \ldots, m\} such that |A| \le d and \{q_i : i\in A\} is algebraically independent. I call (*) an *algebraic dependence decomposition* of f. Algebraic dependence decompositions are not unique. The algorithm used comes from [Raic2012]_. """ from sage.misc.misc import prod decomp = REFDSum() R = self.ring() J = self.algebraic_dependence_certificate() if not J: # No decomposing possible. decomp = REFDSum([self]) else: # Decompose recursively. p = self.numerator() qs = self.denominator_list() m = len(qs) g = J.gens()[0]     # An annihilating polynomial for the qs. new_vars = J.ring().gens() # Note that each new_vars[j] corresponds to qs[j] such that # g([q**e for q,e in qs]) = 0. # Assuming here that g.parent() has negdeglex term order # so that g.lt() is indeed the monomial we want below. # Use g to rewrite r into a sum of REFDs, # each with < m distinct denominator factors. gg = (g.lt() - g)/(g.lc()) numers = map(prod, zip(gg.coefficients(), gg.monomials())) e = list(g.lt().exponents())[0:m] denoms = [(new_vars[j], e[0][j] + 1) for j in range(m)] # Write r in terms of new_vars, # cancel factors in the denominator, and combine like terms. iteration1_temp = REFDSum([REFD(a, denoms) for a in numers]).\ combine_like_terms() # Substitute in qs. qpowsub = dict([(new_vars[j],qs[j][0]**qs[j][1]) for j in range(m)]) iteration1 = REFDSum() for r in iteration1_temp: num1 = p*g.parent()(r.numerator()).subs(qpowsub) denoms1 = [] for q,e in r.denominator_list(): j = new_vars.index(q) denoms1.append((qs[j][0], qs[j][1]*e)) iteration1.append(REFD(num1, denoms1)) # Now decompose each REFD of iteration1. for r in iteration1: decomp.extend(r.algebraic_dependence_decomposition()) # Simplify and return result. return decomp.combine_like_terms().whole_and_parts() def leinartas_decomposition(self): r""" Return a Leinartas decomposition of self as a REFDSum instance. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: R. = PolynomialRing(QQ) sage: f = 1/x + 1/y + 1/(x*y + 1) sage: decomp = REFD(quotient=f).leinartas_decomposition() sage: print decomp [(0, []), (1, [(x*y + 1, 1)]), (x + y, [(y, 1), (x, 1)])] sage: print decomp.sum().quotient() == f True sage: for r in decomp: ...       L = r.nullstellensatz_certificate() ...       print L is None ...       J = r.algebraic_dependence_certificate() ...       Z = J.ring().ideal()    # The zero ideal ...       print J == Z True True True True True True :: sage: R.= PolynomialRing(GF(2, 'a')) sage: f = 1/(x * y * z * (x*y + z)) sage: decomp = REFD(quotient=f).leinartas_decomposition() sage: print decomp [(0, []), (1, [(z, 2), (x*y + z, 1)]), (1, [(z, 2), (y, 1), (x, 1)])] sage: print decomp.sum().quotient() == f True NOTE:: Let f = p/q be a rational expression where p and q lie in a d-variate polynomial ring K[X] for some field K. Let q_1^{e_1} \cdots q_m^{e_m} be the unique factorization of q in K[X] into irreducible factors and let V_i be the algebraic variety \{x\in L^d: q_i(x) = 0\} of q_i over the algebraic closure L of K. By [Raic2012]_, f can be written as (*)   \sum p_A/\prod_{i \in A} q_i^{b_i}, where the b_i are positive integers, the p_A are in K[X], and the sum is taken over all subsets A \subseteq \{1, \ldots, m\} such that (1) |A| \le d, (2) \cap_{i\in A} T_i \neq \emptyset, and (3) \{q_i : i\in A\} is algebraically independent. In particular, any rational expression in d variables can be represented as a sum of rational expressions whose denominators each contain at most d distinct irreducible factors. I call (*) a *Leinartas decomposition* of f. Leinartas decompositions are not unique. The algorithm used comes from [Raic2012]_. """ d = len(self.ring().gens()) if d == 1: # Sage's lift() function doesn't work in # univariate polynomial rings. # So nullstellensatz_decomposition() won't work. # So use algebraic_dependence_decomposition(), # which is sufficient. temp = REFDSum([self]) else: temp = self.nullstellensatz_decomposition() decomp = REFDSum() for r in temp: decomp.extend(r.algebraic_dependence_decomposition()) # Simplify and return result. return decomp.combine_like_terms().whole_and_parts() class REFDSum(list): r""" A list representing the sum of REFD objects with distinct denominator lists. AUTHORS: - Alexander Raichev (2012-06-25) """ def __str__(self): return str([r.list() for r in self]) def __eq__(self, other): return sorted(self) == sorted(other) def __ne__(self, other): return not (self == other) def ring(self): r""" Return the ring of self. """ if self: return self[0].ring() else: return None def whole_and_parts(self): r""" Rewrite self as a REFDSum of a (possibly zero) polynomial REFD followed by reduced rational expression REFDs. Only useful for multivariate decompositions. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: from sage.categories.quotient_fields import REFDSum sage: R = PolynomialRing(QQ, 'x, y') sage: x, y = R.gens() sage: f = x**2 + 3*y + 1/x + 1/y sage: f = REFD(quotient=f) sage: print f (x^3*y + 3*x*y^2 + x + y, [(y, 1), (x, 1)]) sage: print REFDSum([f]).whole_and_parts() [(x^2 + 3*y, []), (x + y, [(y, 1), (x, 1)])] """ whole = 0 parts = [] R = self.ring() for r in self: # Since r has already passed through REFD.__init__()'s reducing # procedure, r is already in lowest terms. # Check if can write r as a mixed fraction: whole + fraction. p = r.numerator() q = r.denominator() if q == 1: # r is already a whole polynomial whole += p else: # r is a fraction num = R(r.numerator())  # Coerce into R in case its a constant. denom = r.denominator() a, b = num.quo_rem(denom) whole += a parts.append(REFD(b, r.denominator_list(), reduce_=False)) return REFDSum([REFD(whole, ())] + parts) def combine_like_terms(self): r""" Combine terms in self with the same denominator. Only useful for multivariate decompositions. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: from sage.categories.quotient_fields import REFDSum sage: R.= PolynomialRing(QQ) sage: f = REFD(quotient=1/(x * y * (x*y + 1))) sage: g = REFD(quotient=x/(x * y * (x*y + 1))) sage: s = REFDSum([f, g, f]) sage: t = s.combine_like_terms() sage: print s [(1, [(y, 1), (x, 1), (x*y + 1, 1)]), (1, [(y, 1), (x*y + 1, 1)]), (1, [(y, 1), (x, 1), (x*y + 1, 1)])] sage: print t [(1, [(y, 1), (x*y + 1, 1)]), (2, [(y, 1), (x, 1), (x*y + 1, 1)])] """ if not self: return self # Combine like terms. refds = sorted(self) new_refds = [] temp = refds[0] for f in refds[1:]: if  temp.denominator_list() == f.denominator_list(): # Add f to temp. num = temp.numerator() + f.numerator() temp = REFD(num, temp.denominator_list()) else: # Append temp to new_refds and update temp. new_refds.append(temp) temp = f new_refds.append(temp) return REFDSum(new_refds) def sum(self): r""" Return the sum of the REFDs in self as a REFD. EXAMPLES:: sage: from sage.categories.quotient_fields import REFD sage: from sage.categories.quotient_fields import REFDSum sage: R. = PolynomialRing(QQ) sage: qs = (x, 1), (y, 1), (x*y + 1, 1) sage: f = REFD(quotient=1/(x * y * (x*y + 1))) sage: g = REFD(quotient=x*y/(x * y * (x*y + 1))) sage: print REFDSum([f, g]).sum() (1, [(y, 1), (x, 1)]) """ if not self: return self # Compute the sum's numerator and denominator. summy = sum((f.quotient() for f in self)) num = summy.numerator() denom = summy.denominator() if denom == 1: # Done return REFD(num, []) # Compute the irreducible factors of denom. # Could use the factor() command here, but that might be slow. factors = [] for f in self: factors.extend([q for q,e in f.denominator_list()]) # Eliminate repeats from factors and sort factors = sorted(list(set(factors))) # The irreducible factors of denom lie in factors. # Use this fact to build a denominator list for denom. denom_list = [] for q in factors: e = 0 quo, rem = denom.quo_rem(q) while rem == 0: e += 1 denom = quo quo, rem = denom.quo_rem(q) if e > 0: denom_list.append((q, e)) return REFD(num, denom_list, reduce_=False) No newline at end of file