I have a working example shown below to help explain my question.
R.<x> = GF(2)[];
I = R.ideal([x^2 - 1])
S.<u> = R.quotient_ring(I)
f = x^2 + x
g = x
print(f.quo_rem(g))
f = u^2 + u
g = u
print(f.quo_rem(g))
The first quo_rem call works fine, but the second does not. Is polynomial division not supported for quotient rings, or am I calling it the wrong way? I'm specifically interested in doing this for multivariate polynomials, but I can't seem to get it to work for a single variable either.
Related to the question, I am interested in understanding this division in quotient rings better. If the operation isn't supported in sage, any sort of references/algorithms on how to do such a calculation would be very helpful and appreciated.