def METHOD_NAME(expr, G, *gens, **args):
    """
    Simplifies a rational expression ``expr`` modulo the prime ideal
    generated by ``G``.  ``G`` should be a Gröbner basis of the
    ideal.
    >>> eq = (x + y**5 + y)/(x - y)
    >>> ratsimpmodprime(eq, [x*y**5 - x - y], x, y, order='lex')
    (x**2 + x*y + x + y)/(x**2 - x*y)
    If ``polynomial`` is False, the algorithm computes a rational
    simplification which minimizes the sum of the total degrees of
    the numerator and the denominator.
    If ``polynomial`` is True, this function just brings numerator and
    denominator into a canonical form. This is much faster, but has
    potentially worse results.
    References
    ==========
    M. Monagan, R. Pearce, Rational Simplification Modulo a Polynomial
    Ideal,
    http://citeseer.ist.psu.edu/viewdoc/summary?doi=10.1.1.163.6984
    (specifically, the second algorithm)
    """
    from ..matrices import Matrix, zeros
    from ..solvers.solvers import minsolve_linear_system
    quick = args.pop('quick', True)
    polynomial = args.pop('polynomial', False)
    # usual preparation of polynomials:
    num, denom = cancel(expr).as_numer_denom()
    polys, opt = parallel_poly_from_expr([num, denom] + G, *gens, **args)
    domain = opt.domain
    opt.domain = domain.field
    # compute only once
    leading_monomials = [g.LM(opt.order) for g in polys[2:]]
    tested = set()
    def staircase(n):
        """
        Compute all monomials with degree less than ``n`` that are
        not divisible by any element of ``leading_monomials``.
        """
        if n == 0:
            return [1]
        S = []
        for mi in combinations_with_replacement(range(len(opt.gens)), n):
            m = [0]*len(opt.gens)
            for i in mi:
                m[i] += 1
            if all(not lmg.divides(m) for lmg in leading_monomials):
                S.append(m)
        return [Monomial(s).as_expr(*opt.gens) for s in S] + staircase(n - 1)
    def _ratsimpmodprime(a, b, allsol, N=0, D=0):
        r"""
        Computes a rational simplification of ``a/b`` which minimizes
        the sum of the total degrees of the numerator and the denominator.
        The algorithm proceeds by looking at ``a * d - b * c`` modulo
        the ideal generated by ``G`` for some ``c`` and ``d`` with degree
        less than ``a`` and ``b`` respectively.
        The coefficients of ``c`` and ``d`` are indeterminates and thus
        the coefficients of the normalform of ``a * d - b * c`` are
        linear polynomials in these indeterminates.
        If these linear polynomials, considered as system of
        equations, have a nontrivial solution, then `\frac{a}{b}
        \equiv \frac{c}{d}` modulo the ideal generated by ``G``. So,
        by construction, the degree of ``c`` and ``d`` is less than
        the degree of ``a`` and ``b``, so a simpler representation
        has been found.
        After a simpler representation has been found, the algorithm
        tries to reduce the degree of the numerator and denominator
        and returns the result afterwards.
        As an extension, if quick=False, we look at all possible degrees such
        that the total degree is less than *or equal to* the best current
        solution. We retain a list of all solutions of minimal degree, and try
        to find the best one at the end.
        """
        c, d = a, b
        steps = 0
        maxdeg = a.total_degree() + b.total_degree()
        if quick:
            bound = maxdeg - 1
        else:
            bound = maxdeg
        while N + D <= bound:
            if (N, D) in tested:
                break
            tested.add((N, D))
            M1 = staircase(N)
            M2 = staircase(D)
            Cs = symbols(f'c:{len(M1):d}', cls=Dummy)
            Ds = symbols(f'd:{len(M2):d}', cls=Dummy)
            ng = Cs + Ds
            c_hat = Poly(
                sum(Cs[i] * M1[i] for i in range(len(M1))), *(opt.gens + ng))
            d_hat = Poly(
                sum(Ds[i] * M2[i] for i in range(len(M2))), *(opt.gens + ng))
            r = reduced(a * d_hat - b * c_hat, G, *(opt.gens + ng),
                        order=opt.order, polys=True)[1]
            S = Poly(r, *opt.gens).coeffs()
            Sv = Cs + Ds
            Sm = zeros(len(S), len(Sv) + 1)
            Sm[:, :-1] = Matrix([[a.diff(b) for b in Sv] for a in S])
            sol = minsolve_linear_system(Sm, *Sv, quick=True)
            if sol and not all(s == 0 for s in sol.values()):
                c = c_hat.subs(sol)
                d = d_hat.subs(sol)
                # The "free" variables occuring before as parameters
                # might still be in the substituted c, d, so set them
                # to the value chosen before:
                c = c.subs(dict(zip(Cs + Ds, [1] * (len(Cs) + len(Ds)))))
                d = d.subs(dict(zip(Cs + Ds, [1] * (len(Cs) + len(Ds)))))
                c = Poly(c, *opt.gens)
                d = Poly(d, *opt.gens)
                if d == 0:
                    raise ValueError('Ideal not prime?')
                allsol.append((c_hat, d_hat, S, Cs + Ds))
                if N + D != maxdeg:
                    allsol = [allsol[-1]]
                break
            steps += 1
            N += 1
            D += 1
        if steps > 0:
            c, d, allsol = _ratsimpmodprime(c, d, allsol, N, D - steps)
            c, d, allsol = _ratsimpmodprime(c, d, allsol, N - steps, D)
        return c, d, allsol
    # preprocessing. this improves performance a bit when deg(num)
    # and deg(denom) are large:
    num = reduced(num, G, *opt.gens, order=opt.order)[1]
    denom = reduced(denom, G, *opt.gens, order=opt.order)[1]
    if polynomial:
        return (num/denom).cancel()
    c, d, allsol = _ratsimpmodprime(
        Poly(num, *opt.gens, domain=opt.domain), Poly(denom, *opt.gens, domain=opt.domain), [])
    if not quick and allsol:
        # Looking for best minimal solution.
        newsol = []
        for c_hat, d_hat, S, ng in allsol:
            Sm = zeros(len(S), len(ng) + 1)
            Sm[:, :-1] = Matrix([[a.diff(b) for b in ng] for a in S])
            sol = minsolve_linear_system(Sm, *ng)
            newsol.append((c_hat.subs(sol), d_hat.subs(sol)))
        c, d = min(newsol, key=lambda x: len(x[0].terms()) + len(x[1].terms()))
    if not domain.is_Field:
        cn, c = c.clear_denoms(convert=True)
        dn, d = d.clear_denoms(convert=True)
    else:
        cn, dn = 1, 1
    cf, c, d = cancel((c, d), *opt.gens, order=opt.order)  # canonicalize signs
    r = cf*Rational(cn, dn)
    return (c*r.denominator)/(d*r.numerator)