diff options
Diffstat (limited to 'vendor/gmp-6.3.0/tests/mpn/t-hgcd.c')
-rw-r--r-- | vendor/gmp-6.3.0/tests/mpn/t-hgcd.c | 406 |
1 files changed, 406 insertions, 0 deletions
diff --git a/vendor/gmp-6.3.0/tests/mpn/t-hgcd.c b/vendor/gmp-6.3.0/tests/mpn/t-hgcd.c new file mode 100644 index 0000000..07f83e9 --- /dev/null +++ b/vendor/gmp-6.3.0/tests/mpn/t-hgcd.c @@ -0,0 +1,406 @@ +/* Test mpn_hgcd. + +Copyright 1991, 1993, 1994, 1996, 1997, 2000-2004 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 <stdio.h> +#include <stdlib.h> + +#include "gmp-impl.h" +#include "tests.h" + +static mp_size_t one_test (mpz_t, mpz_t, int); +static void debug_mp (mpz_t, int); + +#define MIN_OPERAND_SIZE 2 + +/* Fixed values, for regression testing of mpn_hgcd. */ +struct value { int res; const char *a; const char *b; }; +static const struct value hgcd_values[] = { +#if GMP_NUMB_BITS == 32 + { 5, + "0x1bddff867272a9296ac493c251d7f46f09a5591fe", + "0xb55930a2a68a916450a7de006031068c5ddb0e5c" }, + { 4, + "0x2f0ece5b1ee9c15e132a01d55768dc13", + "0x1c6f4fd9873cdb24466e6d03e1cc66e7" }, + { 3, "0x7FFFFC003FFFFFFFFFC5", "0x3FFFFE001FFFFFFFFFE3"}, +#endif + { -1, NULL, NULL } +}; + +struct hgcd_ref +{ + mpz_t m[2][2]; +}; + +static void hgcd_ref_init (struct hgcd_ref *); +static void hgcd_ref_clear (struct hgcd_ref *); +static int hgcd_ref (struct hgcd_ref *, mpz_t, mpz_t); +static int hgcd_ref_equal (const struct hgcd_matrix *, const struct hgcd_ref *); + +int +main (int argc, char **argv) +{ + mpz_t op1, op2, temp1, temp2; + int i, j, chain_len; + gmp_randstate_ptr rands; + mpz_t bs; + unsigned long size_range; + + tests_start (); + rands = RANDS; + + mpz_init (bs); + mpz_init (op1); + mpz_init (op2); + mpz_init (temp1); + mpz_init (temp2); + + for (i = 0; hgcd_values[i].res >= 0; i++) + { + mp_size_t res; + + mpz_set_str (op1, hgcd_values[i].a, 0); + mpz_set_str (op2, hgcd_values[i].b, 0); + + res = one_test (op1, op2, -1-i); + if (res != hgcd_values[i].res) + { + fprintf (stderr, "ERROR in test %d\n", -1-i); + fprintf (stderr, "Bad return code from hgcd\n"); + fprintf (stderr, "op1="); debug_mp (op1, -16); + fprintf (stderr, "op2="); debug_mp (op2, -16); + fprintf (stderr, "expected: %d\n", hgcd_values[i].res); + fprintf (stderr, "hgcd: %d\n", (int) res); + abort (); + } + } + + for (i = 0; i < 15; i++) + { + /* Generate plain operands with unknown gcd. These types of operands + have proven to trigger certain bugs in development versions of the + gcd code. */ + + mpz_urandomb (bs, rands, 32); + size_range = mpz_get_ui (bs) % 13 + 2; + + mpz_urandomb (bs, rands, size_range); + mpz_rrandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE); + mpz_urandomb (bs, rands, size_range); + mpz_rrandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE); + + if (mpz_cmp (op1, op2) < 0) + mpz_swap (op1, op2); + + if (mpz_size (op1) > 0) + one_test (op1, op2, i); + + /* Generate a division chain backwards, allowing otherwise + unlikely huge quotients. */ + + mpz_set_ui (op1, 0); + mpz_urandomb (bs, rands, 32); + mpz_urandomb (bs, rands, mpz_get_ui (bs) % 16 + 1); + mpz_rrandomb (op2, rands, mpz_get_ui (bs)); + mpz_add_ui (op2, op2, 1); + +#if 0 + chain_len = 1000000; +#else + mpz_urandomb (bs, rands, 32); + chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * GCD_DC_THRESHOLD / 256); +#endif + + for (j = 0; j < chain_len; j++) + { + mpz_urandomb (bs, rands, 32); + mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1); + mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1); + mpz_add_ui (temp2, temp2, 1); + mpz_mul (temp1, op2, temp2); + mpz_add (op1, op1, temp1); + + /* Don't generate overly huge operands. */ + if (SIZ (op1) > 3 * GCD_DC_THRESHOLD) + break; + + mpz_urandomb (bs, rands, 32); + mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1); + mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1); + mpz_add_ui (temp2, temp2, 1); + mpz_mul (temp1, op1, temp2); + mpz_add (op2, op2, temp1); + + /* Don't generate overly huge operands. */ + if (SIZ (op2) > 3 * GCD_DC_THRESHOLD) + break; + } + if (mpz_cmp (op1, op2) < 0) + mpz_swap (op1, op2); + + if (mpz_size (op1) > 0) + one_test (op1, op2, i); + } + + mpz_clear (bs); + mpz_clear (op1); + mpz_clear (op2); + mpz_clear (temp1); + mpz_clear (temp2); + + tests_end (); + exit (0); +} + +static void +debug_mp (mpz_t x, int base) +{ + mpz_out_str (stderr, base, x); fputc ('\n', stderr); +} + +static int +mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize); + +static mp_size_t +one_test (mpz_t a, mpz_t b, int i) +{ + struct hgcd_matrix hgcd; + struct hgcd_ref ref; + + mpz_t ref_r0; + mpz_t ref_r1; + mpz_t hgcd_r0; + mpz_t hgcd_r1; + + mp_size_t res[2]; + mp_size_t asize; + mp_size_t bsize; + + mp_size_t hgcd_init_scratch; + mp_size_t hgcd_scratch; + + mp_ptr hgcd_init_tp; + mp_ptr hgcd_tp; + + asize = a->_mp_size; + bsize = b->_mp_size; + + ASSERT (asize >= bsize); + + hgcd_init_scratch = MPN_HGCD_MATRIX_INIT_ITCH (asize); + hgcd_init_tp = refmpn_malloc_limbs (hgcd_init_scratch); + mpn_hgcd_matrix_init (&hgcd, asize, hgcd_init_tp); + + hgcd_scratch = mpn_hgcd_itch (asize); + hgcd_tp = refmpn_malloc_limbs (hgcd_scratch); + +#if 0 + fprintf (stderr, + "one_test: i = %d asize = %d, bsize = %d\n", + i, a->_mp_size, b->_mp_size); + + gmp_fprintf (stderr, + "one_test: i = %d\n" + " a = %Zx\n" + " b = %Zx\n", + i, a, b); +#endif + hgcd_ref_init (&ref); + + mpz_init_set (ref_r0, a); + mpz_init_set (ref_r1, b); + res[0] = hgcd_ref (&ref, ref_r0, ref_r1); + + mpz_init_set (hgcd_r0, a); + mpz_init_set (hgcd_r1, b); + if (bsize < asize) + { + _mpz_realloc (hgcd_r1, asize); + MPN_ZERO (hgcd_r1->_mp_d + bsize, asize - bsize); + } + res[1] = mpn_hgcd (hgcd_r0->_mp_d, + hgcd_r1->_mp_d, + asize, + &hgcd, hgcd_tp); + + if (res[0] != res[1]) + { + fprintf (stderr, "ERROR in test %d\n", i); + fprintf (stderr, "Different return value from hgcd and hgcd_ref\n"); + fprintf (stderr, "op1="); debug_mp (a, -16); + fprintf (stderr, "op2="); debug_mp (b, -16); + fprintf (stderr, "hgcd_ref: %ld\n", (long) res[0]); + fprintf (stderr, "mpn_hgcd: %ld\n", (long) res[1]); + abort (); + } + if (res[0] > 0) + { + if (!hgcd_ref_equal (&hgcd, &ref) + || !mpz_mpn_equal (ref_r0, hgcd_r0->_mp_d, res[1]) + || !mpz_mpn_equal (ref_r1, hgcd_r1->_mp_d, res[1])) + { + fprintf (stderr, "ERROR in test %d\n", i); + fprintf (stderr, "mpn_hgcd and hgcd_ref returned different values\n"); + fprintf (stderr, "op1="); debug_mp (a, -16); + fprintf (stderr, "op2="); debug_mp (b, -16); + abort (); + } + } + + refmpn_free_limbs (hgcd_init_tp); + refmpn_free_limbs (hgcd_tp); + hgcd_ref_clear (&ref); + mpz_clear (ref_r0); + mpz_clear (ref_r1); + mpz_clear (hgcd_r0); + mpz_clear (hgcd_r1); + + return res[0]; +} + +static void +hgcd_ref_init (struct hgcd_ref *hgcd) +{ + unsigned i; + for (i = 0; i<2; i++) + { + unsigned j; + for (j = 0; j<2; j++) + mpz_init (hgcd->m[i][j]); + } +} + +static void +hgcd_ref_clear (struct hgcd_ref *hgcd) +{ + unsigned i; + for (i = 0; i<2; i++) + { + unsigned j; + for (j = 0; j<2; j++) + mpz_clear (hgcd->m[i][j]); + } +} + + +static int +sdiv_qr (mpz_t q, mpz_t r, mp_size_t s, const mpz_t a, const mpz_t b) +{ + mpz_fdiv_qr (q, r, a, b); + if (mpz_size (r) <= s) + { + mpz_add (r, r, b); + mpz_sub_ui (q, q, 1); + } + + return (mpz_sgn (q) > 0); +} + +static int +hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b) +{ + mp_size_t n = MAX (mpz_size (a), mpz_size (b)); + mp_size_t s = n/2 + 1; + mp_size_t asize; + mp_size_t bsize; + mpz_t q; + int res; + + if (mpz_size (a) <= s || mpz_size (b) <= s) + return 0; + + res = mpz_cmp (a, b); + if (res < 0) + { + mpz_sub (b, b, a); + if (mpz_size (b) <= s) + return 0; + + mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 0); + mpz_set_ui (hgcd->m[1][0], 1); mpz_set_ui (hgcd->m[1][1], 1); + } + else if (res > 0) + { + mpz_sub (a, a, b); + if (mpz_size (a) <= s) + return 0; + + mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 1); + mpz_set_ui (hgcd->m[1][0], 0); mpz_set_ui (hgcd->m[1][1], 1); + } + else + return 0; + + mpz_init (q); + + for (;;) + { + ASSERT (mpz_size (a) > s); + ASSERT (mpz_size (b) > s); + + if (mpz_cmp (a, b) > 0) + { + if (!sdiv_qr (q, a, s, a, b)) + break; + mpz_addmul (hgcd->m[0][1], q, hgcd->m[0][0]); + mpz_addmul (hgcd->m[1][1], q, hgcd->m[1][0]); + } + else + { + if (!sdiv_qr (q, b, s, b, a)) + break; + mpz_addmul (hgcd->m[0][0], q, hgcd->m[0][1]); + mpz_addmul (hgcd->m[1][0], q, hgcd->m[1][1]); + } + } + + mpz_clear (q); + + asize = mpz_size (a); + bsize = mpz_size (b); + return MAX (asize, bsize); +} + +static int +mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize) +{ + mp_srcptr ap = a->_mp_d; + mp_size_t asize = a->_mp_size; + + MPN_NORMALIZE (bp, bsize); + return asize == bsize && mpn_cmp (ap, bp, asize) == 0; +} + +static int +hgcd_ref_equal (const struct hgcd_matrix *hgcd, const struct hgcd_ref *ref) +{ + unsigned i; + + for (i = 0; i<2; i++) + { + unsigned j; + + for (j = 0; j<2; j++) + if (!mpz_mpn_equal (ref->m[i][j], hgcd->p[i][j], hgcd->n)) + return 0; + } + + return 1; +} |