from sage.all import *

def multidiv(f, gs):
    r"""
    Divide the multivariate polynomial ``f`` by the ordered list of
    polynomials ``gs`` from the same ring.

    INPUT:

    - ``f`` -- A multivariate polynomial, the "numerator," that we try
      to express as a combination of the ``gs``.

    - ``gs`` -- An ordered list of multipolynomials, the "denominators,"
      in the same ring as ``f``.

    OUTPUT:

    A pair, the first component of which is the list of "quotients," and
    the second of which is the remainder. All quotients and the
    remainder will live in the same ring as ``f`` and the ``gs``. If the
    ordered list of "quotients" is ``qs``, then ``f`` is the remainder
    plus the sum of all ``qs[i]*gs[i]``.

    The remainder is either zero, or a linear combination of monomials
    none of which are divisible by any of the leading terms of the ``gs``.
    Moreover if ``qs[i]*gs[i]`` is nonzero, then the multidegree of ``f``
    is greater than or equal to the multidegree of ``qs[i]*gs[i]``.

    SETUP::

        sage: from mjo.polynomial import multidiv

    ALGORITHM:

    The algorithm from Theorem 3 in  Section 2.3 of Cox, Little, and O'Shea
    is used almost verbatim.

    REFERENCES:

    - David A. Cox, John Little, Donal O'Shea. Ideals, Varieties, and
      Algorithms: An Introduction to Computational Algebraic Geometry
      and Commutative Algebra (Fourth Edition). Springer Undergraduate
      Texts in Mathematics. Springer International Publishing, Switzerland,
      2015. :doi:`10.1007/978-3-319-16721-3`.

    EXAMPLES:

    Example 1 in Section 2.3 of Cox, Little, and O'Shea::

        sage: R = PolynomialRing(QQ, 'x,y', order='lex')
        sage: x,y = R.gens()
        sage: f = x*y^2 + 1
        sage: gs = [ x*y + 1, y + 1 ]
        sage: (qs, r) = multidiv(f, gs)
        sage: (qs, r)
        ([y, -1], 2)
        sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
        True
        sage: not any( g.lt().divides(m) for m in r.monomials()
        ....:          for g in gs )
        True

    Example 2 in Section 2.3 of Cox, Little, and O'Shea::

        sage: R = PolynomialRing(QQ, 'x,y', order='lex')
        sage: x,y = R.gens()
        sage: f = x^2*y + x*y^2 + y^2
        sage: gs = [ x*y - 1, y^2 - 1 ]
        sage: (qs, r) = multidiv(f, gs)
        sage: (qs, r)
        ([x + y, 1], x + y + 1)
        sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
        True
        sage: not any( g.lt().divides(m) for m in r.monomials()
        ....:          for g in gs )
        True

    TESTS:

    If we get a zero remainder, then the numerator should belong to
    the ideal generated by the denominators::

        sage: set_random_seed()
        sage: R = PolynomialRing(QQ, 'x,y,z')
        sage: x,y,z = R.gens()
        sage: s = ZZ.random_element(1,5).abs()
        sage: gs = [ R.random_element() for idx in range(s) ]
        sage: f = R.random_element(ZZ.random_element(10).abs())
        sage: (qs, r) = multidiv(f,gs)
        sage: r != 0 or f in R.ideal(gs)
        True

    The numerator is always the sum of the remainder and the quotients
    times the denominators, and the remainder's monomials aren't divisible
    by the leading term of any denominator::

        sage: set_random_seed()
        sage: R = PolynomialRing(QQ, 'x,y,z')
        sage: x,y,z = R.gens()
        sage: s = ZZ.random_element(1,5).abs()
        sage: gs = [ R.random_element() for idx in range(s) ]
        sage: f = R.random_element(ZZ.random_element(10).abs())
        sage: (qs, r) = multidiv(f,gs)
        sage: r + sum( qs[i]*gs[i] for i in range(len(gs)) ) == f
        True
        sage: r == 0 or (not any( g.lt().divides(m) for m in r.monomials()
        ....:                     for g in gs ))
        True

    """
    R = f.parent()
    s = len(gs)

    p = f
    r = R.zero()
    qs = [R.zero()]*s

    while p != R.zero():
        i = 0
        division_occurred = false
        while i < s:
            # If gs[i].lt() divides p.lt(), then this remainder will
            # be zero and the quotient will be in R (and not the
            # fraction ring, which is important).
            (factor, lt_r) = p.lt().quo_rem(gs[i].lt())
            if lt_r.is_zero():
                qs[i] += factor
                p -= factor*gs[i]
                division_occurred = true
                # Don't increment "i" here because we want to try
                # again with this "denominator" g[i]. We might
                # get another factor out of it, but we know that
                # we can't get another factor out of an *earlier*
                # denominator g[i-k] for some k.
            else:
                i += 1

        if not division_occurred:
            r += p.lt()
            p -= p.lt()

    return (qs,r)