Computing xgcd in Sage

File: /ext/sage/sage-8.9_1804/local/lib/python2.7/site-packages/sage/arith/misc.py
Signature : xgcd(a, b)
Docstring :
Return a triple "(g,s,t)" such that g = s * a+t * b =
gcd(a,b).

Note: One exception is if a and b are not in a principal ideal
  domain (see
  https://en.wikipedia.org/wiki/Principal_ideal_domain), e.g., they
  are both polynomials over the integers. Then this function can't
  in general return "(g,s,t)" as above, since they need not exist.
  Instead, over the integers, we first multiply g by a divisor of
  the resultant of a/g and b/g, up to sign.

INPUT:

* "a, b" - integers or more generally, element of a ring for
  which the xgcd make sense (e.g. a field or univariate
  polynomials).

OUTPUT:

* "g, s, t" - such that g = s * a + t * b

Note: There is no guarantee that the returned cofactors (s and t)
  are minimal.

EXAMPLES:

   sage: xgcd(56, 44)
   (4, 4, -5)
   sage: 4*56 + (-5)*44
   4

   sage: g, a, b = xgcd(5/1, 7/1); g, a, b
   (1, 3, -2)
   sage: a*(5/1) + b*(7/1) == g
   True

   sage: x = polygen(QQ)
   sage: xgcd(x^3 - 1, x^2 - 1)
   (x - 1, 1, -x)

   sage: K.<g> = NumberField(x^2-3)
   sage: g.xgcd(g+2)
   (1, 1/3*g, 0)

   sage: R.<a,b> = K[]
   sage: S.<y> = R.fraction_field()[]
   sage: xgcd(y^2, a*y+b)
   (1, a^2/b^2, ((-a)/b^2)*y + 1/b)
   sage: xgcd((b+g)*y^2, (a+g)*y+b)
   (1, (a^2 + (2*g)*a + 3)/(b^3 + (g)*b^2), ((-a + (-g))/b^2)*y + 1/b)

Here is an example of a xgcd for two polynomials over the integers,
where the linear combination is not the gcd but the gcd multiplied
by the resultant:

   sage: R.<x> = ZZ[]
   sage: gcd(2*x*(x-1), x^2)
   x
   sage: xgcd(2*x*(x-1), x^2)
   (2*x, -1, 2)
   sage: (2*(x-1)).resultant(x)
   2

Tests with numpy and gmpy2 types:

   sage: from numpy import int8
   sage: xgcd(4,int8(8))
   (4, 1, 0)
   sage: xgcd(int8(4),int8(8))
   (4, 1, 0)
   sage: from gmpy2 import mpz
   sage: xgcd(mpz(4), mpz(8))
   (4, 1, 0)
   sage: xgcd(4, mpz(8))
   (4, 1, 0)

   File: /ext/sage/sage-8.9_1804/local/lib/python2.7/site-packages/sage/arith/misc.py
   Source:
   def xgcd(a, b):
    r"""
    Return a triple ``(g,s,t)`` such that `g = s\cdot a+t\cdot b = \gcd(a,b)`.

    .. NOTE::

       One exception is if `a` and `b` are not in a principal ideal domain (see
       :wikipedia:`Principal_ideal_domain`), e.g., they are both polynomials
       over the integers. Then this function can't in general return ``(g,s,t)``
       as above, since they need not exist.  Instead, over the integers, we
       first multiply `g` by a divisor of the resultant of `a/g` and `b/g`, up
       to sign.

    INPUT:

    -  ``a, b`` - integers or more generally, element of a ring for which the
       xgcd make sense (e.g. a field or univariate polynomials).

    OUTPUT:

    -  ``g, s, t`` - such that `g = s\cdot a + t\cdot b`

    .. NOTE::

       There is no guarantee that the returned cofactors (s and t) are
       minimal.

    EXAMPLES::

        sage: xgcd(56, 44)
        (4, 4, -5)
        sage: 4*56 + (-5)*44
        4

        sage: g, a, b = xgcd(5/1, 7/1); g, a, b
        (1, 3, -2)
        sage: a*(5/1) + b*(7/1) == g
        True

        sage: x = polygen(QQ)
        sage: xgcd(x^3 - 1, x^2 - 1)
        (x - 1, 1, -x)

        sage: K.<g> = NumberField(x^2-3)
        sage: g.xgcd(g+2)
        (1, 1/3*g, 0)

        sage: R.<a,b> = K[]
        sage: S.<y> = R.fraction_field()[]
        sage: xgcd(y^2, a*y+b)
        (1, a^2/b^2, ((-a)/b^2)*y + 1/b)
        sage: xgcd((b+g)*y^2, (a+g)*y+b)
        (1, (a^2 + (2*g)*a + 3)/(b^3 + (g)*b^2), ((-a + (-g))/b^2)*y + 1/b)

    Here is an example of a xgcd for two polynomials over the integers, where the linear
    combination is not the gcd but the gcd multiplied by the resultant::

        sage: R.<x> = ZZ[]
        sage: gcd(2*x*(x-1), x^2)
        x
        sage: xgcd(2*x*(x-1), x^2)
        (2*x, -1, 2)
        sage: (2*(x-1)).resultant(x)
        2

    Tests with numpy and gmpy2 types::

        sage: from numpy import int8
        sage: xgcd(4,int8(8))
        (4, 1, 0)
        sage: xgcd(int8(4),int8(8))
        (4, 1, 0)
        sage: from gmpy2 import mpz
        sage: xgcd(mpz(4), mpz(8))
        (4, 1, 0)
        sage: xgcd(4, mpz(8))
        (4, 1, 0)

    TESTS:

    We check that :trac:`3330` has been fixed::

        sage: R.<a,b> = NumberField(x^2-3,'g').extension(x^2-7,'h')[]
        sage: h = R.base_ring().gen()
        sage: S.<y> = R.fraction_field()[]
        sage: xgcd(y^2, a*h*y+b)
        (1, 7*a^2/b^2, (((-h)*a)/b^2)*y + 1/b)
    """
    try:
        return a.xgcd(b)
    except AttributeError:
        a = py_scalar_to_element(a)
        b = py_scalar_to_element(b)
    except TypeError:
        b = py_scalar_to_element(b)
    return a.xgcd(b)

   File: /ext/sage/sage-8.9_1804/local/lib/python2.7/site-packages/sage/rings/integer.pyx
   Source:
       def xgcd(self, Integer n):
        r"""
        Return the extended gcd of this element and ``n``.

        INPUT:

        - ``n`` -- an integer

        OUTPUT:

        A triple ``(g, s, t)`` such that ``g`` is the non-negative gcd of
        ``self`` and ``n``, and ``s`` and ``t`` are cofactors satisfying the
        Bezout identity

        .. MATH::

            g = s \cdot \mathrm{self} + t \cdot n.

        .. NOTE::

            There is no guarantee that the cofactors will be minimal. If you
            need the cofactors to be minimal use :meth:`_xgcd`. Also, using
            :meth:`_xgcd` directly might be faster in some cases, see
            :trac:`13628`.

        EXAMPLES::

            sage: 6.xgcd(4)
            (2, 1, -1)

        """
        return self._xgcd(n)

   File: /ext/sage/sage-8.9_1804/local/lib/python2.7/site-packages/sage/rings/integer.pyx
   Source:
       def _xgcd(self, Integer n, bint minimal=0):
        r"""
        Return the extended gcd of ``self`` and ``n``.

        INPUT:

        - ``n`` -- an integer
        - ``minimal`` -- a boolean (default: ``False``), whether to compute
          minimal cofactors (see below)

        OUTPUT:

        A triple ``(g, s, t)`` such that ``g`` is the non-negative gcd of
        ``self`` and ``n``, and ``s`` and ``t`` are cofactors satisfying the
        Bezout identity

        .. MATH::

            g = s \cdot \mathrm{self} + t \cdot n.

        .. NOTE::

            If ``minimal`` is ``False``, then there is no guarantee that the
            returned cofactors will be minimal in any sense; the only guarantee
            is that the Bezout identity will be satisfied (see examples below).

            If ``minimal`` is ``True``, the cofactors will satisfy the following
            conditions. If either ``self`` or ``n`` are zero, the trivial
            solution is returned. If both ``self`` and ``n`` are nonzero, the
            function returns the unique solution such that `0 \leq s < |n|/g`
            (which then must also satisfy
            `0 \leq |t| \leq |\mbox{\rm self}|/g`).

        EXAMPLES::

            sage: 5._xgcd(7)
            (1, 3, -2)
            sage: 5*3 + 7*-2
            1
            sage: g,s,t = 58526524056._xgcd(101294172798)
            sage: g
            22544886
            sage: 58526524056 * s + 101294172798 * t
            22544886

        Try ``minimal`` option with various edge cases::

            sage: 5._xgcd(0, minimal=True)
            (5, 1, 0)
            sage: (-5)._xgcd(0, minimal=True)
            (5, -1, 0)
            sage: 0._xgcd(5, minimal=True)
            (5, 0, 1)
            sage: 0._xgcd(-5, minimal=True)
            (5, 0, -1)
            sage: 0._xgcd(0, minimal=True)
            (0, 1, 0)

        Output may differ with and without the ``minimal`` option::

            sage: 5._xgcd(6)
            (1, -1, 1)
            sage: 5._xgcd(6, minimal=True)
            (1, 5, -4)

        Exhaustive tests, checking minimality conditions::

            sage: for a in srange(-20, 20):
            ....:   for b in srange(-20, 20):
            ....:     if a == 0 or b == 0: continue
            ....:     g, s, t = a._xgcd(b)
            ....:     assert g > 0
            ....:     assert a % g == 0 and b % g == 0
            ....:     assert a*s + b*t == g
            ....:     g, s, t = a._xgcd(b, minimal=True)
            ....:     assert g > 0
            ....:     assert a % g == 0 and b % g == 0
            ....:     assert a*s + b*t == g
            ....:     assert s >= 0 and s < abs(b)/g
            ....:     assert abs(t) <= abs(a)/g

        AUTHORS:

        - David Harvey (2007-12-26): added minimality option
        """
        cdef Integer g = PY_NEW(Integer)
        cdef Integer s = PY_NEW(Integer)
        cdef Integer t = PY_NEW(Integer)

        sig_on()
        mpz_gcdext(g.value, s.value, t.value, self.value, n.value)
        sig_off()

        # Note: the GMP documentation for mpz_gcdext (or mpn_gcdext for that
        # matter) makes absolutely no claims about any minimality conditions
        # satisfied by the returned cofactors. They guarantee a non-negative
        # gcd, but that's it. So we have to do some work ourselves.

        if not minimal:
            return g, s, t

        # handle degenerate cases n == 0 and self == 0

        if not mpz_sgn(n.value):
            mpz_set_ui(t.value, 0)
            mpz_abs(g.value, self.value)
            mpz_set_si(s.value, 1 if mpz_sgn(self.value) >= 0 else -1)
            return g, s, t

        if not mpz_sgn(self.value):
            mpz_set_ui(s.value, 0)
            mpz_abs(g.value, n.value)
            mpz_set_si(t.value, 1 if mpz_sgn(n.value) >= 0 else -1)
            return g, s, t

        # both n and self are nonzero, so we need to do a division and
        # make the appropriate adjustment

        cdef mpz_t u1, u2
        mpz_init(u1)
        mpz_init(u2)
        mpz_divexact(u1, n.value, g.value)
        mpz_divexact(u2, self.value, g.value)
        if mpz_sgn(u1) > 0:
            mpz_fdiv_qr(u1, s.value, s.value, u1)
        else:
            mpz_cdiv_qr(u1, s.value, s.value, u1)
        mpz_addmul(t.value, u1, u2)
        mpz_clear(u2)
        mpz_clear(u1)

        return g, s, t