mjo/polynomial.py: improve multidiv performance a bit.

[sage.d.git] / mjo / polynomial.py

diff --git a/mjo/polynomial.py b/mjo/polynomial.py
index 3e90ba5bf3cf3862360ed5d5e5e3ccdcd9c17647..7414bdf16b33e8178c746a3edf5b6fa4dd3f985b 100644 (file)
--- a/mjo/polynomial.py
+++ b/mjo/polynomial.py
@@ -51,8 +51,14 @@ def multidiv(f, gs):
         sage: x,y = R.gens()
         sage: f = x*y^2 + 1
         sage: gs = [ x*y + 1, y + 1 ]
-        sage: multidiv(f, gs)
+        sage: (qs, r) = multidiv(f, gs)
+        sage: (qs, r)
         ([y, -1], 2)
+        sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
+        True
+        sage: not any( g.lt().divides(m) for m in r.monomials()
+        ....:          for g in gs )
+        True
 
     Example 2 in Section 2.3 of Cox, Little, and O'Shea::
 
@@ -60,12 +66,46 @@ def multidiv(f, gs):
         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)
+        sage: (qs, r) = multidiv(f, gs)
+        sage: (qs, r)
         ([x + y, 1], x + y + 1)
+        sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
+        True
+        sage: not any( g.lt().divides(m) for m in r.monomials()
+        ....:          for g in gs )
+        True
 
     TESTS:
 
-    Derp.
+    If we get a zero remainder, then the numerator should belong to
+    the ideal generated by the denominators::
+
+        sage: set_random_seed()
+        sage: R = PolynomialRing(QQ, 'x,y,z')
+        sage: x,y,z = R.gens()
+        sage: s = ZZ.random_element(1,5).abs()
+        sage: gs = [ R.random_element() for idx in range(s) ]
+        sage: f = R.random_element(ZZ.random_element(10).abs())
+        sage: (qs, r) = multidiv(f,gs)
+        sage: r != 0 or f in R.ideal(gs)
+        True
+
+    The numerator is always the sum of the remainder and the quotients
+    times the denominators, and the remainder's monomials aren't divisible
+    by the leading term of any denominator::
+
+        sage: set_random_seed()
+        sage: R = PolynomialRing(QQ, 'x,y,z')
+        sage: x,y,z = R.gens()
+        sage: s = ZZ.random_element(1,5).abs()
+        sage: gs = [ R.random_element() for idx in range(s) ]
+        sage: f = R.random_element(ZZ.random_element(10).abs())
+        sage: (qs, r) = multidiv(f,gs)
+        sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
+        True
+        sage: r == 0 or (not any( g.lt().divides(m) for m in r.monomials()
+        ....:                     for g in gs ))
+        True
 
     """
     R = f.parent()
@@ -76,8 +116,9 @@ def multidiv(f, gs):
     qs = [R.zero()]*s
 
     while p != R.zero():
-        for i in range(0,s):
-            division_occurred = false
+        i = 0
+        division_occurred = false
+        while i < s:
             # 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).
@@ -86,7 +127,13 @@ def multidiv(f, gs):
                 qs[i] += factor
                 p -= factor*gs[i]
                 division_occurred = true
-                break
+                # Don't increment "i" here because we want to try
+                # again with this "denominator" g[i]. We might
+                # get another factor out of it, but we know that
+                # we can't get another factor out of an *earlier*
+                # denominator g[i-k] for some k.
+            else:
+                i += 1
 
         if not division_occurred:
             r += p.lt()