/* 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;
}