[Gmp-commit] /var/hg/gmp: mini-gmp: Fix bug in gcdext canonicalization, and s...
mercurial at gmplib.org
mercurial at gmplib.org
Sun Feb 18 20:22:42 CET 2024
details: /var/hg/gmp/rev/7ecb3b2beea1
changeset: 18465:7ecb3b2beea1
user: Niels Möller <nisse at lysator.liu.se>
date: Sun Feb 18 20:20:57 2024 +0100
description:
mini-gmp: Fix bug in gcdext canonicalization, and strengthen related tests.
diffstat:
ChangeLog | 6 +++
mini-gmp/ChangeLog | 10 +++++
mini-gmp/mini-gmp.c | 16 +++++++-
mini-gmp/tests/t-gcd.c | 92 ++++++++++++++++++++++++++++++++-----------------
tests/mpz/t-gcd.c | 12 ++++--
5 files changed, 97 insertions(+), 39 deletions(-)
diffs (251 lines):
diff -r 1c566525476a -r 7ecb3b2beea1 ChangeLog
--- a/ChangeLog Wed Dec 27 19:35:36 2023 +0100
+++ b/ChangeLog Sun Feb 18 20:20:57 2024 +0100
@@ -1,3 +1,9 @@
+2024-02-18 Niels Möller <nisse at lysator.liu.se>
+
+ * tests/mpz/t-gcd.c (gcdext_valid_p): Stricter check that cofactors
+ match the documented behavior. Adapted from the corresponding
+ mini-gmp test.
+
2023-12-27 Niels Möller <nisse at lysator.liu.se>
* tune/tuneup.c (tune_mul): Updated to measure all
diff -r 1c566525476a -r 7ecb3b2beea1 mini-gmp/ChangeLog
--- a/mini-gmp/ChangeLog Wed Dec 27 19:35:36 2023 +0100
+++ b/mini-gmp/ChangeLog Sun Feb 18 20:20:57 2024 +0100
@@ -1,3 +1,13 @@
+2024-02-18 Niels Möller <nisse at lysator.liu.se>
+
+ * mini-gmp.c (mpz_gcdext): Fix logic to canonicalize cofactors.
+ But reported by Guido Vranken.
+
+ * tests/t-gcd.c (gcdext_valid_p): Stricter check that cofactors
+ match the documented behavior.
+ (test_one): New function, extracted from testmain.
+ (testmain): Add exhaustive tests for all inputs |a|, |b| <= 30.
+
2023-07-20 Niels Möller <nisse at lysator.liu.se>
* mini-gmp.c (gmp_umullo_limb): New macro, to avoid unintended
diff -r 1c566525476a -r 7ecb3b2beea1 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 20:20:57 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 generally satify
+
+ |s| < |u| / 2g and |t| < |v| / 2g,
+
+ with some 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 asymmetric 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 -r 7ecb3b2beea1 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 20:20:57 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 -r 7ecb3b2beea1 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 20:20:57 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)
More information about the gmp-commit
mailing list