CoCalc Public Filessupport / 2020-01-27-110136-xgcd.sagewsOpen in with one click!
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 # 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