/* 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.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_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);

    }

  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;

}