1 /* Test mpz_gcd, mpz_gcdext, and mpz_gcd_ui.
2
3 Copyright 1991, 1993, 1994, 1996, 1997, 2000-2005, 2008, 2009, 2012 Free
4 Software Foundation, Inc.
5
6 This file is part of the GNU MP Library test suite.
7
8 The GNU MP Library test suite is free software; you can redistribute it
9 and/or modify it under the terms of the GNU General Public License as
10 published by the Free Software Foundation; either version 3 of the License,
11 or (at your option) any later version.
12
13 The GNU MP Library test suite is distributed in the hope that it will be
14 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
16 Public License for more details.
17
18 You should have received a copy of the GNU General Public License along with
19 the GNU MP Library test suite. If not, see https://www.gnu.org/licenses/. */
20
21 #include <stdio.h>
22 #include <stdlib.h>
23
24 #include "gmp-impl.h"
25 #include "tests.h"
26
27 void one_test (mpz_t, mpz_t, mpz_t, int);
28 void debug_mp (mpz_t, int);
29
30 static int gcdext_valid_p (const mpz_t, const mpz_t, const mpz_t, const mpz_t);
31
32 /* Keep one_test's variables global, so that we don't need
33 to reinitialize them for each test. */
34 mpz_t gcd1, gcd2, s, temp1, temp2, temp3;
35
36 #define MAX_SCHOENHAGE_THRESHOLD HGCD_REDUCE_THRESHOLD
37
38 /* Define this to make all operands be large enough for Schoenhage gcd
39 to be used. */
40 #ifndef WHACK_SCHOENHAGE
41 #define WHACK_SCHOENHAGE 0
42 #endif
43
44 #if WHACK_SCHOENHAGE
45 #define MIN_OPERAND_BITSIZE (MAX_SCHOENHAGE_THRESHOLD * GMP_NUMB_BITS)
46 #else
47 #define MIN_OPERAND_BITSIZE 1
48 #endif
49
50
51 void
check_data(void)52 check_data (void)
53 {
54 static const struct {
55 const char *a;
56 const char *b;
57 const char *want;
58 } data[] = {
59 /* This tickled a bug in gmp 4.1.2 mpn/x86/k6/gcd_finda.asm. */
60 { "0x3FFC000007FFFFFFFFFF00000000003F83FFFFFFFFFFFFFFF80000000000000001",
61 "0x1FFE0007FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFC000000000000000000000001",
62 "5" }
63 };
64
65 mpz_t a, b, got, want;
66 int i;
67
68 mpz_inits (a, b, got, want, NULL);
69
70 for (i = 0; i < numberof (data); i++)
71 {
72 mpz_set_str_or_abort (a, data[i].a, 0);
73 mpz_set_str_or_abort (b, data[i].b, 0);
74 mpz_set_str_or_abort (want, data[i].want, 0);
75 mpz_gcd (got, a, b);
76 MPZ_CHECK_FORMAT (got);
77 if (mpz_cmp (got, want) != 0)
78 {
79 printf ("mpz_gcd wrong on data[%d]\n", i);
80 printf (" a %s\n", data[i].a);
81 printf (" b %s\n", data[i].b);
82 mpz_trace (" a", a);
83 mpz_trace (" b", b);
84 mpz_trace (" want", want);
85 mpz_trace (" got ", got);
86 abort ();
87 }
88 }
89
90 mpz_clears (a, b, got, want, NULL);
91 }
92
93 void
make_chain_operands(mpz_t ref,mpz_t a,mpz_t b,gmp_randstate_t rs,int nb1,int nb2,int chain_len)94 make_chain_operands (mpz_t ref, mpz_t a, mpz_t b, gmp_randstate_t rs, int nb1, int nb2, int chain_len)
95 {
96 mpz_t bs, temp1, temp2;
97 int j;
98
99 mpz_inits (bs, temp1, temp2, NULL);
100
101 /* Generate a division chain backwards, allowing otherwise unlikely huge
102 quotients. */
103
104 mpz_set_ui (a, 0);
105 mpz_urandomb (bs, rs, 32);
106 mpz_urandomb (bs, rs, mpz_get_ui (bs) % nb1 + 1);
107 mpz_rrandomb (b, rs, mpz_get_ui (bs));
108 mpz_add_ui (b, b, 1);
109 mpz_set (ref, b);
110
111 for (j = 0; j < chain_len; j++)
112 {
113 mpz_urandomb (bs, rs, 32);
114 mpz_urandomb (bs, rs, mpz_get_ui (bs) % nb2 + 1);
115 mpz_rrandomb (temp2, rs, mpz_get_ui (bs) + 1);
116 mpz_add_ui (temp2, temp2, 1);
117 mpz_mul (temp1, b, temp2);
118 mpz_add (a, a, temp1);
119
120 mpz_urandomb (bs, rs, 32);
121 mpz_urandomb (bs, rs, mpz_get_ui (bs) % nb2 + 1);
122 mpz_rrandomb (temp2, rs, mpz_get_ui (bs) + 1);
123 mpz_add_ui (temp2, temp2, 1);
124 mpz_mul (temp1, a, temp2);
125 mpz_add (b, b, temp1);
126 }
127
128 mpz_clears (bs, temp1, temp2, NULL);
129 }
130
131 /* Test operands from a table of seed data. This variant creates the operands
132 using plain ol' mpz_rrandomb. This is a hack for better coverage of the gcd
133 code, which depends on that the random number generators give the exact
134 numbers we expect. */
135 void
check_kolmo1(void)136 check_kolmo1 (void)
137 {
138 static const struct {
139 unsigned int seed;
140 int nb;
141 const char *want;
142 } data[] = {
143 { 59618, 38208, "5"},
144 { 76521, 49024, "3"},
145 { 85869, 54976, "1"},
146 { 99449, 63680, "1"},
147 {112453, 72000, "1"}
148 };
149
150 gmp_randstate_t rs;
151 mpz_t bs, a, b, want;
152 int i, unb, vnb, nb;
153
154 gmp_randinit_default (rs);
155
156 mpz_inits (bs, a, b, want, NULL);
157
158 for (i = 0; i < numberof (data); i++)
159 {
160 nb = data[i].nb;
161
162 gmp_randseed_ui (rs, data[i].seed);
163
164 mpz_urandomb (bs, rs, 32);
165 unb = mpz_get_ui (bs) % nb;
166 mpz_urandomb (bs, rs, 32);
167 vnb = mpz_get_ui (bs) % nb;
168
169 mpz_rrandomb (a, rs, unb);
170 mpz_rrandomb (b, rs, vnb);
171
172 mpz_set_str_or_abort (want, data[i].want, 0);
173
174 one_test (a, b, want, -1);
175 }
176
177 mpz_clears (bs, a, b, want, NULL);
178 gmp_randclear (rs);
179 }
180
181 /* Test operands from a table of seed data. This variant creates the operands
182 using a division chain. This is a hack for better coverage of the gcd
183 code, which depends on that the random number generators give the exact
184 numbers we expect. */
185 void
check_kolmo2(void)186 check_kolmo2 (void)
187 {
188 static const struct {
189 unsigned int seed;
190 int nb, chain_len;
191 } data[] = {
192 { 917, 15, 5 },
193 { 1032, 18, 6 },
194 { 1167, 18, 6 },
195 { 1174, 18, 6 },
196 { 1192, 18, 6 },
197 };
198
199 gmp_randstate_t rs;
200 mpz_t bs, a, b, want;
201 int i;
202
203 gmp_randinit_default (rs);
204
205 mpz_inits (bs, a, b, want, NULL);
206
207 for (i = 0; i < numberof (data); i++)
208 {
209 gmp_randseed_ui (rs, data[i].seed);
210 make_chain_operands (want, a, b, rs, data[i].nb, data[i].nb, data[i].chain_len);
211 one_test (a, b, want, -1);
212 }
213
214 mpz_clears (bs, a, b, want, NULL);
215 gmp_randclear (rs);
216 }
217
218 int
main(int argc,char ** argv)219 main (int argc, char **argv)
220 {
221 mpz_t op1, op2, ref;
222 int i, chain_len;
223 gmp_randstate_ptr rands;
224 mpz_t bs;
225 unsigned long bsi, size_range;
226 long int reps = 200;
227
228 tests_start ();
229 TESTS_REPS (reps, argv, argc);
230
231 rands = RANDS;
232
233 mpz_inits (bs, op1, op2, ref, gcd1, gcd2, temp1, temp2, temp3, s, NULL);
234
235 check_data ();
236 check_kolmo1 ();
237 check_kolmo2 ();
238
239 /* Testcase to exercise the u0 == u1 case in mpn_gcdext_lehmer_n. */
240 mpz_set_ui (op2, GMP_NUMB_MAX); /* FIXME: Huge limb doesn't always fit */
241 mpz_mul_2exp (op1, op2, 100);
242 mpz_add (op1, op1, op2);
243 mpz_mul_ui (op2, op2, 2);
244 one_test (op1, op2, NULL, -1);
245
246 for (i = 0; i < reps; i++)
247 {
248 /* Generate plain operands with unknown gcd. These types of operands
249 have proven to trigger certain bugs in development versions of the
250 gcd code. The "hgcd->row[3].rsize > M" ASSERT is not triggered by
251 the division chain code below, but that is most likely just a result
252 of that other ASSERTs are triggered before it. */
253
254 mpz_urandomb (bs, rands, 32);
255 size_range = mpz_get_ui (bs) % 17 + 2;
256
257 mpz_urandomb (bs, rands, size_range);
258 mpz_rrandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE);
259 mpz_urandomb (bs, rands, size_range);
260 mpz_rrandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE);
261
262 mpz_urandomb (bs, rands, 8);
263 bsi = mpz_get_ui (bs);
264
265 if ((bsi & 0x3c) == 4)
266 mpz_mul (op1, op1, op2); /* make op1 a multiple of op2 */
267 else if ((bsi & 0x3c) == 8)
268 mpz_mul (op2, op1, op2); /* make op2 a multiple of op1 */
269
270 if ((bsi & 1) != 0)
271 mpz_neg (op1, op1);
272 if ((bsi & 2) != 0)
273 mpz_neg (op2, op2);
274
275 one_test (op1, op2, NULL, i);
276
277 /* Generate a division chain backwards, allowing otherwise unlikely huge
278 quotients. */
279
280 mpz_urandomb (bs, rands, 32);
281 chain_len = mpz_get_ui (bs) % LOG2C (GMP_NUMB_BITS * MAX_SCHOENHAGE_THRESHOLD);
282 mpz_urandomb (bs, rands, 32);
283 chain_len = mpz_get_ui (bs) % (1 << chain_len) / 32;
284
285 make_chain_operands (ref, op1, op2, rands, 16, 12, chain_len);
286
287 one_test (op1, op2, ref, i);
288 }
289
290 /* Check that we can use NULL as first argument of mpz_gcdext. */
291 mpz_set_si (op1, -10);
292 mpz_set_si (op2, 0);
293 mpz_gcdext (NULL, temp1, temp2, op1, op2);
294 ASSERT_ALWAYS (mpz_cmp_si (temp1, -1) == 0);
295 ASSERT_ALWAYS (mpz_cmp_si (temp2, 0) == 0);
296 mpz_set_si (op2, 6);
297 mpz_gcdext (NULL, temp1, temp2, op1, op2);
298 ASSERT_ALWAYS (mpz_cmp_si (temp1, 1) == 0);
299 ASSERT_ALWAYS (mpz_cmp_si (temp2, 2) == 0);
300
301 mpz_clears (bs, op1, op2, ref, gcd1, gcd2, temp1, temp2, temp3, s, NULL);
302
303 tests_end ();
304 exit (0);
305 }
306
307 void
debug_mp(mpz_t x,int base)308 debug_mp (mpz_t x, int base)
309 {
310 mpz_out_str (stderr, base, x); fputc ('\n', stderr);
311 }
312
313 void
one_test(mpz_t op1,mpz_t op2,mpz_t ref,int i)314 one_test (mpz_t op1, mpz_t op2, mpz_t ref, int i)
315 {
316 /*
317 printf ("%d %d %d\n", SIZ (op1), SIZ (op2), ref != NULL ? SIZ (ref) : 0);
318 fflush (stdout);
319 */
320
321 /*
322 fprintf (stderr, "op1="); debug_mp (op1, -16);
323 fprintf (stderr, "op2="); debug_mp (op2, -16);
324 */
325
326 mpz_gcdext (gcd1, s, NULL, op1, op2);
327 MPZ_CHECK_FORMAT (gcd1);
328 MPZ_CHECK_FORMAT (s);
329
330 if (ref && mpz_cmp (ref, gcd1) != 0)
331 {
332 fprintf (stderr, "ERROR in test %d\n", i);
333 fprintf (stderr, "mpz_gcdext returned incorrect result\n");
334 fprintf (stderr, "op1="); debug_mp (op1, -16);
335 fprintf (stderr, "op2="); debug_mp (op2, -16);
336 fprintf (stderr, "expected result:\n"); debug_mp (ref, -16);
337 fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd1, -16);
338 abort ();
339 }
340
341 if (!gcdext_valid_p(op1, op2, gcd1, s))
342 {
343 fprintf (stderr, "ERROR in test %d\n", i);
344 fprintf (stderr, "mpz_gcdext returned invalid result\n");
345 fprintf (stderr, "op1="); debug_mp (op1, -16);
346 fprintf (stderr, "op2="); debug_mp (op2, -16);
347 fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd1, -16);
348 fprintf (stderr, "s="); debug_mp (s, -16);
349 abort ();
350 }
351
352 mpz_gcd (gcd2, op1, op2);
353 MPZ_CHECK_FORMAT (gcd2);
354
355 if (mpz_cmp (gcd2, gcd1) != 0)
356 {
357 fprintf (stderr, "ERROR in test %d\n", i);
358 fprintf (stderr, "mpz_gcd returned incorrect result\n");
359 fprintf (stderr, "op1="); debug_mp (op1, -16);
360 fprintf (stderr, "op2="); debug_mp (op2, -16);
361 fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16);
362 fprintf (stderr, "mpz_gcd returns:\n"); debug_mp (gcd2, -16);
363 abort ();
364 }
365
366 /* This should probably move to t-gcd_ui.c */
367 if (mpz_fits_ulong_p (op1) || mpz_fits_ulong_p (op2))
368 {
369 if (mpz_fits_ulong_p (op1))
370 mpz_gcd_ui (gcd2, op2, mpz_get_ui (op1));
371 else
372 mpz_gcd_ui (gcd2, op1, mpz_get_ui (op2));
373 if (mpz_cmp (gcd2, gcd1))
374 {
375 fprintf (stderr, "ERROR in test %d\n", i);
376 fprintf (stderr, "mpz_gcd_ui returned incorrect result\n");
377 fprintf (stderr, "op1="); debug_mp (op1, -16);
378 fprintf (stderr, "op2="); debug_mp (op2, -16);
379 fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16);
380 fprintf (stderr, "mpz_gcd_ui returns:\n"); debug_mp (gcd2, -16);
381 abort ();
382 }
383 }
384
385 mpz_gcdext (gcd2, temp1, temp2, op1, op2);
386 MPZ_CHECK_FORMAT (gcd2);
387 MPZ_CHECK_FORMAT (temp1);
388 MPZ_CHECK_FORMAT (temp2);
389
390 mpz_mul (temp1, temp1, op1);
391 mpz_mul (temp2, temp2, op2);
392 mpz_add (temp1, temp1, temp2);
393
394 if (mpz_cmp (gcd1, gcd2) != 0
395 || mpz_cmp (gcd2, temp1) != 0)
396 {
397 fprintf (stderr, "ERROR in test %d\n", i);
398 fprintf (stderr, "mpz_gcdext returned incorrect result\n");
399 fprintf (stderr, "op1="); debug_mp (op1, -16);
400 fprintf (stderr, "op2="); debug_mp (op2, -16);
401 fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16);
402 fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd2, -16);
403 abort ();
404 }
405 }
406
407 /* Called when g is supposed to be gcd(a,b), and g = s a + t b, for some t.
408 Uses temp1, temp2 and temp3. */
409 static int
gcdext_valid_p(const mpz_t a,const mpz_t b,const mpz_t g,const mpz_t s)410 gcdext_valid_p (const mpz_t a, const mpz_t b, const mpz_t g, const mpz_t s)
411 {
412 /* It's not clear that gcd(0,0) is well defined, but we allow it and require that
413 gcd(0,0) = 0. */
414 if (mpz_sgn (g) < 0)
415 return 0;
416
417 if (mpz_sgn (a) == 0)
418 {
419 /* Must have g == abs (b). Any value for s is in some sense "correct",
420 but it makes sense to require that s == 0. */
421 return mpz_cmpabs (g, b) == 0 && mpz_sgn (s) == 0;
422 }
423 else if (mpz_sgn (b) == 0)
424 {
425 /* Must have g == abs (a), s == sign (a) */
426 return mpz_cmpabs (g, a) == 0 && mpz_cmp_si (s, mpz_sgn (a)) == 0;
427 }
428
429 if (mpz_sgn (g) <= 0)
430 return 0;
431
432 mpz_tdiv_qr (temp1, temp3, a, g);
433 if (mpz_sgn (temp3) != 0)
434 return 0;
435
436 mpz_tdiv_qr (temp2, temp3, b, g);
437 if (mpz_sgn (temp3) != 0)
438 return 0;
439
440 /* Require that 2 |s| < |b/g|, or |s| == 1. */
441 if (mpz_cmpabs_ui (s, 1) > 0)
442 {
443 mpz_mul_2exp (temp3, s, 1);
444 if (mpz_cmpabs (temp3, temp2) >= 0)
445 return 0;
446 }
447
448 /* Compute the other cofactor. */
449 mpz_mul(temp2, s, a);
450 mpz_sub(temp2, g, temp2);
451 mpz_tdiv_qr(temp2, temp3, temp2, b);
452
453 if (mpz_sgn (temp3) != 0)
454 return 0;
455
456 /* Require that 2 |t| < |a/g| or |t| == 1*/
457 if (mpz_cmpabs_ui (temp2, 1) > 0)
458 {
459 mpz_mul_2exp (temp2, temp2, 1);
460 if (mpz_cmpabs (temp2, temp1) >= 0)
461 return 0;
462 }
463 return 1;
464 }
465