gcdCoefficents improvements

M2/Macaulay2/d/actors3.d

@@ -2393,6 +2393,29 @@ gcd(x:Expr,y:Expr):Expr := (
 gcdfun(e:Expr):Expr := accumulate(plus0,plus1,gcd,e);
 setupfun("gcd0",gcdfun);
 
+gcdCoefficients(e:Expr):Expr := (
+    when e
+    is a:Sequence do (
+	if length(a) == 2 then (
+	    when a.0
+	    is x:ZZcell do (
+		when a.1
+		is y:ZZcell do (
+		    g := newZZmutable();
+		    s := newZZmutable();
+		    t := newZZmutable();
+		    Ccode(void, "mpz_gcdext(", g, ", ", s, ", ", t, ", ", x.v,
+			", ", y.v, ")");
+		    seq(
+			Expr(ZZcell(moveToZZandclear(g))),
+			Expr(ZZcell(moveToZZandclear(s))),
+			Expr(ZZcell(moveToZZandclear(t)))))
+		else WrongArgZZ(2))
+	    else WrongArgZZ(1))
+	else WrongNumArgs(2))
+    else WrongNumArgs(2));
+setupfun("gcdCoefficients0", gcdCoefficients);
+
 binomial(e:Expr):Expr := (
     when e
     is a:Sequence do (

M2/Macaulay2/m2/factor.m2

@@ -39,13 +39,15 @@ gcd(RingElement,RingElement) := RingElement => (r,s) -> (
      else notImplemented()))
 gcd RingElement := identity
 
+gcdCoefficients(ZZ, RingElement) := (r, s) -> gcdCoefficients(r_(ring s), s)
+gcdCoefficients(RingElement, ZZ) := (r, s) -> gcdCoefficients(r, s_(ring r))
 gcdCoefficients(RingElement,RingElement) := (f,g) -> (	    -- ??
      R := ring f;
      if R =!= ring g then error "expected elements of the same ring";
      if not isPolynomialRing R then error "expected a polynomial ring";
      if not isField coefficientRing R then error "expected a polynomial ring over a field";
      if numgens R > 1 then error "expected a polynomial ring in at most one variable";
-     toList apply(rawExtendedGCD(raw f, raw g), r -> new R from r))
+     apply(rawExtendedGCD(raw f, raw g), r -> new R from r))
 
 lcm(ZZ,RingElement) := (r,s) -> lcm(promote(abs r,ring s),s)
 lcm(RingElement,ZZ) := (r,s) -> lcm(promote(abs s,ring r),r)

M2/Macaulay2/m2/quotient.m2

@@ -21,18 +21,7 @@ saturate = method(Options => options quotient)
 -- moved to packages/Saturation.m2 in October 2020
 annihilator = method(Options => {Strategy => null}) -- Intersection or Quotient
 
-gcdCoefficients(ZZ,ZZ) := (a,b) -> (
-     m := {a,1,0};
-     n := {b,0,1};
-     if a<0 then m=-m;
-     if b<0 then n=-n;
-     if a>b then (k :=m;m=n;n=k);
-     while m#0 > 0 do (
-	  t := n#0 // m#0;
-	  n = n - apply(m,y -> t * y);
-	  (k=m;m=n;n=k);
-	  );
-     n);
+gcdCoefficients(ZZ,ZZ) := gcdCoefficients0
 
 mod = (i,n) -> i * 1_(ZZ/n)
 

M2/Macaulay2/packages/Divisor.m2

@@ -979,9 +979,8 @@ bezoutNumbers := (l1) -> (
 	if (mySize == 1) then ({1}) else (
 		l2 := take(l1, 2);
 		l3 := take(l1, -(mySize - 2));
-		temp := gcdCoefficients toSequence l2;
-		tempCoeff := take(temp, -2);
-		tempGCD := temp#0;
+		(tempGCD, r, s) := gcdCoefficients toSequence l2;
+		tempCoeff := {r, s};
 		recursiveList := bezoutNumbers(prepend(tempGCD, l3));
 		recursiveList2 := take(recursiveList, -(mySize-2));
 		newCoeff := recursiveList#0;

M2/Macaulay2/packages/Macaulay2Doc/doc_arithmetic.m2

@@ -72,15 +72,42 @@ document {
 	  "conjugate 3"
 	  }
      }
-document {
-     Key => {gcdCoefficients,(gcdCoefficients, RingElement, RingElement),(gcdCoefficients, ZZ, ZZ)},
-     Headline => "gcd with coefficients",
-     TT "gcdCoefficients(a,b)", " -- returns ", TT "{d,r,s}", " so that
-     ", TT"a*r + b*s", " is the greatest common divisor ", TT "d", " of ", TT "a", "
-     and ", TT "b", ".",
-     PARA{},
-     "Works for integers or elements of polynomial rings in onve variable.",
-     SeeAlso => "gcd"}
+
+doc ///
+  Key
+    gcdCoefficients
+    (gcdCoefficients, RingElement, RingElement)
+    (gcdCoefficients, RingElement, ZZ)
+    (gcdCoefficients, ZZ, RingElement)
+    (gcdCoefficients, ZZ, ZZ)
+  Headline
+    greatest common divisor with coefficients
+  Usage
+    gcdCoefficients(a, b)
+  Inputs
+    a:{RingElement, ZZ}
+    b:{RingElement, ZZ}
+  Description
+    Text
+      This returns a sequence of the form @CODE "(d, r, s)"@, where $d=\gcd(a, b)$
+      and $r$ and $s$ are the minimal Bézout coefficients satisfying the
+      equation $d = ar + bs$.
+
+      It works for integers or elements of polynomial rings in one variable.
+    Example
+      (d, r, s) = gcdCoefficients(46, 240)
+      gcd(46, 240)
+      46 * r + 240 * s
+      R = ZZ/2[x]
+      f = x^8 + x^4 + x^3 + x + 1
+      g = x^6 + x^4 + x + 1
+      (d, r, s) = gcdCoefficients(f, g)
+      gcd(f, g)
+      f * r + g * s
+  SeeAlso
+    gcd
+///
+
 document {
      Key => mod,
      Headline => "reduce modulo an integer",

M2/Macaulay2/packages/TriangularSets.m2

@@ -719,8 +719,7 @@ myHermite = (A,c,negativepivot,badcol) -> (
             assign(i1,A#i2,c#i2,{});
             assign(i2,A2,c2,b2);
         )else(
-            D := gcdCoefficients(x,y);
-            d:=D#0; r:=D#1; s:=D#2;
+            (d,r,s) := gcdCoefficients(x,y);
             (A1,c1,b1):=addRows(i1,r,i2,s,j);
             (A2,c2,b2)=addRows(i1,y//d,i2,-x//d,j);
             assign(i1,A1,c1,b1);

M2/Macaulay2/tests/engine/test-smith.m2

@@ -29,7 +29,7 @@ reduceColGCD = (p,i,j,col) -> (
 	  g := gcd(e,pivotelem);
 	  a := (-pivotelem)//g;
 	  b := e//g;
-	  (GCD,c,d) := toSequence gcdCoefficients(e,pivotelem);
+	  (GCD,c,d) := gcdCoefficients(e,pivotelem);
 	  (a,b,c,d) = (raw (a_ZZ), raw(b_ZZ), raw(c_ZZ), raw(d_ZZ));
 	  rawMatrixColumnOperation2(p,col,j,a,b,c,d,false);
      ))
@@ -41,7 +41,7 @@ reduceRowGCD = (p,i,j,row) -> (
 	  g := gcd(e,pivotelem);
 	  a := (-pivotelem)//g;
 	  b := e//g;
-	  (GCD,c,d) := toSequence gcdCoefficients(e,pivotelem);
+	  (GCD,c,d) := gcdCoefficients(e,pivotelem);
 	  (a,b,c,d) = (raw (a_ZZ), raw(b_ZZ), raw(c_ZZ), raw(d_ZZ));
 	  rawMatrixRowOperation2(p,row,i,a,b,c,d,false);
      ))