tests/mpz/t-gcd_ui.c (check_ui_factors): New check triggering a possible bug.

1 files changed, 94 insertions, 0 deletions

diff --git a/tests/mpz/t-gcd_ui.c b/tests/mpz/t-gcd_ui.c
index 03b63941c..3f56a97e6 100644
--- a/tests/mpz/t-gcd_ui.c
+++ b/tests/mpz/t-gcd_ui.c
@@ -50,12 +50,106 @@ check_ui_range (void)
   mpz_clear (x);
 }
 
+static void
+check_ui_factors (void)
+{
+#define NUM_FACTORS 9
+  static const char* factors[NUM_FACTORS] = {
+    "641", "274177", "3", "5", "17", "257", "65537",
+    "59649589127497217", "1238926361552897" };
+  unsigned long  got;
+  mpz_t  x, b, d, f, g;
+  int  i, j;
+  gmp_randstate_ptr rands;
+
+  if (GMP_NUMB_BITS < 5 || GMP_NUMB_BITS == 8
+      || GMP_NUMB_BITS == 16 || GMP_NUMB_BITS > 511)
+    {
+      printf ("No usable factors for 2^%i+1.\n", GMP_NUMB_BITS);
+      return;
+    }
+
+  mpz_init (x);
+  mpz_init (d);
+  mpz_init (f);
+  mpz_init (g);
+
+  mpz_setbit (x, GMP_NUMB_BITS);
+  mpz_add_ui (x, x, 1);
+
+  for (i = 0; i < NUM_FACTORS; ++i)
+    {
+      mpz_set_str (f, factors[i], 10);
+      if (mpz_divisible_p (x, f))
+	{
+	  mpz_mul_2exp (f, f, 1);
+	  /* d is an odd multiple of the factor f, exactly filling a limb. */
+	  mpz_sub (d, x, f);
+	  /* f = 2^GMP_NUMB_BITS mod d. */
+	  mpz_sub_ui (f, f, 1);
+	  break;
+	}
+    }
+
+  mpz_gcd (g, f, d);
+  if (mpz_even_p (d) || mpz_cmp (d, f) <= 0 || mpz_cmp_ui (g, 1) != 0)
+    {
+      printf ("No usable factor found.\n");
+      abort ();
+    }
+
+  rands = RANDS;
+  mpz_mul_ui (x, d, gmp_urandomm_ui (rands, 30000) + 1);
+
+  mpz_init (b);
+  mpz_setbit (b, GMP_NUMB_BITS - 1);
+  for (j = 0; j < 4; ++j)
+    {
+      mpz_add (x, x, b);
+
+      for (i = 1; i >= -1; --i)
+	{
+	  if (mpz_fits_ulong_p (d)
+	      && ((got = mpz_gcd_ui (NULL, x, mpz_get_ui (d)))
+		  != (i != 0 ? 1 : mpz_get_ui (d))))
+	    {
+	      printf ("mpz_gcd_ui (f, kV+%i*2^%i, V): error (j = %i)\n", i, GMP_NUMB_BITS - 1, j);
+	      printf ("   return %#lx\n", got);
+	      printf ("   should be %#lx\n", (i != 0 ? 1 : mpz_get_ui (d)));
+	      abort ();
+	    }
+
+	  mpz_gcd (g, x, d);
+	  if ((mpz_cmp_ui (g, 1) == 0) != (i != 0))
+	    {
+	      printf ("mpz_gcd (f, kV+%i*2^%i, V): error (j = %i)\n", i, GMP_NUMB_BITS - 1, j);
+	      printf ("   should%s be one.\n",(i != 0 ? "" : " not"));
+	      abort ();
+	    }
+
+	  mpz_sub (x, x, b);
+	}
+      /* Back to the original x. */
+      mpz_addmul_ui (x, b, 2);
+      mpz_mul (b, b, f);
+      mpz_mod (b, b, d);
+    }
+
+  mpz_clear (g);
+  mpz_clear (x);
+  mpz_clear (f);
+  mpz_clear (d);
+  mpz_clear (b);
+}
+
+
 int
 main (void)
 {
   tests_start ();
 
   check_ui_range ();
+  check_ui_factors ();
 
   tests_end ();
   exit (0);