The problem is that in such rings the quotient and remainder are not well-defined in general.
Consider the following example: divide $f = x^2 y + xy^2 + y^2$ by $y^2-1$ in the ring $\mathbb{Q}[x,y]/\langle xy-1\rangle$.
On the one hand, $$f = (x+1)(y^2-1) + 2x + 1 + x(xy-1),$$ so $(x+1)$ is "a quotient" and $2x+1$ is "a remainder" (and $x(xy-1)$ is zero in this ring).
On the other hand, $$f = 1\cdot(y^2-1) + x + y + 1 + (x+y)(xy-1),$$ so $1$ is "a quotient" and $x+y+1$ is "a remainder" (and $(x+y)(xy-1)$ is zero in this ring).
The problem is basically that division in a polynomial ring by several polynomials ($y^2-1$ and $xy-1$ in the example) is a priori not well-defined. For one, the problem is underdetermined. Choosing a monomial ordering and an ordering of the polynomials to be divided by allows you to define an unambiguous (multivariate) division algorithm. Unfortunately its output will depend on the orderings chosen. Furthermore, while a remainder of zero will mean that the input belongs to the ideal generated by the polynomials to be divided by, the converse is not necessarily true. Groebner bases were invented to address these issues. See for example the book Ideals, Varieties, and Algorithms by Cox, Little, and O'Shea. In fact I took the example from there; I just (slightly artificially) rephrased it to be about a quotient ring. If you are only interested in "a remainder" (any), then you can use the division algorithm of Theorem 3 in that book. If you want more, you have to read more. Good luck.
SageMath can find the remainder using the multivariate division algorithm:
sage: R.<x,y> = PolynomialRing(QQ,2,order='lex')
sage: (x^2*y+x*y^2+y^2).reduce([y^2-1,x*y-1])
2*x + 1
sage: R.<x,y> = PolynomialRing(QQ,2,order='invlex')
sage: (x^2*y+x*y^2+y^2).reduce([y^2-1,x*y-1])
y + x + 1
To find the quotient(s), I'm not sure if there is built-in method in SageMath. But it's not difficult to implement the division algorithm in SageMath, because things like monomial orderings and getting the leading term/coefficient, etc. are already there. There is also a trick to get the quotient(s) by adding extra variables to track what happens in the reduction. It goes something like this:
sage: R.<x,y,c1,c2> = PolynomialRing(QQ,4,order='lex(2),lex(2)')
sage: rem = (x^2*y+x*y^2+y^2).reduce([y^2-1-c1,x*y-1-c2])
sage: rem.coefficient(c1), rem.coefficient(c2)
(x + 1, x)
sage: rem.subs({c1:0,c2:0})
2*x + 1
sage: R.<x,y,c1,c2> = PolynomialRing(QQ,4,order='invlex(2),invlex(2)')
sage: rem = (x^2*y+x*y^2+y^2).reduce([y^2-1-c1,x*y-1-c2])
sage: rem.coefficient(c1), rem.coefficient(c2)
(1, y + x)
sage: rem.subs({c1:0,c2:0})
y + x + 1
