mini-gmp mpz_gcdext Bézout coefficients do not match documentation

Niels Möller nisse at lysator.liu.se
Sun Feb 18 10:10:41 CET 2024


marco.bodrato at tutanota.com writes:

> s==0 is not a special case, it's one of the cases with smaller |s|.

Right, s == 0 is not a special case in itself, what's special is the
case g = a = b, then we have too candidate cofactor pairs, s=1, t=0,
s'=0, t'= 1. And we have the asymmetric preference for the first.

> What about;
>
> /* Arrange so that |s| < |u| / 2g */
>   mpz_add (s1, s0, s1);
>   {
>     int cmp = mpz_cmpabs (s0, s1);
>     if (cmp >= 0)
>       {
>         mpz_swap (s0, s1);
>         mpz_sub (t1, t0, t1);
>         if (cmp != 0 || mpz_cmpabs (t0, t1) > 0)
>           mpz_swap (t0, t1);
>       }
>   }

I think swapping of s and t must go together? I've tried out this
similar variant 

  mpz_add (s1, s0, s1);
  mpz_sub (t1, t0, t1);
  cmp = mpz_cmpabs (s0, s1);
  if (cmp > 0 || (cmp == 0 && mpz_cmpabs (t0, t1) > 0))
    {
      mpz_swap (s0, s1);
      mpz_swap (t0, t1);
    }

Seems to work fine, and it's a rather clear way to express the "lexical"
preference that we want the cofactor pair with the smallest of |s|, |s'|, but
in case of a tie, we choose the pair with smallest |t|, |t'|.

Updated version of the patch appended below.

Regrads,
/Niels

diff -r 1c566525476a mini-gmp/mini-gmp.c
--- a/mini-gmp/mini-gmp.c	Wed Dec 27 19:35:36 2023 +0100
+++ b/mini-gmp/mini-gmp.c	Sun Feb 18 10:02:56 2024 +0100
@@ -2809,6 +2809,7 @@
   mpz_t tu, tv, s0, s1, t0, t1;
   mp_bitcnt_t uz, vz, gz;
   mp_bitcnt_t power;
+  int cmp;
 
   if (u->_mp_size == 0)
     {
@@ -2960,12 +2961,21 @@
       mpz_tdiv_q_2exp (t0, t0, 1);
     }
 
-  /* Arrange so that |s| < |u| / 2g */
+  /* Choose small cofactors (they should satify
+
+       |s| < |u| / 2g and |t| < |v| / 2g,
+
+     except for documented exceptions). Always choose the smallest s,
+     if there are two choices for s with same absolute value, choose
+     the one with smallest corresponding t (this asymmetryic condition
+     is needed to prefer s = 0, |t| = 1 when g = |a| = |b|). */
   mpz_add (s1, s0, s1);
-  if (mpz_cmpabs (s0, s1) > 0)
+  mpz_sub (t1, t0, t1);
+  cmp = mpz_cmpabs (s0, s1);
+  if (cmp > 0 || (cmp == 0 && mpz_cmpabs (t0, t1) > 0))
     {
       mpz_swap (s0, s1);
-      mpz_sub (t0, t0, t1);
+      mpz_swap (t0, t1);
     }
   if (u->_mp_size < 0)
     mpz_neg (s0, s0);
diff -r 1c566525476a mini-gmp/tests/t-gcd.c
--- a/mini-gmp/tests/t-gcd.c	Wed Dec 27 19:35:36 2023 +0100
+++ b/mini-gmp/tests/t-gcd.c	Sun Feb 18 10:02:56 2024 +0100
@@ -55,6 +55,10 @@
   if (mpz_sgn (g) <= 0)
     return 0;
 
+  /* Require that s==0 iff g==abs(b) */
+  if (!mpz_sgn (s) != !mpz_cmpabs (g, b))
+    goto fail;
+
   mpz_init (ta);
   mpz_init (tb);
   mpz_init (r);
@@ -79,20 +83,20 @@
   if (mpz_sgn (r) != 0)
     goto fail;
 
-  /* Require that 2 |s| < |b/g|, or |s| == 1. */
-  if (mpz_cmpabs_ui (s, 1) > 0)
+  /* Require that 2 |s| < |b/g|, or s == sgn(a) */
+  if (mpz_cmp_si (s, mpz_sgn (a)) != 0)
     {
       mpz_mul_2exp (r, s, 1);
-      if (mpz_cmpabs (r, tb) > 0)
+      if (mpz_cmpabs (r, tb) >= 0)
 	goto fail;
     }
 
-  /* Require that 2 |t| < |a/g| or |t| == 1*/
-  if (mpz_cmpabs_ui (t, 1) > 0)
+  /* Require that 2 |t| < |a/g| or t == sgn(b) */
+  if (mpz_cmp_si (t, mpz_sgn (b)) != 0)
     {
       mpz_mul_2exp (r, t, 1);
-      if (mpz_cmpabs (r, ta) > 0)
-	return 0;
+      if (mpz_cmpabs (r, ta) >= 0)
+	goto fail;
     }
 
   mpz_clear (ta);
@@ -102,17 +106,53 @@
   return 1;
 }
 
+static void test_one(const mpz_t a, const mpz_t b)
+{
+  mpz_t g, s, t;
+
+  mpz_init (g);
+  mpz_init (s);
+  mpz_init (t);
+
+  mpz_gcdext (g, s, t, a, b);
+  if (!gcdext_valid_p (a, b, g, s, t))
+    {
+      fprintf (stderr, "mpz_gcdext failed:\n");
+      dump ("a", a);
+      dump ("b", b);
+      dump ("g", g);
+      dump ("s", s);
+      dump ("t", t);
+      abort ();
+    }
+
+  mpz_gcd (s, a, b);
+  if (mpz_cmp (g, s))
+    {
+      fprintf (stderr, "mpz_gcd failed:\n");
+      dump ("a", a);
+      dump ("b", b);
+      dump ("r", g);
+      dump ("ref", s);
+      abort ();
+    }
+
+  mpz_clear (g);
+  mpz_clear (s);
+  mpz_clear (t);
+}
+
 void
 testmain (int argc, char **argv)
 {
   unsigned i;
-  mpz_t a, b, g, s, t;
+  mpz_t a, b, g, s;
+  int ai, bi;
 
   mpz_init (a);
   mpz_init (b);
   mpz_init (g);
   mpz_init (s);
-  mpz_init (t);
 
   for (i = 0; i < COUNT; i++)
     {
@@ -129,6 +169,15 @@
 	}
     }
 
+  /* Exhaustive test of small inputs */
+  for (ai = -30; ai <= 30; ai++)
+    for (bi = -30; bi <= 30; bi++)
+      {
+	mpz_set_si (a, ai);
+	mpz_set_si (b, bi);
+	test_one (a, b);
+      }
+
   for (i = 0; i < COUNT; i++)
     {
       unsigned flags;
@@ -147,32 +196,11 @@
       if (flags & 2)
 	mpz_neg (b, b);
 
-      mpz_gcdext (g, s, t, a, b);
-      if (!gcdext_valid_p (a, b, g, s, t))
-	{
-	  fprintf (stderr, "mpz_gcdext failed:\n");
-	  dump ("a", a);
-	  dump ("b", b);
-	  dump ("g", g);
-	  dump ("s", s);
-	  dump ("t", t);
-	  abort ();
-	}
+      test_one (a, b);
+    }
 
-      mpz_gcd (s, a, b);
-      if (mpz_cmp (g, s))
-	{
-	  fprintf (stderr, "mpz_gcd failed:\n");
-	  dump ("a", a);
-	  dump ("b", b);
-	  dump ("r", g);
-	  dump ("ref", s);
-	  abort ();
-	}
-    }
   mpz_clear (a);
   mpz_clear (b);
   mpz_clear (g);
   mpz_clear (s);
-  mpz_clear (t);
 }
diff -r 1c566525476a tests/mpz/t-gcd.c
--- a/tests/mpz/t-gcd.c	Wed Dec 27 19:35:36 2023 +0100
+++ b/tests/mpz/t-gcd.c	Sun Feb 18 10:02:56 2024 +0100
@@ -432,6 +432,10 @@
   if (mpz_sgn (g) <= 0)
     return 0;
 
+  /* Require that s==0 iff g==abs(b) */
+  if (!mpz_sgn (s) != !mpz_cmpabs (g, b))
+    return 0;
+
   mpz_tdiv_qr (temp1, temp3, a, g);
   if (mpz_sgn (temp3) != 0)
     return 0;
@@ -440,8 +444,8 @@
   if (mpz_sgn (temp3) != 0)
     return 0;
 
-  /* Require that 2 |s| < |b/g|, or |s| == 1. */
-  if (mpz_cmpabs_ui (s, 1) > 0)
+  /* Require that 2 |s| < |b/g|, or s == sgn(a) */
+  if (mpz_cmp_si (s, mpz_sgn (a)) != 0)
     {
       mpz_mul_2exp (temp3, s, 1);
       if (mpz_cmpabs (temp3, temp2) >= 0)
@@ -456,8 +460,8 @@
   if (mpz_sgn (temp3) != 0)
     return 0;
 
-  /* Require that 2 |t| < |a/g| or |t| == 1*/
-  if (mpz_cmpabs_ui (temp2, 1) > 0)
+  /* Require that 2 |t| < |a/g| or t == sgn(b) */
+  if (mpz_cmp_si (temp2, mpz_sgn (b)) != 0)
     {
       mpz_mul_2exp (temp2, temp2, 1);
       if (mpz_cmpabs (temp2, temp1) >= 0)


-- 
Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677.
Internet email is subject to wholesale government surveillance.


More information about the gmp-bugs mailing list