I try to do a simple division of polynomials in two variables using sage. Unfortunately, I get an unexpected result, as is illustrated in the code below. I tried several different ways to instantiate the ring and its variables, but the result stays the same. Using the %-operator to get a rest of a division yields the same results. Any ideas?
R, (x, y) = PolynomialRing(RationalField(), 2, 'xy').objgens()
t = x^2*y^2 + x^2*y - y + 1
f1 = x*y^2 + x
f2 = x*y - y^3
(q, r) = t.quo_rem(f1)
print "quotient:", q, "reminder:", r
assert q == x and r == x^2*y-x^2-y+1
(q, r) = r.quo_rem(f2) # error, expected q = x, r = -x^2+x*y^3-y+1
print "quotient:", q, "reminder:", r
assert q == x and r == -x^2+x*y^3-y+1
Output:
quotient: x reminder: x^2*y - x^2 - y + 1
quotient: 0 reminder: x^2*y - x^2 - y + 1
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
<ipython-input-8-861e8a9d1112> in <module>()
10 (q, r) = r.quo_rem(f2) # error, expected q = x, r = -x^2+x*y^3-y+1
11 print "quotient:", q, "reminder:", r
---> 12 assert q == x and r == -x**Integer(2)+x*y**Integer(3)-y+Integer(1)
AssertionError:
The result of division of multivariable polynomials depends on the chosen order of monomials, as is explained in Wikipedia. For the present, it suffices to say that the quotient from dividing P by Q is nonzero if and only if P contains a monomial that is divisible by the leading monomial of Q with respect to the chosen order.
By default, polynomial rings in Sage use the degree-reverse lexicographic order (degrevlex for short), in which monomials are first sorted by total degree, and then lexicographically within each degree. You can check this by using R.term_order(). And here is the list of term orders.
In the degrevlex order, y^3 precedes x*y because it has higher total degree. And since x^2*y-x^2-y+1 does not have a monomial divisible by y^3, you get zero quotient.
Your expected results are consistent with the lexicographic (lex) order of monomials. So, it would seem that the fix is to prescribe lexicographic term order when constructing R:
R.<x,y> = PolynomialRing(QQ, 2, 'xy', order='lex')
Unfortunately, quo_rem ignores the term order of R and still uses degrevlex. This appears to be a bug, either in the way Sage communicates with Singular, or in Singular itself (the Quotient command is deprecated).
Meanwhile, there is a workaround: use the reduce command instead of quo_rem, as it respects the term order of R. For example:
R.<x,y> = PolynomialRing(QQ, 2, 'xy', order='lex')
a = x^2*y-x^2-y+1
f2 = x*y-y^3
r = a.reduce(Ideal([f2]))
q = (a-r)/f2
This returns r = -x^2 + x*y^3 - y + 1 and q = x.