/* Test mpn_hgcd.

Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 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 http://www.gnu.org/licenses/.  */

#include <stdio.h>
#include <stdlib.h>

#include "gmp.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;
}