Monomial comparison in SymPy with various monomial orderings

Question

Division of polynomials of more than one variable is a lot more complicated than one variable which we learned in high school.

I am trying to implement in SymPy the division algorithm for multi-variable polynomials given in pseudo-code form in theorem 3 of Chapter 2 of Ideals, Varieties, and Algorithms by David Cox, John Little, Donal O’Shea:

Input: f1 , . . . , fs , f
Output: a1 , . . . , as , r
a1 := 0; . . . ; as := 0; r := 0
p := f
WHILE p = 0 DO
    i := 1
    divisionoccurred := false
    WHILE i ≤ s AND divisionoccurred = false DO
        IF LT( f i ) divides LT( p) THEN
            ai := ai + LT( p)/LT( f i )
            p := p − (LT( p)/LT( f i )) f i
            divisionoccurred:= true
        ELSE
            i := i + 1
    IF divisionoccurred = false THEN
        r := r + LT( p)
        p := p − LT( p)

(Note that there is a typo in the book corrected in the errata. The book has the line:

    IF LT( f i ) divides ( p) THEN

which forgets the LT(p).)

This is the question: Given two monomials, m1 and m2, how to write a boolean function in SymPy deciding on whether m2 is divisible by m1. After searching the SymPy descriptions and the internet, I came up with the very awkward and klugey function:

from sympy import *

def is_div(m1,m2,variables):
    """ The inputs m1, m2 are monomials in sympy,
    e.g., x**3 * y**2. The input variables
    is a list of variables for the monomials
    for example [x,y] or [x1,x2,x3].
    
    The function will return True if m2 is divisible by
    the monomial m1. Otherwise it returns false.
    """
    
    # Here we make a tuple of the exponents of the 
    # monomial, e.g., x**3*y**2 -> (3,2)
    
    exp1 = Poly(m1,variables).monoms()[0]
    exp2 = Poly(m2,variables).monoms()[0]
    
    # The following compares all the entries of two tuples:
    
    if all(x <= y for x, y in zip(exp1, exp2)):
        return True
    else:
        return False 

print( is_div(x*y,x**2 * y**2,[x,y]) )
print( is_div(x*y**4,x**2 * y**2,[x,y]) )

I have searched the the internet in every permutation, but this basic question/procedure hasn't been addressed. I have not delved into the source code for the monoms() function.