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
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
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
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
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
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_set_ui (op2, 0);
242 mpz_setbit (op2, GMP_NUMB_BITS);
243 mpz_sub_ui (op2, op2, 1);
244 mpz_mul_2exp (op1, op2, 100);
245 mpz_add (op1, op1, op2);
246 mpz_mul_ui (op2, op2, 2);
247 one_test (op1, op2, NULL, -1);
248
249 for (i = 0; i < reps; i++)
250 {
251 /* Generate plain operands with unknown gcd. These types of operands
252 have proven to trigger certain bugs in development versions of the
253 gcd code. The "hgcd->row[3].rsize > M" ASSERT is not triggered by
254 the division chain code below, but that is most likely just a result
255 of that other ASSERTs are triggered before it. */
256
257 mpz_urandomb (bs, rands, 32);
258 size_range = mpz_get_ui (bs) % 17 + 2;
259
260 mpz_urandomb (bs, rands, size_range);
261 mpz_rrandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE);
262 mpz_urandomb (bs, rands, size_range);
263 mpz_rrandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_BITSIZE);
264
265 mpz_urandomb (bs, rands, 8);
266 bsi = mpz_get_ui (bs);
267
268 if ((bsi & 0x3c) == 4)
269 mpz_mul (op1, op1, op2); /* make op1 a multiple of op2 */
270 else if ((bsi & 0x3c) == 8)
271 mpz_mul (op2, op1, op2); /* make op2 a multiple of op1 */
272
273 if ((bsi & 1) != 0)
274 mpz_neg (op1, op1);
275 if ((bsi & 2) != 0)
276 mpz_neg (op2, op2);
277
278 one_test (op1, op2, NULL, i);
279
280 /* Generate a division chain backwards, allowing otherwise unlikely huge
281 quotients. */
282
283 mpz_urandomb (bs, rands, 32);
284 chain_len = mpz_get_ui (bs) % LOG2C (GMP_NUMB_BITS * MAX_SCHOENHAGE_THRESHOLD);
285 mpz_urandomb (bs, rands, 32);
286 chain_len = mpz_get_ui (bs) % (1 << chain_len) / 32;
287
288 make_chain_operands (ref, op1, op2, rands, 16, 12, chain_len);
289
290 one_test (op1, op2, ref, i);
291 }
292
293 /* Check that we can use NULL as first argument of mpz_gcdext. */
294 mpz_set_si (op1, -10);
295 mpz_set_si (op2, 0);
296 mpz_gcdext (NULL, temp1, temp2, op1, op2);
297 ASSERT_ALWAYS (mpz_cmp_si (temp1, -1) == 0);
298 ASSERT_ALWAYS (mpz_cmp_si (temp2, 0) == 0);
299 mpz_set_si (op2, 6);
300 mpz_gcdext (NULL, temp1, temp2, op1, op2);
301 ASSERT_ALWAYS (mpz_cmp_si (temp1, 1) == 0);
302 ASSERT_ALWAYS (mpz_cmp_si (temp2, 2) == 0);
303
304 mpz_clears (bs, op1, op2, ref, gcd1, gcd2, temp1, temp2, temp3, s, NULL);
305
306 tests_end ();
307 exit (0);
308 }
309
310 void
311 debug_mp (mpz_t x, int base)
312 {
313 mpz_out_str (stderr, base, x); fputc ('\n', stderr);
314 }
315
316 void
317 one_test (mpz_t op1, mpz_t op2, mpz_t ref, int i)
318 {
319 /*
320 printf ("%d %d %d\n", SIZ (op1), SIZ (op2), ref != NULL ? SIZ (ref) : 0);
321 fflush (stdout);
322 */
323
324 /*
325 fprintf (stderr, "op1="); debug_mp (op1, -16);
326 fprintf (stderr, "op2="); debug_mp (op2, -16);
327 */
328
329 mpz_gcdext (gcd1, s, NULL, op1, op2);
330 MPZ_CHECK_FORMAT (gcd1);
331 MPZ_CHECK_FORMAT (s);
332
333 if (ref && mpz_cmp (ref, gcd1) != 0)
334 {
335 fprintf (stderr, "ERROR in test %d\n", i);
336 fprintf (stderr, "mpz_gcdext returned incorrect result\n");
337 fprintf (stderr, "op1="); debug_mp (op1, -16);
338 fprintf (stderr, "op2="); debug_mp (op2, -16);
339 fprintf (stderr, "expected result:\n"); debug_mp (ref, -16);
340 fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd1, -16);
341 abort ();
342 }
343
344 if (!gcdext_valid_p(op1, op2, gcd1, s))
345 {
346 fprintf (stderr, "ERROR in test %d\n", i);
347 fprintf (stderr, "mpz_gcdext returned invalid result\n");
348 fprintf (stderr, "op1="); debug_mp (op1, -16);
349 fprintf (stderr, "op2="); debug_mp (op2, -16);
350 fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd1, -16);
351 fprintf (stderr, "s="); debug_mp (s, -16);
352 abort ();
353 }
354
355 mpz_gcd (gcd2, op1, op2);
356 MPZ_CHECK_FORMAT (gcd2);
357
358 if (mpz_cmp (gcd2, gcd1) != 0)
359 {
360 fprintf (stderr, "ERROR in test %d\n", i);
361 fprintf (stderr, "mpz_gcd returned incorrect result\n");
362 fprintf (stderr, "op1="); debug_mp (op1, -16);
363 fprintf (stderr, "op2="); debug_mp (op2, -16);
364 fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16);
365 fprintf (stderr, "mpz_gcd returns:\n"); debug_mp (gcd2, -16);
366 abort ();
367 }
368
369 /* This should probably move to t-gcd_ui.c */
370 if (mpz_fits_ulong_p (op1) || mpz_fits_ulong_p (op2))
371 {
372 if (mpz_fits_ulong_p (op1))
373 mpz_gcd_ui (gcd2, op2, mpz_get_ui (op1));
374 else
375 mpz_gcd_ui (gcd2, op1, mpz_get_ui (op2));
376 if (mpz_cmp (gcd2, gcd1))
377 {
378 fprintf (stderr, "ERROR in test %d\n", i);
379 fprintf (stderr, "mpz_gcd_ui returned incorrect result\n");
380 fprintf (stderr, "op1="); debug_mp (op1, -16);
381 fprintf (stderr, "op2="); debug_mp (op2, -16);
382 fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16);
383 fprintf (stderr, "mpz_gcd_ui returns:\n"); debug_mp (gcd2, -16);
384 abort ();
385 }
386 }
387
388 mpz_gcdext (gcd2, temp1, temp2, op1, op2);
389 MPZ_CHECK_FORMAT (gcd2);
390 MPZ_CHECK_FORMAT (temp1);
391 MPZ_CHECK_FORMAT (temp2);
392
393 mpz_mul (temp1, temp1, op1);
394 mpz_mul (temp2, temp2, op2);
395 mpz_add (temp1, temp1, temp2);
396
397 if (mpz_cmp (gcd1, gcd2) != 0
398 || mpz_cmp (gcd2, temp1) != 0)
399 {
400 fprintf (stderr, "ERROR in test %d\n", i);
401 fprintf (stderr, "mpz_gcdext returned incorrect result\n");
402 fprintf (stderr, "op1="); debug_mp (op1, -16);
403 fprintf (stderr, "op2="); debug_mp (op2, -16);
404 fprintf (stderr, "expected result:\n"); debug_mp (gcd1, -16);
405 fprintf (stderr, "mpz_gcdext returns:\n");debug_mp (gcd2, -16);
406 abort ();
407 }
408 }
409
410 /* Called when g is supposed to be gcd(a,b), and g = s a + t b, for some t.
411 Uses temp1, temp2 and temp3. */
412 static int
413 gcdext_valid_p (const mpz_t a, const mpz_t b, const mpz_t g, const mpz_t s)
414 {
415 /* It's not clear that gcd(0,0) is well defined, but we allow it and require that
416 gcd(0,0) = 0. */
417 if (mpz_sgn (g) < 0)
418 return 0;
419
420 if (mpz_sgn (a) == 0)
421 {
422 /* Must have g == abs (b). Any value for s is in some sense "correct",
423 but it makes sense to require that s == 0. */
424 return mpz_cmpabs (g, b) == 0 && mpz_sgn (s) == 0;
425 }
426 else if (mpz_sgn (b) == 0)
427 {
428 /* Must have g == abs (a), s == sign (a) */
429 return mpz_cmpabs (g, a) == 0 && mpz_cmp_si (s, mpz_sgn (a)) == 0;
430 }
431
432 if (mpz_sgn (g) <= 0)
433 return 0;
434
435 mpz_tdiv_qr (temp1, temp3, a, g);
436 if (mpz_sgn (temp3) != 0)
437 return 0;
438
439 mpz_tdiv_qr (temp2, temp3, b, g);
440 if (mpz_sgn (temp3) != 0)
441 return 0;
442
443 /* Require that 2 |s| < |b/g|, or |s| == 1. */
444 if (mpz_cmpabs_ui (s, 1) > 0)
445 {
446 mpz_mul_2exp (temp3, s, 1);
447 if (mpz_cmpabs (temp3, temp2) >= 0)
448 return 0;
449 }
450
451 /* Compute the other cofactor. */
452 mpz_mul(temp2, s, a);
453 mpz_sub(temp2, g, temp2);
454 mpz_tdiv_qr(temp2, temp3, temp2, b);
455
456 if (mpz_sgn (temp3) != 0)
457 return 0;
458
459 /* Require that 2 |t| < |a/g| or |t| == 1*/
460 if (mpz_cmpabs_ui (temp2, 1) > 0)
461 {
462 mpz_mul_2exp (temp2, temp2, 1);
463 if (mpz_cmpabs (temp2, temp1) >= 0)
464 return 0;
465 }
466 return 1;
467 }