mjo/polynomial.py

   1 from sage.all import *
   2
   3 def multidiv(f, gs):
   4     r"""
   5     Divide the multivariate polynomial ``f`` by the ordered list of
   6     polynomials ``gs`` from the same ring.
   7
   8     INPUT:
   9
  10     - ``f`` -- A multivariate polynomial, the "numerator," that we try
  11       to express as a combination of the ``gs``.
  12
  13     - ``gs`` -- An ordered list of multipolynomials, the "denominators,"
  14       in the same ring as ``f``.
  15
  16     OUTPUT:
  17
  18     A pair, the first component of which is the list of "quotients," and
  19     the second of which is the remainder. All quotients and the
  20     remainder will live in the same ring as ``f`` and the ``gs``. If the
  21     ordered list of "quotients" is ``qs``, then ``f`` is the remainder
  22     plus the sum of all ``qs[i]*gs[i]``.
  23
  24     The remainder is either zero, or a linear combination of monomials
  25     none of which are divisible by any of the leading terms of the ``gs``.
  26     Moreover if ``qs[i]*gs[i]`` is nonzero, then the multidegree of ``f``
  27     is greater than or equal to the multidegree of ``qs[i]*gs[i]``.
  28
  29     SETUP::
  30
  31         sage: from mjo.polynomial import multidiv
  32
  33     ALGORITHM:
  34
  35     The algorithm from Theorem 3 in  Section 2.3 of Cox, Little, and O'Shea
  36     is used almost verbatim.
  37
  38     REFERENCES:
  39
  40     - David A. Cox, John Little, Donal O'Shea. Ideals, Varieties, and
  41       Algorithms: An Introduction to Computational Algebraic Geometry
  42       and Commutative Algebra (Fourth Edition). Springer Undergraduate
  43       Texts in Mathematics. Springer International Publishing, Switzerland,
  44       2015. :doi:`10.1007/978-3-319-16721-3`.
  45
  46     EXAMPLES:
  47
  48     Example 1 in Section 2.3 of Cox, Little, and O'Shea::
  49
  50         sage: R = PolynomialRing(QQ, 'x,y', order='lex')
  51         sage: x,y = R.gens()
  52         sage: f = x*y^2 + 1
  53         sage: gs = [ x*y + 1, y + 1 ]
  54         sage: (qs, r) = multidiv(f, gs)
  55         sage: (qs, r)
  56         ([y, -1], 2)
  57         sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
  58         True
  59         sage: not any( g.lt().divides(m) for m in r.monomials()
  60         ....:          for g in gs )
  61         True
  62
  63     Example 2 in Section 2.3 of Cox, Little, and O'Shea::
  64
  65         sage: R = PolynomialRing(QQ, 'x,y', order='lex')
  66         sage: x,y = R.gens()
  67         sage: f = x^2*y + x*y^2 + y^2
  68         sage: gs = [ x*y - 1, y^2 - 1 ]
  69         sage: (qs, r) = multidiv(f, gs)
  70         sage: (qs, r)
  71         ([x + y, 1], x + y + 1)
  72         sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
  73         True
  74         sage: not any( g.lt().divides(m) for m in r.monomials()
  75         ....:          for g in gs )
  76         True
  77
  78     A solution ``g`` to Exercise 6 in Section 2.3 of Cox, Little, and
  79     O'Shea that lives in the ideal generated by ``f1`` and ``f2`` but
  80     which has nonzero remainder after division::
  81
  82         sage: R = PolynomialRing(QQ, 'x,y', order='deglex')
  83         sage: x,y = R.gens()
  84         sage: f1 = 2*x*y^2 - x
  85         sage: f2 = 3*x^2*y - y - 1
  86         sage: I = R.ideal(f1,f2)
  87         sage: g = 2*y*f2
  88         sage: g in I
  89         True
  90         sage: (qs,r) = multidiv(g,[f1,f2])
  91         sage: r.is_zero()
  92         False
  93
  94     Two solutions ``g`` to Exercise 7 in Section 2.3 of Cox, Little, and
  95     O'Shea that live in the ideal generated by ``f1``, ``f2``, and ``f3``
  96     but which have nonzero remainders after division::
  97
  98         sage: R = PolynomialRing(QQ, 'x,y,z', order='deglex')
  99         sage: x,y,z = R.gens()
 100         sage: f1 = x^4*y^2 - z
 101         sage: f2 = x^3*y^3 - 1
 102         sage: f3 = x^2*y^4 - 2*z
 103         sage: I = R.ideal(f1,f2,f3)
 104         sage: g = x^2*f3
 105         sage: g in I
 106         True
 107         sage: (qs, r) = multidiv(g, [f1,f2,f3])
 108         sage: r.is_zero()
 109         False
 110         sage: g = x*f2
 111         sage: g in I
 112         True
 113         sage: (qs, r) = multidiv(g, [f1,f2,f3])
 114         sage: r.is_zero()
 115         False
 116
 117     TESTS:
 118
 119     If we get a zero remainder, then the numerator should belong to
 120     the ideal generated by the denominators::
 121
 122         sage: set_random_seed()
 123         sage: R = PolynomialRing(QQ, 'x,y,z')
 124         sage: x,y,z = R.gens()
 125         sage: s = ZZ.random_element(1,5).abs()
 126         sage: gs = [ R.random_element() for idx in range(s) ]
 127         sage: f = R.random_element(ZZ.random_element(10).abs())
 128         sage: (qs, r) = multidiv(f,gs)
 129         sage: r != 0 or f in R.ideal(gs)
 130         True
 131
 132     The numerator is always the sum of the remainder and the quotients
 133     times the denominators, and the remainder's monomials aren't divisible
 134     by the leading term of any denominator::
 135
 136         sage: set_random_seed()
 137         sage: R = PolynomialRing(QQ, 'x,y,z')
 138         sage: x,y,z = R.gens()
 139         sage: s = ZZ.random_element(1,5).abs()
 140         sage: gs = [ R.random_element() for idx in range(s) ]
 141         sage: f = R.random_element(ZZ.random_element(10).abs())
 142         sage: (qs, r) = multidiv(f,gs)
 143         sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
 144         True
 145         sage: r == 0 or (not any( g.lt().divides(m) for m in r.monomials()
 146         ....:                     for g in gs ))
 147         True
 148
 149     """
 150     R = f.parent()
 151     s = len(gs)
 152
 153     p = f
 154     r = R.zero()
 155     qs = [R.zero()]*s
 156
 157     while p != R.zero():
 158         i = 0
 159         division_occurred = false
 160         while i < s:
 161             # If gs[i].lt() divides p.lt(), then this remainder will
 162             # be zero and the quotient will be in R (and not the
 163             # fraction ring, which is important).
 164             (factor, lt_r) = p.lt().quo_rem(gs[i].lt())
 165             if lt_r.is_zero():
 166                 qs[i] += factor
 167                 p -= factor*gs[i]
 168                 division_occurred = true
 169                 # Don't increment "i" here because we want to try
 170                 # again with this "denominator" g[i]. We might
 171                 # get another factor out of it, but we know that
 172                 # we can't get another factor out of an *earlier*
 173                 # denominator g[i-k] for some k.
 174             else:
 175                 i += 1
 176
 177         if not division_occurred:
 178             r += p.lt()
 179             p -= p.lt()
 180
 181     return (qs,r)