 CoCalc Public Filessupport / 2020-01-27-110136-xgcd.sagews
Author: William A. Stein
Description: Jupyter notebook support/2015-06-04-141749-bokeh.ipynb

# Computing xgcd in Sage

 xgcd(10983, 981)

(3, 46, -515)
xgcd?

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)

xgcd??

   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)


10983.xgcd??

   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)


10983._xgcd??

   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

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_clear(u2)
mpz_clear(u1)

return g, s, t