From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001 From: Duncan Wilkie Date: Sat, 18 Nov 2023 06:11:09 -0600 Subject: Initial commit. --- gmp-6.3.0/tests/mpz/t-gcd.c | 467 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 467 insertions(+) create mode 100644 gmp-6.3.0/tests/mpz/t-gcd.c (limited to 'gmp-6.3.0/tests/mpz/t-gcd.c') diff --git a/gmp-6.3.0/tests/mpz/t-gcd.c b/gmp-6.3.0/tests/mpz/t-gcd.c new file mode 100644 index 0000000..46b278a --- /dev/null +++ b/gmp-6.3.0/tests/mpz/t-gcd.c @@ -0,0 +1,467 @@ +/* Test mpz_gcd, mpz_gcdext, and mpz_gcd_ui. + +Copyright 1991, 1993, 1994, 1996, 1997, 2000-2005, 2008, 2009, 2012 Free +Software Foundation, Inc. + +This file is part of the GNU MP Library test suite. + +The GNU MP Library test suite is free software; you can redistribute it +and/or modify it under the terms of the GNU General Public License as +published by the Free Software Foundation; either version 3 of the License, +or (at your option) any later version. + +The GNU MP Library test suite is distributed in the hope that it will be +useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +Public License for more details. + +You should have received a copy of the GNU General Public License along with +the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */ + +#include +#include + +#include "gmp-impl.h" +#include "tests.h" + +void one_test (mpz_t, mpz_t, mpz_t, int); +void debug_mp (mpz_t, int); + +static int gcdext_valid_p (const mpz_t, const mpz_t, const mpz_t, const mpz_t); + +/* Keep one_test's variables global, so that we don't need + to reinitialize them for each test. */ +mpz_t gcd1, gcd2, s, temp1, temp2, temp3; + +#define MAX_SCHOENHAGE_THRESHOLD HGCD_REDUCE_THRESHOLD + +/* Define this to make all operands be large enough for Schoenhage gcd + to be used. */ +#ifndef WHACK_SCHOENHAGE +#define WHACK_SCHOENHAGE 0 +#endif + +#if WHACK_SCHOENHAGE +#define MIN_OPERAND_BITSIZE (MAX_SCHOENHAGE_THRESHOLD * GMP_NUMB_BITS) +#else +#define MIN_OPERAND_BITSIZE 1 +#endif + + +void +check_data (void) +{ + static const struct { + const char *a; + const char *b; + const char *want; + } data[] = { + /* This tickled a bug in gmp 4.1.2 mpn/x86/k6/gcd_finda.asm. */ + { "0x3FFC000007FFFFFFFFFF00000000003F83FFFFFFFFFFFFFFF80000000000000001", + "0x1FFE0007FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFC000000000000000000000001", + "5" } + }; + + mpz_t a, b, got, want; + int i; + + mpz_inits (a, b, got, want, NULL); + + for (i = 0; i < numberof (data); i++) + { + mpz_set_str_or_abort (a, data[i].a, 0); + mpz_set_str_or_abort (b, data[i].b, 0); + mpz_set_str_or_abort (want, data[i].want, 0); + mpz_gcd (got, a, b); + MPZ_CHECK_FORMAT (got); + if (mpz_cmp (got, want) != 0) + { + printf ("mpz_gcd wrong on data[%d]\n", i); + printf (" a %s\n", data[i].a); + printf (" b %s\n", data[i].b); + mpz_trace (" a", a); + mpz_trace (" b", b); + mpz_trace (" want", want); + mpz_trace (" got ", got); + abort (); + } + } + + mpz_clears (a, b, got, want, NULL); +} + +void +make_chain_operands (mpz_t ref, mpz_t a, mpz_t b, gmp_randstate_t rs, int nb1, int nb2, int chain_len) +{ + mpz_t bs, temp1, temp2; + int j; + + mpz_inits (bs, temp1, temp2, NULL); + + /* Generate a division chain backwards, allowing otherwise unlikely huge + quotients. */ + + mpz_set_ui (a, 0); + mpz_urandomb (bs, rs, 32); + mpz_urandomb (bs, rs, mpz_get_ui (bs) % nb1 + 1); + mpz_rrandomb (b, rs, mpz_get_ui (bs)); + mpz_add_ui (b, b, 1); + mpz_set (ref, b); + + for (j = 0; j < chain_len; j++) + { + mpz_urandomb (bs, rs, 32); + mpz_urandomb (bs, rs, mpz_get_ui (bs) % nb2 + 1); + mpz_rrandomb (temp2, rs, mpz_get_ui (bs) + 1); + mpz_add_ui (temp2, temp2, 1); + mpz_mul (temp1, b, temp2); + mpz_add (a, a, temp1); + + mpz_urandomb (bs, rs, 32); + mpz_urandomb (bs, rs, mpz_get_ui (bs) % nb2 + 1); + mpz_rrandomb (temp2, rs, mpz_get_ui (bs) + 1); + mpz_add_ui (temp2, temp2, 1); + mpz_mul (temp1, a, temp2); + mpz_add (b, b, temp1); + } + + mpz_clears (bs, temp1, temp2, NULL); +} + +/* Test operands from a table of seed data. This variant creates the operands + using plain ol' mpz_rrandomb. This is a hack for better coverage of the gcd + code, which depends on that the random number generators give the exact + numbers we expect. */ +void +check_kolmo1 (void) +{ + static const struct { + unsigned int seed; + int nb; + const char *want; + } data[] = { + { 59618, 38208, "5"}, + { 76521, 49024, "3"}, + { 85869, 54976, "1"}, + { 99449, 63680, "1"}, + {112453, 72000, "1"} + }; + + gmp_randstate_t rs; + mpz_t bs, a, b, want; + int i, unb, vnb, nb; + + gmp_randinit_default (rs); + + mpz_inits (bs, a, b, want, NULL); + + for (i = 0; i < numberof (data); i++) + { + nb = data[i].nb; + + gmp_randseed_ui (rs, data[i].seed); + + mpz_urandomb (bs, rs, 32); + unb = mpz_get_ui (bs) % nb; + mpz_urandomb (bs, rs, 32); + vnb = mpz_get_ui (bs) % nb; + + mpz_rrandomb (a, rs, unb); + mpz_rrandomb (b, rs, vnb); + + mpz_set_str_or_abort (want, data[i].want, 0); + + one_test (a, b, want, -1); + } + + mpz_clears (bs, a, b, want, NULL); + gmp_randclear (rs); +} + +/* Test operands from a table of seed data. This variant creates the operands + using a division chain. This is a hack for better coverage of the gcd + code, which depends on that the random number generators give the exact + numbers we expect. */ +void +check_kolmo2 (void) +{ + static const struct { + unsigned int seed; + int nb, chain_len; + } data[] = { + { 917, 15, 5 }, + { 1032, 18, 6 }, + { 1167, 18, 6 }, + { 1174, 18, 6 }, + { 1192, 18, 6 }, + }; + + gmp_randstate_t rs; + mpz_t bs, a, b, want; + int i; + + gmp_randinit_default (rs); + + mpz_inits (bs, a, b, want, NULL); + + for (i = 0; i < numberof (data); i++) + { + gmp_randseed_ui (rs, data[i].seed); + make_chain_operands (want, a, b, rs, data[i].nb, data[i].nb, data[i].chain_len); + one_test (a, b, want, -1); + } + + mpz_clears (bs, a, b, want, NULL); + gmp_randclear (rs); +} + +int +main (int argc, char **argv) +{ + mpz_t op1, op2, ref; + int i, chain_len; + gmp_randstate_ptr rands; + mpz_t bs; + unsigned long bsi, size_range; + long int reps = 200; + + tests_start (); + TESTS_REPS (reps, argv, argc); + + rands = RANDS; + + mpz_inits (bs, op1, op2, ref, gcd1, gcd2, temp1, temp2, temp3, s, NULL); + + check_data (); + check_kolmo1 (); + check_kolmo2 (); + + /* Testcase to exercise the u0 == u1 case in mpn_gcdext_lehmer_n. */ + /* mpz_set_ui (op2, GMP_NUMB_MAX); */ /* FIXME: Huge limb doesn't always fit */ + mpz_set_ui (op2, 0); + mpz_setbit (op2, GMP_NUMB_BITS); + mpz_sub_ui (op2, op2, 1); + mpz_mul_2exp (op1, op2, 100); + mpz_add (op1, op1, op2); + mpz_mul_ui (op2, op2, 2); + one_test (op1, op2, NULL, -1); + + for (i = 0; i < reps; i++) + { + /* Generate plain operands with unknown gcd. These types of operands + have proven to trigger certain bugs in development versions of the + gcd code. The "hgcd->row[3].rsize > M" ASSERT is not triggered by + the division chain code below, but that is most likely just a result + of that other ASSERTs are triggered before it. */ + + mpz_urandomb (bs, rands, 32); + size_range = mpz_get_ui (bs) % 17 + 2; + + mpz_urandomb (bs, rands, size_range); + mpz_rrandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE); + mpz_urandomb (bs, rands, size_range); + mpz_rrandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE); + + mpz_urandomb (bs, rands, 8); + bsi = mpz_get_ui (bs); + + if ((bsi & 0x3c) == 4) + mpz_mul (op1, op1, op2); /* make op1 a multiple of op2 */ + else if ((bsi & 0x3c) == 8) + mpz_mul (op2, op1, op2); /* make op2 a multiple of op1 */ + + if ((bsi & 1) != 0) + mpz_neg (op1, op1); + if ((bsi & 2) != 0) + mpz_neg (op2, op2); + + one_test (op1, op2, NULL, i); + + /* Generate a division chain backwards, allowing otherwise unlikely huge + quotients. */ + + mpz_urandomb (bs, rands, 32); + chain_len = mpz_get_ui (bs) % LOG2C (GMP_NUMB_BITS * MAX_SCHOENHAGE_THRESHOLD); + mpz_urandomb (bs, rands, 32); + chain_len = mpz_get_ui (bs) % (1 << chain_len) / 32; + + make_chain_operands (ref, op1, op2, rands, 16, 12, chain_len); + + one_test (op1, op2, ref, i); + } + + /* Check that we can use NULL as first argument of mpz_gcdext. */ + mpz_set_si (op1, -10); + mpz_set_si (op2, 0); + mpz_gcdext (NULL, temp1, temp2, op1, op2); + ASSERT_ALWAYS (mpz_cmp_si (temp1, -1) == 0); + ASSERT_ALWAYS (mpz_cmp_si (temp2, 0) == 0); + mpz_set_si (op2, 6); + mpz_gcdext (NULL, temp1, temp2, op1, op2); + ASSERT_ALWAYS (mpz_cmp_si (temp1, 1) == 0); + ASSERT_ALWAYS (mpz_cmp_si (temp2, 2) == 0); + + mpz_clears (bs, op1, op2, ref, gcd1, gcd2, temp1, temp2, temp3, s, NULL); + + tests_end (); + exit (0); +} + +void +debug_mp (mpz_t x, int base) +{ + mpz_out_str (stderr, base, x); fputc ('\n', stderr); +} + +void +one_test (mpz_t op1, mpz_t op2, mpz_t ref, int i) +{ + /* + printf ("%d %d %d\n", SIZ (op1), SIZ (op2), ref != NULL ? SIZ (ref) : 0); + fflush (stdout); + */ + + /* + fprintf (stderr, "op1="); debug_mp (op1, -16); + fprintf (stderr, "op2="); debug_mp (op2, -16); + */ + + mpz_gcdext (gcd1, s, NULL, op1, op2); + MPZ_CHECK_FORMAT (gcd1); + MPZ_CHECK_FORMAT (s); + + if (ref && mpz_cmp (ref, gcd1) != 0) + { + fprintf (stderr, "ERROR in test %d\n", i); + fprintf (stderr, "mpz_gcdext returned incorrect result\n"); + fprintf (stderr, "op1="); debug_mp (op1, -16); + fprintf (stderr, "op2="); debug_mp (op2, -16); + fprintf (stderr, "expected result:\n"); debug_mp (ref, -16); + fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd1, -16); + abort (); + } + + if (!gcdext_valid_p(op1, op2, gcd1, s)) + { + fprintf (stderr, "ERROR in test %d\n", i); + fprintf (stderr, "mpz_gcdext returned invalid result\n"); + fprintf (stderr, "op1="); debug_mp (op1, -16); + fprintf (stderr, "op2="); debug_mp (op2, -16); + fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd1, -16); + fprintf (stderr, "s="); debug_mp (s, -16); + abort (); + } + + mpz_gcd (gcd2, op1, op2); + MPZ_CHECK_FORMAT (gcd2); + + if (mpz_cmp (gcd2, gcd1) != 0) + { + fprintf (stderr, "ERROR in test %d\n", i); + fprintf (stderr, "mpz_gcd returned incorrect result\n"); + fprintf (stderr, "op1="); debug_mp (op1, -16); + fprintf (stderr, "op2="); debug_mp (op2, -16); + fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16); + fprintf (stderr, "mpz_gcd returns:\n"); debug_mp (gcd2, -16); + abort (); + } + + /* This should probably move to t-gcd_ui.c */ + if (mpz_fits_ulong_p (op1) || mpz_fits_ulong_p (op2)) + { + if (mpz_fits_ulong_p (op1)) + mpz_gcd_ui (gcd2, op2, mpz_get_ui (op1)); + else + mpz_gcd_ui (gcd2, op1, mpz_get_ui (op2)); + if (mpz_cmp (gcd2, gcd1)) + { + fprintf (stderr, "ERROR in test %d\n", i); + fprintf (stderr, "mpz_gcd_ui returned incorrect result\n"); + fprintf (stderr, "op1="); debug_mp (op1, -16); + fprintf (stderr, "op2="); debug_mp (op2, -16); + fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16); + fprintf (stderr, "mpz_gcd_ui returns:\n"); debug_mp (gcd2, -16); + abort (); + } + } + + mpz_gcdext (gcd2, temp1, temp2, op1, op2); + MPZ_CHECK_FORMAT (gcd2); + MPZ_CHECK_FORMAT (temp1); + MPZ_CHECK_FORMAT (temp2); + + mpz_mul (temp1, temp1, op1); + mpz_mul (temp2, temp2, op2); + mpz_add (temp1, temp1, temp2); + + if (mpz_cmp (gcd1, gcd2) != 0 + || mpz_cmp (gcd2, temp1) != 0) + { + fprintf (stderr, "ERROR in test %d\n", i); + fprintf (stderr, "mpz_gcdext returned incorrect result\n"); + fprintf (stderr, "op1="); debug_mp (op1, -16); + fprintf (stderr, "op2="); debug_mp (op2, -16); + fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16); + fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd2, -16); + abort (); + } +} + +/* Called when g is supposed to be gcd(a,b), and g = s a + t b, for some t. + Uses temp1, temp2 and temp3. */ +static int +gcdext_valid_p (const mpz_t a, const mpz_t b, const mpz_t g, const mpz_t s) +{ + /* It's not clear that gcd(0,0) is well defined, but we allow it and require that + gcd(0,0) = 0. */ + if (mpz_sgn (g) < 0) + return 0; + + if (mpz_sgn (a) == 0) + { + /* Must have g == abs (b). Any value for s is in some sense "correct", + but it makes sense to require that s == 0. */ + return mpz_cmpabs (g, b) == 0 && mpz_sgn (s) == 0; + } + else if (mpz_sgn (b) == 0) + { + /* Must have g == abs (a), s == sign (a) */ + return mpz_cmpabs (g, a) == 0 && mpz_cmp_si (s, mpz_sgn (a)) == 0; + } + + if (mpz_sgn (g) <= 0) + return 0; + + mpz_tdiv_qr (temp1, temp3, a, g); + if (mpz_sgn (temp3) != 0) + return 0; + + mpz_tdiv_qr (temp2, temp3, b, g); + if (mpz_sgn (temp3) != 0) + return 0; + + /* Require that 2 |s| < |b/g|, or |s| == 1. */ + if (mpz_cmpabs_ui (s, 1) > 0) + { + mpz_mul_2exp (temp3, s, 1); + if (mpz_cmpabs (temp3, temp2) >= 0) + return 0; + } + + /* Compute the other cofactor. */ + mpz_mul(temp2, s, a); + mpz_sub(temp2, g, temp2); + mpz_tdiv_qr(temp2, temp3, temp2, b); + + if (mpz_sgn (temp3) != 0) + return 0; + + /* Require that 2 |t| < |a/g| or |t| == 1*/ + if (mpz_cmpabs_ui (temp2, 1) > 0) + { + mpz_mul_2exp (temp2, temp2, 1); + if (mpz_cmpabs (temp2, temp1) >= 0) + return 0; + } + return 1; +} -- cgit v1.2.3