Viewing file:  t-gcd.c (3.74 KB) -rw-r--r--
Select action/file-type:
 (+) |  (+) |  (+) | Code (+) | Session (+) |  (+) | SDB (+) |  (+) |  (+) |  (+) |  (+) |  (+) |
/*
Copyright 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 <limits.h>
#include <stdlib.h>
#include <stdio.h>
#include "testutils.h"
#define MAXBITS 400
#define COUNT 10000
/* Called when g is supposed to be gcd(a,b), and g = s a + t b. */
static int
gcdext_valid_p (const mpz_t a, const mpz_t b,
const mpz_t g, const mpz_t s, const mpz_t t)
{
mpz_t ta, tb, r;
/* 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, t = sgn (b)*/
return mpz_cmpabs (g, b) == 0
&& mpz_sgn (s) == 0 && mpz_cmp_si (t, mpz_sgn (b)) == 0;
}
else if (mpz_sgn (b) == 0)
{
/* Must have g == abs (a), s == sign (a), t = 0 */
return mpz_cmpabs (g, a) == 0
&& mpz_cmp_si (s, mpz_sgn (a)) == 0 && mpz_sgn (t) == 0;
}
if (mpz_sgn (g) <= 0)
return 0;
mpz_init (ta);
mpz_init (tb);
mpz_init (r);
mpz_mul (ta, s, a);
mpz_mul (tb, t, b);
mpz_add (ta, ta, tb);
if (mpz_cmp (ta, g) != 0)
{
fail:
mpz_clear (ta);
mpz_clear (tb);
mpz_clear (r);
return 0;
}
mpz_tdiv_qr (ta, r, a, g);
if (mpz_sgn (r) != 0)
goto fail;
mpz_tdiv_qr (tb, r, b, g);
if (mpz_sgn (r) != 0)
goto fail;
/* Require that 2 |s| < |b/g|, or |s| == 1. */
if (mpz_cmpabs_ui (s, 1) > 0)
{
mpz_mul_2exp (r, s, 1);
if (mpz_cmpabs (r, tb) > 0)
goto fail;
}
/* Require that 2 |t| < |a/g| or |t| == 1*/
if (mpz_cmpabs_ui (t, 1) > 0)
{
mpz_mul_2exp (r, t, 1);
if (mpz_cmpabs (r, ta) > 0)
return 0;
}
mpz_clear (ta);
mpz_clear (tb);
mpz_clear (r);
return 1;
}
void
testmain (int argc, char **argv)
{
unsigned i;
mpz_t a, b, g, s, t;
mpz_init (a);
mpz_init (b);
mpz_init (g);
mpz_init (s);
mpz_init (t);
for (i = 0; i < COUNT; i++)
{
mini_random_op3 (OP_GCD, MAXBITS, a, b, s);
mpz_gcd (g, a, b);
if (mpz_cmp (g, s))
{
fprintf (stderr, "mpz_gcd failed:\n");
dump ("a", a);
dump ("b", b);
dump ("r", g);
dump ("ref", s);
abort ();
}
}
for (i = 0; i < COUNT; i++)
{
unsigned flags;
mini_urandomb (a, 32);
flags = mpz_get_ui (a);
mini_rrandomb (a, MAXBITS);
mini_rrandomb (b, MAXBITS);
if (flags % 37 == 0)
mpz_mul (a, a, b);
if (flags % 37 == 1)
mpz_mul (b, a, b);
if (flags & 1)
mpz_neg (a, a);
if (flags & 2)
mpz_neg (b, b);
mpz_gcdext (g, s, t, a, b);
if (!gcdext_valid_p (a, b, g, s, t))
{
fprintf (stderr, "mpz_gcdext failed:\n");
dump ("a", a);
dump ("b", b);
dump ("g", g);
dump ("s", s);
dump ("t", t);
abort ();
}
mpz_gcd (s, a, b);
if (mpz_cmp (g, s))
{
fprintf (stderr, "mpz_gcd failed:\n");
dump ("a", a);
dump ("b", b);
dump ("r", g);
dump ("ref", s);
abort ();
}
}
mpz_clear (a);
mpz_clear (b);
mpz_clear (g);
mpz_clear (s);
mpz_clear (t);
}
|
