1 /* Test mpn_hgcd.
2
3 Copyright 1991, 1993, 1994, 1996, 1997, 2000-2004 Free Software Foundation,
4 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 static mp_size_t one_test (mpz_t, mpz_t, int);
28 static void debug_mp (mpz_t, int);
29
30 #define MIN_OPERAND_SIZE 2
31
32 /* Fixed values, for regression testing of mpn_hgcd. */
33 struct value { int res; const char *a; const char *b; };
34 static const struct value hgcd_values[] = {
35 #if GMP_NUMB_BITS == 32
36 { 5,
37 "0x1bddff867272a9296ac493c251d7f46f09a5591fe",
38 "0xb55930a2a68a916450a7de006031068c5ddb0e5c" },
39 { 4,
40 "0x2f0ece5b1ee9c15e132a01d55768dc13",
41 "0x1c6f4fd9873cdb24466e6d03e1cc66e7" },
42 { 3, "0x7FFFFC003FFFFFFFFFC5", "0x3FFFFE001FFFFFFFFFE3"},
43 #endif
44 { -1, NULL, NULL }
45 };
46
47 struct hgcd_ref
48 {
49 mpz_t m[2][2];
50 };
51
52 static void hgcd_ref_init (struct hgcd_ref *);
53 static void hgcd_ref_clear (struct hgcd_ref *);
54 static int hgcd_ref (struct hgcd_ref *, mpz_t, mpz_t);
55 static int hgcd_ref_equal (const struct hgcd_matrix *, const struct hgcd_ref *);
56
57 int
58 main (int argc, char **argv)
59 {
60 mpz_t op1, op2, temp1, temp2;
61 int i, j, chain_len;
62 gmp_randstate_ptr rands;
63 mpz_t bs;
64 unsigned long size_range;
65
66 tests_start ();
67 rands = RANDS;
68
69 mpz_init (bs);
70 mpz_init (op1);
71 mpz_init (op2);
72 mpz_init (temp1);
73 mpz_init (temp2);
74
75 for (i = 0; hgcd_values[i].res >= 0; i++)
76 {
77 mp_size_t res;
78
79 mpz_set_str (op1, hgcd_values[i].a, 0);
80 mpz_set_str (op2, hgcd_values[i].b, 0);
81
82 res = one_test (op1, op2, -1-i);
83 if (res != hgcd_values[i].res)
84 {
85 fprintf (stderr, "ERROR in test %d\n", -1-i);
86 fprintf (stderr, "Bad return code from hgcd\n");
87 fprintf (stderr, "op1="); debug_mp (op1, -16);
88 fprintf (stderr, "op2="); debug_mp (op2, -16);
89 fprintf (stderr, "expected: %d\n", hgcd_values[i].res);
90 fprintf (stderr, "hgcd: %d\n", (int) res);
91 abort ();
92 }
93 }
94
95 for (i = 0; i < 15; i++)
96 {
97 /* Generate plain operands with unknown gcd. These types of operands
98 have proven to trigger certain bugs in development versions of the
99 gcd code. */
100
101 mpz_urandomb (bs, rands, 32);
102 size_range = mpz_get_ui (bs) % 13 + 2;
103
104 mpz_urandomb (bs, rands, size_range);
105 mpz_rrandomb (op1, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE);
106 mpz_urandomb (bs, rands, size_range);
107 mpz_rrandomb (op2, rands, mpz_get_ui (bs) + MIN_OPERAND_SIZE);
108
109 if (mpz_cmp (op1, op2) < 0)
110 mpz_swap (op1, op2);
111
112 if (mpz_size (op1) > 0)
113 one_test (op1, op2, i);
114
115 /* Generate a division chain backwards, allowing otherwise
116 unlikely huge quotients. */
117
118 mpz_set_ui (op1, 0);
119 mpz_urandomb (bs, rands, 32);
120 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 16 + 1);
121 mpz_rrandomb (op2, rands, mpz_get_ui (bs));
122 mpz_add_ui (op2, op2, 1);
123
124 #if 0
125 chain_len = 1000000;
126 #else
127 mpz_urandomb (bs, rands, 32);
128 chain_len = mpz_get_ui (bs) % (GMP_NUMB_BITS * GCD_DC_THRESHOLD / 256);
129 #endif
130
131 for (j = 0; j < chain_len; j++)
132 {
133 mpz_urandomb (bs, rands, 32);
134 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
135 mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
136 mpz_add_ui (temp2, temp2, 1);
137 mpz_mul (temp1, op2, temp2);
138 mpz_add (op1, op1, temp1);
139
140 /* Don't generate overly huge operands. */
141 if (SIZ (op1) > 3 * GCD_DC_THRESHOLD)
142 break;
143
144 mpz_urandomb (bs, rands, 32);
145 mpz_urandomb (bs, rands, mpz_get_ui (bs) % 12 + 1);
146 mpz_rrandomb (temp2, rands, mpz_get_ui (bs) + 1);
147 mpz_add_ui (temp2, temp2, 1);
148 mpz_mul (temp1, op1, temp2);
149 mpz_add (op2, op2, temp1);
150
151 /* Don't generate overly huge operands. */
152 if (SIZ (op2) > 3 * GCD_DC_THRESHOLD)
153 break;
154 }
155 if (mpz_cmp (op1, op2) < 0)
156 mpz_swap (op1, op2);
157
158 if (mpz_size (op1) > 0)
159 one_test (op1, op2, i);
160 }
161
162 mpz_clear (bs);
163 mpz_clear (op1);
164 mpz_clear (op2);
165 mpz_clear (temp1);
166 mpz_clear (temp2);
167
168 tests_end ();
169 exit (0);
170 }
171
172 static void
173 debug_mp (mpz_t x, int base)
174 {
175 mpz_out_str (stderr, base, x); fputc ('\n', stderr);
176 }
177
178 static int
179 mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize);
180
181 static mp_size_t
182 one_test (mpz_t a, mpz_t b, int i)
183 {
184 struct hgcd_matrix hgcd;
185 struct hgcd_ref ref;
186
187 mpz_t ref_r0;
188 mpz_t ref_r1;
189 mpz_t hgcd_r0;
190 mpz_t hgcd_r1;
191
192 mp_size_t res[2];
193 mp_size_t asize;
194 mp_size_t bsize;
195
196 mp_size_t hgcd_init_scratch;
197 mp_size_t hgcd_scratch;
198
199 mp_ptr hgcd_init_tp;
200 mp_ptr hgcd_tp;
201
202 asize = a->_mp_size;
203 bsize = b->_mp_size;
204
205 ASSERT (asize >= bsize);
206
207 hgcd_init_scratch = MPN_HGCD_MATRIX_INIT_ITCH (asize);
208 hgcd_init_tp = refmpn_malloc_limbs (hgcd_init_scratch);
209 mpn_hgcd_matrix_init (&hgcd, asize, hgcd_init_tp);
210
211 hgcd_scratch = mpn_hgcd_itch (asize);
212 hgcd_tp = refmpn_malloc_limbs (hgcd_scratch);
213
214 #if 0
215 fprintf (stderr,
216 "one_test: i = %d asize = %d, bsize = %d\n",
217 i, a->_mp_size, b->_mp_size);
218
219 gmp_fprintf (stderr,
220 "one_test: i = %d\n"
221 " a = %Zx\n"
222 " b = %Zx\n",
223 i, a, b);
224 #endif
225 hgcd_ref_init (&ref);
226
227 mpz_init_set (ref_r0, a);
228 mpz_init_set (ref_r1, b);
229 res[0] = hgcd_ref (&ref, ref_r0, ref_r1);
230
231 mpz_init_set (hgcd_r0, a);
232 mpz_init_set (hgcd_r1, b);
233 if (bsize < asize)
234 {
235 _mpz_realloc (hgcd_r1, asize);
236 MPN_ZERO (hgcd_r1->_mp_d + bsize, asize - bsize);
237 }
238 res[1] = mpn_hgcd (hgcd_r0->_mp_d,
239 hgcd_r1->_mp_d,
240 asize,
241 &hgcd, hgcd_tp);
242
243 if (res[0] != res[1])
244 {
245 fprintf (stderr, "ERROR in test %d\n", i);
246 fprintf (stderr, "Different return value from hgcd and hgcd_ref\n");
247 fprintf (stderr, "op1="); debug_mp (a, -16);
248 fprintf (stderr, "op2="); debug_mp (b, -16);
249 fprintf (stderr, "hgcd_ref: %ld\n", (long) res[0]);
250 fprintf (stderr, "mpn_hgcd: %ld\n", (long) res[1]);
251 abort ();
252 }
253 if (res[0] > 0)
254 {
255 if (!hgcd_ref_equal (&hgcd, &ref)
256 || !mpz_mpn_equal (ref_r0, hgcd_r0->_mp_d, res[1])
257 || !mpz_mpn_equal (ref_r1, hgcd_r1->_mp_d, res[1]))
258 {
259 fprintf (stderr, "ERROR in test %d\n", i);
260 fprintf (stderr, "mpn_hgcd and hgcd_ref returned different values\n");
261 fprintf (stderr, "op1="); debug_mp (a, -16);
262 fprintf (stderr, "op2="); debug_mp (b, -16);
263 abort ();
264 }
265 }
266
267 refmpn_free_limbs (hgcd_init_tp);
268 refmpn_free_limbs (hgcd_tp);
269 hgcd_ref_clear (&ref);
270 mpz_clear (ref_r0);
271 mpz_clear (ref_r1);
272 mpz_clear (hgcd_r0);
273 mpz_clear (hgcd_r1);
274
275 return res[0];
276 }
277
278 static void
279 hgcd_ref_init (struct hgcd_ref *hgcd)
280 {
281 unsigned i;
282 for (i = 0; i<2; i++)
283 {
284 unsigned j;
285 for (j = 0; j<2; j++)
286 mpz_init (hgcd->m[i][j]);
287 }
288 }
289
290 static void
291 hgcd_ref_clear (struct hgcd_ref *hgcd)
292 {
293 unsigned i;
294 for (i = 0; i<2; i++)
295 {
296 unsigned j;
297 for (j = 0; j<2; j++)
298 mpz_clear (hgcd->m[i][j]);
299 }
300 }
301
302
303 static int
304 sdiv_qr (mpz_t q, mpz_t r, mp_size_t s, const mpz_t a, const mpz_t b)
305 {
306 mpz_fdiv_qr (q, r, a, b);
307 if (mpz_size (r) <= s)
308 {
309 mpz_add (r, r, b);
310 mpz_sub_ui (q, q, 1);
311 }
312
313 return (mpz_sgn (q) > 0);
314 }
315
316 static int
317 hgcd_ref (struct hgcd_ref *hgcd, mpz_t a, mpz_t b)
318 {
319 mp_size_t n = MAX (mpz_size (a), mpz_size (b));
320 mp_size_t s = n/2 + 1;
321 mp_size_t asize;
322 mp_size_t bsize;
323 mpz_t q;
324 int res;
325
326 if (mpz_size (a) <= s || mpz_size (b) <= s)
327 return 0;
328
329 res = mpz_cmp (a, b);
330 if (res < 0)
331 {
332 mpz_sub (b, b, a);
333 if (mpz_size (b) <= s)
334 return 0;
335
336 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 0);
337 mpz_set_ui (hgcd->m[1][0], 1); mpz_set_ui (hgcd->m[1][1], 1);
338 }
339 else if (res > 0)
340 {
341 mpz_sub (a, a, b);
342 if (mpz_size (a) <= s)
343 return 0;
344
345 mpz_set_ui (hgcd->m[0][0], 1); mpz_set_ui (hgcd->m[0][1], 1);
346 mpz_set_ui (hgcd->m[1][0], 0); mpz_set_ui (hgcd->m[1][1], 1);
347 }
348 else
349 return 0;
350
351 mpz_init (q);
352
353 for (;;)
354 {
355 ASSERT (mpz_size (a) > s);
356 ASSERT (mpz_size (b) > s);
357
358 if (mpz_cmp (a, b) > 0)
359 {
360 if (!sdiv_qr (q, a, s, a, b))
361 break;
362 mpz_addmul (hgcd->m[0][1], q, hgcd->m[0][0]);
363 mpz_addmul (hgcd->m[1][1], q, hgcd->m[1][0]);
364 }
365 else
366 {
367 if (!sdiv_qr (q, b, s, b, a))
368 break;
369 mpz_addmul (hgcd->m[0][0], q, hgcd->m[0][1]);
370 mpz_addmul (hgcd->m[1][0], q, hgcd->m[1][1]);
371 }
372 }
373
374 mpz_clear (q);
375
376 asize = mpz_size (a);
377 bsize = mpz_size (b);
378 return MAX (asize, bsize);
379 }
380
381 static int
382 mpz_mpn_equal (const mpz_t a, mp_srcptr bp, mp_size_t bsize)
383 {
384 mp_srcptr ap = a->_mp_d;
385 mp_size_t asize = a->_mp_size;
386
387 MPN_NORMALIZE (bp, bsize);
388 return asize == bsize && mpn_cmp (ap, bp, asize) == 0;
389 }
390
391 static int
392 hgcd_ref_equal (const struct hgcd_matrix *hgcd, const struct hgcd_ref *ref)
393 {
394 unsigned i;
395
396 for (i = 0; i<2; i++)
397 {
398 unsigned j;
399
400 for (j = 0; j<2; j++)
401 if (!mpz_mpn_equal (ref->m[i][j], hgcd->p[i][j], hgcd->n))
402 return 0;
403 }
404
405 return 1;
406 }