The SymPy function reduced is not an alternative to multivariate polynomial division, it is exactly the same thing. Now some people prefer to say it division only when the set of divisor polynomials are Grobner (so the remainder is unique). The div function is for univariate division. And the reduced is for multivariate, and you can choose your monomial order as well. reduced( f, F, vars, order = ord) where f is the polynomial that you want to divide, F is the (ordered) list/tuple of divisor polynomials, vars you just write your variables in the order of which is more important (which will be used in your monomial order), and ord is name of your monomial order. So for example you could do reduced( x ** 2 + y, [ 2 * y ], x, y, order = 'grlex' ). In your simple example using lexicographic or graded lexicographic won't make any change, but in other examples it can.

You can watch this short video https://youtu.be/PT3wCsV-8-I and just to let you know, you can use other orders that are not by default defined in SymPy, check this video https://youtu.be/-xIiJ1L9Ack.