tests/mpz/t-gcd.c - gmp - Git at Google

 /* 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 <stdio.h>
 #include <stdlib.h>

 #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;
 }
	/* 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 <stdio.h>
	#include <stdlib.h>

	#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;
	}