1 /* Test mpn_hgcd_appr.
2
3 Copyright 1991, 1993, 1994, 1996, 1997, 2000-2004, 2011 Free Software
4 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 #include <string.h>
24
25 #include "gmp-impl.h"
26 #include "tests.h"
27
28 static mp_size_t one_test (mpz_t, mpz_t, int);
29 static void debug_mp (mpz_t, int);
30
31 #define MIN_OPERAND_SIZE 2
32
33 struct hgcd_ref
34 {
35 mpz_t m[2][2];
36 };
37
38 static void hgcd_ref_init (struct hgcd_ref *hgcd);
39 static void hgcd_ref_clear (struct hgcd_ref *hgcd);
40 static int hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b);
41 static int hgcd_ref_equal (const struct hgcd_ref *, const struct hgcd_ref *);
42 static int hgcd_appr_valid_p (mpz_t, mpz_t, mp_size_t, struct hgcd_ref *,
43 mpz_t, mpz_t, mp_size_t, struct hgcd_matrix *);
44
45 static int verbose_flag = 0;
46
47 int
48 main (int argc, char **argv)
49 {
50 mpz_t op1, op2, temp1, temp2;
51 int i, j, chain_len;
52 gmp_randstate_ptr rands;
53 mpz_t bs;
54 unsigned long size_range;
55
56 if (argc > 1)
57 {
58 if (strcmp (argv[1], "-v") == 0)
59 verbose_flag = 1;
60 else
61 {
62 fprintf (stderr, "Invalid argument.\n");
63 return 1;
64 }
65 }
66
67 tests_start ();
68 rands = RANDS;
69
70 mpz_init (bs);
71 mpz_init (op1);
72 mpz_init (op2);
73 mpz_init (temp1);
74 mpz_init (temp2);
75
76 for (i = 0; i < 15; i++)
77 {
78 /* Generate plain operands with unknown gcd. These types of operands
79 have proven to trigger certain bugs in development versions of the
80 gcd code. */
81
82 mpz_urandomb (bs, rands, 32);
83 size_range = mpz_get_ui (bs) % 13 + 2;
84
85 mpz_urandomb (bs, rands, size_range);
86 mpz_urandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE);
87 mpz_urandomb (bs, rands, size_range);
88 mpz_urandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE);
89
90 if (mpz_cmp (op1, op2) < 0)
91 mpz_swap (op1, op2);
92
93 if (mpz_size (op1) > 0)
94 one_test (op1, op2, i);
95
96 /* Generate a division chain backwards, allowing otherwise
97 unlikely huge quotients. */
98
99 mpz_set_ui (op1, 0);
100 mpz_urandomb (bs, rands, 32);
101 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 16 + 1);
102 mpz_rrandomb (op2, rands, mpz_get_ui (bs));
103 mpz_add_ui (op2, op2, 1);
104
105 #if 0
106 chain_len = 1000000;
107 #else
108 mpz_urandomb (bs, rands, 32);
109 chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * GCD_DC_THRESHOLD / 256);
110 #endif
111
112 for (j = 0; j < chain_len; j++)
113 {
114 mpz_urandomb (bs, rands, 32);
115 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
116 mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
117 mpz_add_ui (temp2, temp2, 1);
118 mpz_mul (temp1, op2, temp2);
119 mpz_add (op1, op1, temp1);
120
121 /* Don't generate overly huge operands. */
122 if (SIZ (op1) > 3 * GCD_DC_THRESHOLD)
123 break;
124
125 mpz_urandomb (bs, rands, 32);
126 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
127 mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
128 mpz_add_ui (temp2, temp2, 1);
129 mpz_mul (temp1, op1, temp2);
130 mpz_add (op2, op2, temp1);
131
132 /* Don't generate overly huge operands. */
133 if (SIZ (op2) > 3 * GCD_DC_THRESHOLD)
134 break;
135 }
136 if (mpz_cmp (op1, op2) < 0)
137 mpz_swap (op1, op2);
138
139 if (mpz_size (op1) > 0)
140 one_test (op1, op2, i);
141 }
142
143 mpz_clear (bs);
144 mpz_clear (op1);
145 mpz_clear (op2);
146 mpz_clear (temp1);
147 mpz_clear (temp2);
148
149 tests_end ();
150 exit (0);
151 }
152
153 static void
154 debug_mp (mpz_t x, int base)
155 {
156 mpz_out_str (stderr, base, x); fputc ('\n', stderr);
157 }
158
159 static mp_size_t
160 one_test (mpz_t a, mpz_t b, int i)
161 {
162 struct hgcd_matrix hgcd;
163 struct hgcd_ref ref;
164
165 mpz_t ref_r0;
166 mpz_t ref_r1;
167 mpz_t hgcd_r0;
168 mpz_t hgcd_r1;
169
170 int res[2];
171 mp_size_t asize;
172 mp_size_t bsize;
173
174 mp_size_t hgcd_init_scratch;
175 mp_size_t hgcd_scratch;
176
177 mp_ptr hgcd_init_tp;
178 mp_ptr hgcd_tp;
179 mp_limb_t marker[4];
180
181 asize = a->_mp_size;
182 bsize = b->_mp_size;
183
184 ASSERT (asize >= bsize);
185
186 hgcd_init_scratch = MPN_HGCD_MATRIX_INIT_ITCH (asize);
187 hgcd_init_tp = refmpn_malloc_limbs (hgcd_init_scratch + 2) + 1;
188 mpn_hgcd_matrix_init (&hgcd, asize, hgcd_init_tp);
189
190 hgcd_scratch = mpn_hgcd_appr_itch (asize);
191 hgcd_tp = refmpn_malloc_limbs (hgcd_scratch + 2) + 1;
192
193 mpn_random (marker, 4);
194
195 hgcd_init_tp[-1] = marker[0];
196 hgcd_init_tp[hgcd_init_scratch] = marker[1];
197 hgcd_tp[-1] = marker[2];
198 hgcd_tp[hgcd_scratch] = marker[3];
199
200 #if 0
201 fprintf (stderr,
202 "one_test: i = %d asize = %d, bsize = %d\n",
203 i, a->_mp_size, b->_mp_size);
204
205 gmp_fprintf (stderr,
206 "one_test: i = %d\n"
207 " a = %Zx\n"
208 " b = %Zx\n",
209 i, a, b);
210 #endif
211 hgcd_ref_init (&ref);
212
213 mpz_init_set (ref_r0, a);
214 mpz_init_set (ref_r1, b);
215 res[0] = hgcd_ref (&ref, ref_r0, ref_r1);
216
217 mpz_init_set (hgcd_r0, a);
218 mpz_init_set (hgcd_r1, b);
219 if (bsize < asize)
220 {
221 _mpz_realloc (hgcd_r1, asize);
222 MPN_ZERO (hgcd_r1->_mp_d + bsize, asize - bsize);
223 }
224 res[1] = mpn_hgcd_appr (hgcd_r0->_mp_d,
225 hgcd_r1->_mp_d,
226 asize,
227 &hgcd, hgcd_tp);
228
229 if (hgcd_init_tp[-1] != marker[0]
230 || hgcd_init_tp[hgcd_init_scratch] != marker[1]
231 || hgcd_tp[-1] != marker[2]
232 || hgcd_tp[hgcd_scratch] != marker[3])
233 {
234 fprintf (stderr, "ERROR in test %d\n", i);
235 fprintf (stderr, "scratch space overwritten!\n");
236
237 if (hgcd_init_tp[-1] != marker[0])
238 gmp_fprintf (stderr,
239 "before init_tp: %Mx\n"
240 "expected: %Mx\n",
241 hgcd_init_tp[-1], marker[0]);
242 if (hgcd_init_tp[hgcd_init_scratch] != marker[1])
243 gmp_fprintf (stderr,
244 "after init_tp: %Mx\n"
245 "expected: %Mx\n",
246 hgcd_init_tp[hgcd_init_scratch], marker[1]);
247 if (hgcd_tp[-1] != marker[2])
248 gmp_fprintf (stderr,
249 "before tp: %Mx\n"
250 "expected: %Mx\n",
251 hgcd_tp[-1], marker[2]);
252 if (hgcd_tp[hgcd_scratch] != marker[3])
253 gmp_fprintf (stderr,
254 "after tp: %Mx\n"
255 "expected: %Mx\n",
256 hgcd_tp[hgcd_scratch], marker[3]);
257
258 abort ();
259 }
260
261 if (!hgcd_appr_valid_p (a, b, res[0], &ref, ref_r0, ref_r1,
262 res[1], &hgcd))
263 {
264 fprintf (stderr, "ERROR in test %d\n", i);
265 fprintf (stderr, "Invalid results for hgcd and hgcd_ref\n");
266 fprintf (stderr, "op1="); debug_mp (a, -16);
267 fprintf (stderr, "op2="); debug_mp (b, -16);
268 fprintf (stderr, "hgcd_ref: %ld\n", (long) res[0]);
269 fprintf (stderr, "mpn_hgcd_appr: %ld\n", (long) res[1]);
270 abort ();
271 }
272
273 refmpn_free_limbs (hgcd_init_tp - 1);
274 refmpn_free_limbs (hgcd_tp - 1);
275 hgcd_ref_clear (&ref);
276 mpz_clear (ref_r0);
277 mpz_clear (ref_r1);
278 mpz_clear (hgcd_r0);
279 mpz_clear (hgcd_r1);
280
281 return res[0];
282 }
283
284 static void
285 hgcd_ref_init (struct hgcd_ref *hgcd)
286 {
287 unsigned i;
288 for (i = 0; i<2; i++)
289 {
290 unsigned j;
291 for (j = 0; j<2; j++)
292 mpz_init (hgcd->m[i][j]);
293 }
294 }
295
296 static void
297 hgcd_ref_clear (struct hgcd_ref *hgcd)
298 {
299 unsigned i;
300 for (i = 0; i<2; i++)
301 {
302 unsigned j;
303 for (j = 0; j<2; j++)
304 mpz_clear (hgcd->m[i][j]);
305 }
306 }
307
308 static int
309 sdiv_qr (mpz_t q, mpz_t r, mp_size_t s, const mpz_t a, const mpz_t b)
310 {
311 mpz_fdiv_qr (q, r, a, b);
312 if (mpz_size (r) <= s)
313 {
314 mpz_add (r, r, b);
315 mpz_sub_ui (q, q, 1);
316 }
317
318 return (mpz_sgn (q) > 0);
319 }
320
321 static int
322 hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b)
323 {
324 mp_size_t n = MAX (mpz_size (a), mpz_size (b));
325 mp_size_t s = n/2 + 1;
326 mpz_t q;
327 int res;
328
329 if (mpz_size (a) <= s || mpz_size (b) <= s)
330 return 0;
331
332 res = mpz_cmp (a, b);
333 if (res < 0)
334 {
335 mpz_sub (b, b, a);
336 if (mpz_size (b) <= s)
337 return 0;
338
339 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 0);
340 mpz_set_ui (hgcd->m[1][0], 1); mpz_set_ui (hgcd->m[1][1], 1);
341 }
342 else if (res > 0)
343 {
344 mpz_sub (a, a, b);
345 if (mpz_size (a) <= s)
346 return 0;
347
348 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 1);
349 mpz_set_ui (hgcd->m[1][0], 0); mpz_set_ui (hgcd->m[1][1], 1);
350 }
351 else
352 return 0;
353
354 mpz_init (q);
355
356 for (;;)
357 {
358 ASSERT (mpz_size (a) > s);
359 ASSERT (mpz_size (b) > s);
360
361 if (mpz_cmp (a, b) > 0)
362 {
363 if (!sdiv_qr (q, a, s, a, b))
364 break;
365 mpz_addmul (hgcd->m[0][1], q, hgcd->m[0][0]);
366 mpz_addmul (hgcd->m[1][1], q, hgcd->m[1][0]);
367 }
368 else
369 {
370 if (!sdiv_qr (q, b, s, b, a))
371 break;
372 mpz_addmul (hgcd->m[0][0], q, hgcd->m[0][1]);
373 mpz_addmul (hgcd->m[1][0], q, hgcd->m[1][1]);
374 }
375 }
376
377 mpz_clear (q);
378
379 return 1;
380 }
381
382 static int
383 hgcd_ref_equal (const struct hgcd_ref *A, const struct hgcd_ref *B)
384 {
385 unsigned i;
386
387 for (i = 0; i<2; i++)
388 {
389 unsigned j;
390
391 for (j = 0; j<2; j++)
392 if (mpz_cmp (A->m[i][j], B->m[i][j]) != 0)
393 return 0;
394 }
395
396 return 1;
397 }
398
399 static int
400 hgcd_appr_valid_p (mpz_t a, mpz_t b, mp_size_t res0,
401 struct hgcd_ref *ref, mpz_t ref_r0, mpz_t ref_r1,
402 mp_size_t res1, struct hgcd_matrix *hgcd)
403 {
404 mp_size_t n = MAX (mpz_size (a), mpz_size (b));
405 mp_size_t s = n/2 + 1;
406
407 mp_bitcnt_t dbits, abits, margin;
408 mpz_t appr_r0, appr_r1, t, q;
409 struct hgcd_ref appr;
410
411 if (!res0)
412 {
413 if (!res1)
414 return 1;
415
416 fprintf (stderr, "mpn_hgcd_appr returned 1 when no reduction possible.\n");
417 return 0;
418 }
419
420 /* NOTE: No *_clear calls on error return, since we're going to
421 abort anyway. */
422 mpz_init (t);
423 mpz_init (q);
424 hgcd_ref_init (&appr);
425 mpz_init (appr_r0);
426 mpz_init (appr_r1);
427
428 if (mpz_size (ref_r0) <= s)
429 {
430 fprintf (stderr, "ref_r0 too small!!!: "); debug_mp (ref_r0, 16);
431 return 0;
432 }
433 if (mpz_size (ref_r1) <= s)
434 {
435 fprintf (stderr, "ref_r1 too small!!!: "); debug_mp (ref_r1, 16);
436 return 0;
437 }
438
439 mpz_sub (t, ref_r0, ref_r1);
440 dbits = mpz_sizeinbase (t, 2);
441 if (dbits > s*GMP_NUMB_BITS)
442 {
443 fprintf (stderr, "ref |r0 - r1| too large!!!: "); debug_mp (t, 16);
444 return 0;
445 }
446
447 if (!res1)
448 {
449 mpz_set (appr_r0, a);
450 mpz_set (appr_r1, b);
451 }
452 else
453 {
454 unsigned i;
455
456 for (i = 0; i<2; i++)
457 {
458 unsigned j;
459
460 for (j = 0; j<2; j++)
461 {
462 mp_size_t mn = hgcd->n;
463 MPN_NORMALIZE (hgcd->p[i][j], mn);
464 mpz_realloc (appr.m[i][j], mn);
465 MPN_COPY (PTR (appr.m[i][j]), hgcd->p[i][j], mn);
466 SIZ (appr.m[i][j]) = mn;
467 }
468 }
469 mpz_mul (appr_r0, appr.m[1][1], a);
470 mpz_mul (t, appr.m[0][1], b);
471 mpz_sub (appr_r0, appr_r0, t);
472 if (mpz_sgn (appr_r0) <= 0
473 || mpz_size (appr_r0) <= s)
474 {
475 fprintf (stderr, "appr_r0 too small: "); debug_mp (appr_r0, 16);
476 return 0;
477 }
478
479 mpz_mul (appr_r1, appr.m[1][0], a);
480 mpz_mul (t, appr.m[0][0], b);
481 mpz_sub (appr_r1, t, appr_r1);
482 if (mpz_sgn (appr_r1) <= 0
483 || mpz_size (appr_r1) <= s)
484 {
485 fprintf (stderr, "appr_r1 too small: "); debug_mp (appr_r1, 16);
486 return 0;
487 }
488 }
489
490 mpz_sub (t, appr_r0, appr_r1);
491 abits = mpz_sizeinbase (t, 2);
492 if (abits < dbits)
493 {
494 fprintf (stderr, "|r0 - r1| too small: "); debug_mp (t, 16);
495 return 0;
496 }
497
498 /* We lose one bit each time we discard the least significant limbs.
499 For the lehmer code, that can happen at most s * (GMP_NUMB_BITS)
500 / (GMP_NUMB_BITS - 1) times. For the dc code, we lose an entire
501 limb (or more?) for each level of recursion. */
502
503 margin = (n/2+1) * GMP_NUMB_BITS / (GMP_NUMB_BITS - 1);
504 {
505 mp_size_t rn;
506 for (rn = n; ABOVE_THRESHOLD (rn, HGCD_APPR_THRESHOLD); rn = (rn + 1)/2)
507 margin += GMP_NUMB_BITS;
508 }
509
510 if (verbose_flag && abits > dbits)
511 fprintf (stderr, "n = %u: sbits = %u: ref #(r0-r1): %u, appr #(r0-r1): %u excess: %d, margin: %u\n",
512 (unsigned) n, (unsigned) s*GMP_NUMB_BITS,
513 (unsigned) dbits, (unsigned) abits,
514 (int) (abits - s * GMP_NUMB_BITS), (unsigned) margin);
515
516 if (abits > s*GMP_NUMB_BITS + margin)
517 {
518 fprintf (stderr, "appr |r0 - r1| much larger than minimal (by %u bits, margin = %u bits)\n",
519 (unsigned) (abits - s*GMP_NUMB_BITS), (unsigned) margin);
520 return 0;
521 }
522
523 while (mpz_cmp (appr_r0, ref_r0) > 0 || mpz_cmp (appr_r1, ref_r1) > 0)
524 {
525 ASSERT (mpz_size (appr_r0) > s);
526 ASSERT (mpz_size (appr_r1) > s);
527
528 if (mpz_cmp (appr_r0, appr_r1) > 0)
529 {
530 if (!sdiv_qr (q, appr_r0, s, appr_r0, appr_r1))
531 break;
532 mpz_addmul (appr.m[0][1], q, appr.m[0][0]);
533 mpz_addmul (appr.m[1][1], q, appr.m[1][0]);
534 }
535 else
536 {
537 if (!sdiv_qr (q, appr_r1, s, appr_r1, appr_r0))
538 break;
539 mpz_addmul (appr.m[0][0], q, appr.m[0][1]);
540 mpz_addmul (appr.m[1][0], q, appr.m[1][1]);
541 }
542 }
543
544 if (mpz_cmp (appr_r0, ref_r0) != 0
545 || mpz_cmp (appr_r1, ref_r1) != 0
546 || !hgcd_ref_equal (ref, &appr))
547 {
548 fprintf (stderr, "appr_r0: "); debug_mp (appr_r0, 16);
549 fprintf (stderr, "ref_r0: "); debug_mp (ref_r0, 16);
550
551 fprintf (stderr, "appr_r1: "); debug_mp (appr_r1, 16);
552 fprintf (stderr, "ref_r1: "); debug_mp (ref_r1, 16);
553
554 return 0;
555 }
556 mpz_clear (t);
557 mpz_clear (q);
558 hgcd_ref_clear (&appr);
559 mpz_clear (appr_r0);
560 mpz_clear (appr_r1);
561
562 return 1;
563 }