Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
I've confirmed with my fuzzer that the committed patch resolves the issues. Thanks. On Sun, Feb 18, 2024 at 8:27 PM Niels Möller wrote: > marco.bodr...@tutanota.com writes: > > > The only reason why I prefer my solution is: when cmp<0, there is no > need to compute > > mpz_sub (t1, t0, t1); > > That's certainly a micro optimization, but let's keep things simple. > > I've pushed this fix now, I think there were only comment and ChangeLog > changes since the previous version of the patch. > > Thanks for the review, > /Niels > > -- > Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677. > Internet email is subject to wholesale government surveillance. > ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
marco.bodr...@tutanota.com writes: > The only reason why I prefer my solution is: when cmp<0, there is no need to > compute > mpz_sub (t1, t0, t1); That's certainly a micro optimization, but let's keep things simple. I've pushed this fix now, I think there were only comment and ChangeLog changes since the previous version of the patch. Thanks for the review, /Niels -- Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677. Internet email is subject to wholesale government surveillance. ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
Ciao Niels, 18 feb 2024, 10:10 da ni...@lysator.liu.se: > marco.bodr...@tutanota.com writes: > >> 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 > You are right, it's better with the two swaps together. > 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); > } > The only reason why I prefer my solution is: when cmp<0, there is no need to compute mpz_sub (t1, 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'|. > I agree. Ĝis, m ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
marco.bodr...@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.cWed Dec 27 19:35:36 2023 +0100 +++ b/mini-gmp/tests/t-gcd.cSun 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",
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
Ciao Niels, 17 feb 2024, 21:26 da ni...@lysator.liu.se: > marco.bodr...@tutanota.com writes: > >> Maybe we should also add the test: >> >> /* Require that s==0 iff g==abs(b) */ >> if (!mpz_sgn (s) ^ !mpz_cmpabs (g, b)) >> goto fail; >> > > Good point, that's better than only checking the special case |a| == > |b|. (But maybe more readable with != instead of ^). > Yes, with != it's more readable! > To get mini-gmp to conform, I find no simpler way than special casing s > == 0, like > s==0 is not a special case, it's one of the cases with smaller |s|. 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); } } Ĝis, m ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
marco.bodr...@tutanota.com writes: > Maybe we should also add the test: > > /* Require that s==0 iff g==abs(b) */ > if (!mpz_sgn (s) ^ !mpz_cmpabs (g, b)) > goto fail; Good point, that's better than only checking the special case |a| == |b|. (But maybe more readable with != instead of ^). To get mini-gmp to conform, I find no simpler way than special casing s == 0, like if (mpz_sgn (s0) != 0) { /* Arrange so that |s| < |u| / 2g and |t| < |v| / 2g, if possible. */ mpz_add (s1, s0, s1); mpz_sub (t1, t0, t1); if (mpz_cmpabs (s0, s1) > 0 || mpz_cmpabs (t0, t1) > 0) { mpz_swap (s0, s1); mpz_swap (t0, t1); } } And I think this is the only condition that isn't fully symmetrical, so perhaps not surprising if it needs a special case. Regards, /Niels -- Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677. Internet email is subject to wholesale government surveillance. ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
Ciao, 17 feb 2024, 19:15 da : > Ciao Niels, > > 17 feb 2024, 17:53 da ni...@lysator.liu.se: > >> Niels Möller writes: >> The documented conditions say that gmp should return the same cofactors >> > > The documentation also says "If abs(a) = abs(b), then s = 0, t = sgn(b)." > Maybe we should also add the test: /* Require that s==0 iff g==abs(b) */ if (!mpz_sgn (s) ^ !mpz_cmpabs (g, b)) goto fail; Ĝis, m ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
marco.bodr...@tutanota.com writes: > And if I correctly patched and tested your proposed code. with equal numbers > I get t=0, instead of s=0. Thanks for testing, I have to look into that case, then. > Shat about simply changing the test from > to >= ? > > /* Arrange so that |s| < |u| / 2g */ > mpz_add (s1, s0, s1); > - if (mpz_cmpabs (s0, s1) > 0) > + if (mpz_cmpabs (s0, s1) >= 0) I don't think that's right. My understanding is that if |s| = |s'| (notation from my mail, the variables in the code are reassigned several times, which is a bit confusing), then we ncessarily have |t| != |t'|, and we must choose candidate to get the smaller one. Regards, /Niels -- Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677. Internet email is subject to wholesale government surveillance. ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
Ciao Niels, 17 feb 2024, 17:53 da ni...@lysator.liu.se: > Niels Möller writes: > The documented conditions say that gmp should return the same cofactors > The documentation also says "If abs(a) = abs(b), then s = 0, t = sgn(b)." > > canonicalization logic at the end wasn't right. Appending a patch with a > fix + stricter tests. > And if I correctly patched and tested your proposed code. with equal numbers I get t=0, instead of s=0. > Without the fix, we could make the wrong choice in case > > |s'| = |s| = |b| / 2g > > Since it's a bit subtle, it would be nice with a review of this code > before I commit it. > Shat about simply changing the test from > to >= ? /* Arrange so that |s| < |u| / 2g */ mpz_add (s1, s0, s1); - if (mpz_cmpabs (s0, s1) > 0) + if (mpz_cmpabs (s0, s1) >= 0) { mpz_swap (s0, s1); mpz_sub (t0, t0, t1); Ĝis, m ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
Niels Möller writes: > Guido Vranken writes: > >> Computing extended gcd using mpz_gcdext where a = 1, b = 2: >> >> libgmp: g: 1, s: 1, t: 0 >> mini-gmp: g: 1, s: -1, t: 1 >> This violates the docs: "s = sgn(a) if abs(b) = 2", e.g. s must be 1 >> >> Computing extended gcd using mpz_gcdext where a = 6, b = 4: >> libgmp: g: 2, s: 1, t: -1 >> mini-libgmp: g: 2, s: -1, t: 2 > > I think mini-gmp is wrong here, I'll have to investigate. Thanks for > reporting. The documented conditions say that gmp should return the same cofactors as one would get naturally with Euclid's algorithm. It appears mini-gmp uses binary gcdext (I had forgotten that detail), and the canonicalization logic at the end wasn't right. Appending a patch with a fix + stricter tests. At this place in the algorithm, we have a, b >= 0, and we have a candidate s, t, with -b / g <= s <= 0, 0 <= t <= a / g (not sure at which edges we might have equality), and then we construct another candidate s', t' as s' = s + b / g >= 0, t' = t - a / g <= 0 And we then return s', t' if either |s'| < |s| or |t'| < |t|. Without the fix, we could make the wrong choice in case |s'| = |s| = |b| / 2g Since it's a bit subtle, it would be nice with a review of this code before I commit it. Regards, /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 Sat Feb 17 17:17:21 2024 +0100 @@ -2960,12 +2960,13 @@ mpz_tdiv_q_2exp (t0, t0, 1); } - /* Arrange so that |s| < |u| / 2g */ + /* Arrange so that |s| < |u| / 2g and |t| < |v| / 2g, if possible. */ mpz_add (s1, s0, s1); - if (mpz_cmpabs (s0, s1) > 0) + mpz_sub (t1, t0, t1); + if (mpz_cmpabs (s0, s1) > 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.cWed Dec 27 19:35:36 2023 +0100 +++ b/mini-gmp/tests/t-gcd.cSat Feb 17 17:17:21 2024 +0100 @@ -79,20 +79,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 +102,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 +165,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 +192,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 Sat Feb 17 17:17:21 2024 +0100 @@ -440,8 +440,8 @@ if (mpz_sgn
Re: mini-gmp mpz_gcdext Bézout coefficients do not match documentation
Guido Vranken writes: > Computing extended gcd using mpz_gcdext where a = 1, b = 2: > > libgmp: g: 1, s: 1, t: 0 > mini-gmp: g: 1, s: -1, t: 1 > This violates the docs: "s = sgn(a) if abs(b) = 2", e.g. s must be 1 > > Computing extended gcd using mpz_gcdext where a = 6, b = 4: > libgmp: g: 2, s: 1, t: -1 > mini-libgmp: g: 2, s: -1, t: 2 I think mini-gmp is wrong here, I'll have to investigate. Thanks for reporting. > The docs state: "abs(s) < abs(b) / (2 g) and abs(t) < abs(a) / (2 g)" > I think that should be: "abs(s) <= abs(b) / (2 g) and abs(t) <= abs(a) / (2 > g)" > (< should be <=) These strict equalities don't hold for the documented "exceptional" cases. IT seems both your examples have abs(B) = 2G, so that's the case that applies. It is a bit subtle. And even if one relaxes the inequalities as you suggest, 1 = gcd(1,1) would still be an exception (since in that case, we'd get |s|, |t| < 1/2, i.e., s = t = 0, which isn't reasonable). Can you suggest any doc change to make it clearer? Regards, /Niels -- Niels Möller. PGP key CB4962D070D77D7FCB8BA36271D8F1FF368C6677. Internet email is subject to wholesale government surveillance. ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs
mini-gmp mpz_gcdext Bézout coefficients do not match documentation
Computing extended gcd using mpz_gcdext where a = 1, b = 2: libgmp: g: 1, s: 1, t: 0 mini-gmp: g: 1, s: -1, t: 1 This violates the docs: "s = sgn(a) if abs(b) = 2", e.g. s must be 1 Computing extended gcd using mpz_gcdext where a = 6, b = 4: libgmp: g: 2, s: 1, t: -1 mini-libgmp: g: 2, s: -1, t: 2 The docs state: "abs(s) < abs(b) / (2 g) and abs(t) < abs(a) / (2 g)" I think that should be: "abs(s) <= abs(b) / (2 g) and abs(t) <= abs(a) / (2 g)" (< should be <=) Then, for libgmp this is true: s: 1 <= 4 / 4, t: 1 <= 6 / 4 For mini-gmp this is false for t: s: 1 <= 4 / 4, t: 2 <= 6 / 4 Reproducer: //#include #include "mini-gmp.c" #include static void test( const unsigned long int a_ui, const unsigned long int b_ui) { mpz_t a, b, g, s, t; mpz_init(a); mpz_init(b); mpz_init(g); mpz_init(s); mpz_init(t); mpz_set_ui(a, a_ui); mpz_set_ui(b, b_ui); mpz_gcdext(g, s, t, a, b); printf("g: %s, s: %s, t: %s\n", mpz_get_str(NULL, 10, g), mpz_get_str(NULL, 10, s), mpz_get_str(NULL, 10, t)); } int main() { test(1, 2); test(6, 4); } Tested on Linux x64 using the latest repository code. ___ gmp-bugs mailing list gmp-bugs@gmplib.org https://gmplib.org/mailman/listinfo/gmp-bugs