aboutsummaryrefslogtreecommitdiff
path: root/gmp-6.3.0/tests/mpn/t-hgcd.c
diff options
context:
space:
mode:
Diffstat (limited to 'gmp-6.3.0/tests/mpn/t-hgcd.c')
-rw-r--r--gmp-6.3.0/tests/mpn/t-hgcd.c406
1 files changed, 406 insertions, 0 deletions
diff --git a/gmp-6.3.0/tests/mpn/t-hgcd.c b/gmp-6.3.0/tests/mpn/t-hgcd.c
new file mode 100644
index 0000000..07f83e9
--- /dev/null
+++ b/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;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-09-19 16:54:58 +0000