1 /*
2
3 Copyright 2012, Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library test suite.
6
7 The GNU MP Library test suite is free software; you can redistribute it
8 and/or modify it under the terms of the GNU General Public License as
9 published by the Free Software Foundation; either version 3 of the License,
10 or (at your option) any later version.
11
12 The GNU MP Library test suite is distributed in the hope that it will be
13 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
15 Public License for more details.
16
17 You should have received a copy of the GNU General Public License along with
18 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
19
20 #include <limits.h>
21 #include <stdlib.h>
22 #include <stdio.h>
23
24 #include "testutils.h"
25
26 #define MAXBITS 400
27 #define COUNT 10000
28
29 /* Called when g is supposed to be gcd(a,b), and g = s a + t b. */
30 static int
31 gcdext_valid_p (const mpz_t a, const mpz_t b,
32 const mpz_t g, const mpz_t s, const mpz_t t)
33 {
34 mpz_t ta, tb, r;
35
36 /* It's not clear that gcd(0,0) is well defined, but we allow it and
37 require that gcd(0,0) = 0. */
38 if (mpz_sgn (g) < 0)
39 return 0;
40
41 if (mpz_sgn (a) == 0)
42 {
43 /* Must have g == abs (b). Any value for s is in some sense "correct",
44 but it makes sense to require that s == 0, t = sgn (b)*/
45 return mpz_cmpabs (g, b) == 0
46 && mpz_sgn (s) == 0 && mpz_cmp_si (t, mpz_sgn (b)) == 0;
47 }
48 else if (mpz_sgn (b) == 0)
49 {
50 /* Must have g == abs (a), s == sign (a), t = 0 */
51 return mpz_cmpabs (g, a) == 0
52 && mpz_cmp_si (s, mpz_sgn (a)) == 0 && mpz_sgn (t) == 0;
53 }
54
55 if (mpz_sgn (g) <= 0)
56 return 0;
57
58 mpz_init (ta);
59 mpz_init (tb);
60 mpz_init (r);
61
62 mpz_mul (ta, s, a);
63 mpz_mul (tb, t, b);
64 mpz_add (ta, ta, tb);
65
66 if (mpz_cmp (ta, g) != 0)
67 {
68 fail:
69 mpz_clear (ta);
70 mpz_clear (tb);
71 mpz_clear (r);
72 return 0;
73 }
74 mpz_tdiv_qr (ta, r, a, g);
75 if (mpz_sgn (r) != 0)
76 goto fail;
77
78 mpz_tdiv_qr (tb, r, b, g);
79 if (mpz_sgn (r) != 0)
80 goto fail;
81
82 /* Require that 2 |s| < |b/g|, or |s| == 1. */
83 if (mpz_cmpabs_ui (s, 1) > 0)
84 {
85 mpz_mul_2exp (r, s, 1);
86 if (mpz_cmpabs (r, tb) > 0)
87 goto fail;
88 }
89
90 /* Require that 2 |t| < |a/g| or |t| == 1*/
91 if (mpz_cmpabs_ui (t, 1) > 0)
92 {
93 mpz_mul_2exp (r, t, 1);
94 if (mpz_cmpabs (r, ta) > 0)
95 return 0;
96 }
97
98 mpz_clear (ta);
99 mpz_clear (tb);
100 mpz_clear (r);
101
102 return 1;
103 }
104
105 void
106 testmain (int argc, char **argv)
107 {
108 unsigned i;
109 mpz_t a, b, g, s, t;
110
111 mpz_init (a);
112 mpz_init (b);
113 mpz_init (g);
114 mpz_init (s);
115 mpz_init (t);
116
117 for (i = 0; i < COUNT; i++)
118 {
119 mini_random_op3 (OP_GCD, MAXBITS, a, b, s);
120 mpz_gcd (g, a, b);
121 if (mpz_cmp (g, s))
122 {
123 fprintf (stderr, "mpz_gcd failed:\n");
124 dump ("a", a);
125 dump ("b", b);
126 dump ("r", g);
127 dump ("ref", s);
128 abort ();
129 }
130 }
131
132 for (i = 0; i < COUNT; i++)
133 {
134 unsigned flags;
135 mini_urandomb (a, 32);
136 flags = mpz_get_ui (a);
137 mini_rrandomb (a, MAXBITS);
138 mini_rrandomb (b, MAXBITS);
139
140 if (flags % 37 == 0)
141 mpz_mul (a, a, b);
142 if (flags % 37 == 1)
143 mpz_mul (b, a, b);
144
145 if (flags & 1)
146 mpz_neg (a, a);
147 if (flags & 2)
148 mpz_neg (b, b);
149
150 mpz_gcdext (g, s, t, a, b);
151 if (!gcdext_valid_p (a, b, g, s, t))
152 {
153 fprintf (stderr, "mpz_gcdext failed:\n");
154 dump ("a", a);
155 dump ("b", b);
156 dump ("g", g);
157 dump ("s", s);
158 dump ("t", t);
159 abort ();
160 }
161
162 mpz_gcd (s, a, b);
163 if (mpz_cmp (g, s))
164 {
165 fprintf (stderr, "mpz_gcd failed:\n");
166 dump ("a", a);
167 dump ("b", b);
168 dump ("r", g);
169 dump ("ref", s);
170 abort ();
171 }
172 }
173 mpz_clear (a);
174 mpz_clear (b);
175 mpz_clear (g);
176 mpz_clear (s);
177 mpz_clear (t);
178 }