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?
I found a good alternative: reduced(x**2+y, [2*y]) gives the desired ([1/2], x^2).