Add new mjo.polynomial module with the "multidiv" function.

This new module implements the generalized polynomial division
algorithm in Section 2.3 of Cox, Little, and O'Shea.

diff --git a/mjo/all.py b/mjo/all.py
index 4962a15912fb10c2fc6997335ef91fa2ee607f25..9cdb23978eb1b7e25a7e2e957a136d820096081c 100644 (file)
--- a/mjo/all.py
+++ b/mjo/all.py
@@ -7,4 +7,5 @@ from mjo.cone.all import *
 from mjo.interpolation import *
 from mjo.misc import *
 from mjo.orthogonal_polynomials import *
+from mjo.polynomial import *
 from mjo.symbol_sequence import *

diff --git a/mjo/polynomial.py b/mjo/polynomial.py
new file mode 100644 (file)
index 0000000..3e90ba5
--- /dev/null
+++ b/mjo/polynomial.py
@@ -0,0 +1,95 @@
+from sage.all import *
+
+def multidiv(f, gs):
+    r"""
+    Divide the multivariate polynomial ``f`` by the ordered list of
+    polynomials ``gs`` from the same ring.
+
+    INPUT:
+
+    - ``f`` -- A multivariate polynomial, the "numerator," that we try
+      to express as a combination of the ``gs``.
+
+    - ``gs`` -- An ordered list of multipolynomials, the "denominators,"
+      in the same ring as ``f``.
+
+    OUTPUT:
+
+    A pair, the first component of which is the list of "quotients," and
+    the second of which is the remainder. All quotients and the
+    remainder will live in the same ring as ``f`` and the ``gs``. If the
+    ordered list of "quotients" is ``qs``, then ``f`` is the remainder
+    plus the sum of all ``qs[i]*gs[i]``.
+
+    The remainder is either zero, or a linear combination of monomials
+    none of which are divisible by any of the leading terms of the ``gs``.
+    Moreover if ``qs[i]*gs[i]`` is nonzero, then the multidegree of ``f``
+    is greater than or equal to the multidegree of ``qs[i]*gs[i]``.
+
+    SETUP::
+
+        sage: from mjo.polynomial import multidiv
+
+    ALGORITHM:
+
+    The algorithm from Theorem 3 in  Section 2.3 of Cox, Little, and O'Shea
+    is used almost verbatim.
+
+    REFERENCES:
+
+    - David A. Cox, John Little, Donal O'Shea. Ideals, Varieties, and
+      Algorithms: An Introduction to Computational Algebraic Geometry
+      and Commutative Algebra (Fourth Edition). Springer Undergraduate
+      Texts in Mathematics. Springer International Publishing, Switzerland,
+      2015. :doi:`10.1007/978-3-319-16721-3`.
+
+    EXAMPLES:
+
+    Example 1 in Section 2.3 of Cox, Little, and O'Shea::
+
+        sage: R = PolynomialRing(QQ, 'x,y', order='lex')
+        sage: x,y = R.gens()
+        sage: f = x*y^2 + 1
+        sage: gs = [ x*y + 1, y + 1 ]
+        sage: multidiv(f, gs)
+        ([y, -1], 2)
+
+    Example 2 in Section 2.3 of Cox, Little, and O'Shea::
+
+        sage: R = PolynomialRing(QQ, 'x,y', order='lex')
+        sage: x,y = R.gens()
+        sage: f = x^2*y + x*y^2 + y^2
+        sage: gs = [ x*y - 1, y^2 - 1 ]
+        sage: multidiv(f, gs)
+        ([x + y, 1], x + y + 1)
+
+    TESTS:
+
+    Derp.
+
+    """
+    R = f.parent()
+    s = len(gs)
+
+    p = f
+    r = R.zero()
+    qs = [R.zero()]*s
+
+    while p != R.zero():
+        for i in range(0,s):
+            division_occurred = false
+            # If gs[i].lt() divides p.lt(), then this remainder will
+            # be zero and the quotient will be in R (and not the
+            # fraction ring, which is important).
+            (factor, lt_r) = p.lt().quo_rem(gs[i].lt())
+            if lt_r.is_zero():
+                qs[i] += factor
+                p -= factor*gs[i]
+                division_occurred = true
+                break
+
+        if not division_occurred:
+            r += p.lt()
+            p -= p.lt()
+
+    return (qs,r)