Using div() on multivariate polynomials in Sympy: incorrect remainder?

Question

I want to compute remainders in Python for multivariate polynomials and I found that div() from sympy should do the trick (I also need sympy for Gröbner computations). But the problem I keep finding is that it seems that div() only checks the leading term for division because

q, r = div(x**2 + y, 2*x)

gives r=y, while

q, r = div(x**2 + y, 2*y) 

gives r=x**2+y.

I want to do things like Ideal Membership, hence finding the remainder of some polynomial f on division by G={g_1,...,g_s}, where by above I now cannot rely upon div().

Whilst working with Sage I don't get this problem (using (x^2+y)%y gives x^2), but I'm more familier with Python and prefer to do this via Python.

Can someone please tell me if I'm doing something wrong? Or does someone know a better function to use for remainders?

Answer 1

I found a good alternative: reduced(x**2+y, [2*y]) gives the desired ([1/2], x^2).