1 /*
2     Copyright (C) 2018 Daniel Schultz
3 
4     This file is part of FLINT.
5 
6     FLINT is free software: you can redistribute it and/or modify it under
7     the terms of the GNU Lesser General Public License (LGPL) as published
8     by the Free Software Foundation; either version 2.1 of the License, or
9     (at your option) any later version.  See <http://www.gnu.org/licenses/>.
10 */
11 
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include "fmpz_mpoly.h"
15 
gcd_check(fmpz_mpoly_t g,fmpz_mpoly_t a,fmpz_mpoly_t b,const fmpz_mpoly_t gdiv,fmpz_mpoly_ctx_t ctx,slong i,slong j,const char * name)16 void gcd_check(
17     fmpz_mpoly_t g,
18     fmpz_mpoly_t a,
19     fmpz_mpoly_t b,
20     const fmpz_mpoly_t gdiv,
21     fmpz_mpoly_ctx_t ctx,
22     slong i,
23     slong j,
24     const char * name)
25 {
26     int res;
27     fmpz_mpoly_t ca, cb, cg;
28 
29     fmpz_mpoly_init(ca, ctx);
30     fmpz_mpoly_init(cb, ctx);
31     fmpz_mpoly_init(cg, ctx);
32 
33     res = fmpz_mpoly_gcd(g, a, b, ctx);
34     fmpz_mpoly_assert_canonical(g, ctx);
35 
36     if (!res)
37     {
38         flint_printf("Check gcd can be computed\n"
39                                          "i = %wd, j = %wd, %s\n", i, j, name);
40         flint_abort();
41     }
42 
43     if (!fmpz_mpoly_is_zero(gdiv, ctx))
44     {
45         if (!fmpz_mpoly_divides(ca, g, gdiv, ctx))
46         {
47             printf("FAIL\n");
48             flint_printf("Check divisor of gcd\n"
49                                          "i = %wd, j = %wd, %s\n", i, j, name);
50             flint_abort();
51         }
52     }
53 
54     if (fmpz_mpoly_is_zero(g, ctx))
55     {
56         if (!fmpz_mpoly_is_zero(a, ctx) || !fmpz_mpoly_is_zero(b, ctx))
57         {
58             printf("FAIL\n");
59             flint_printf("Check zero gcd only results from zero inputs\n"
60                                          "i = %wd, j = %wd, %s\n", i, j, name);
61             flint_abort();
62         }
63         goto cleanup;
64     }
65 
66     if (fmpz_sgn(g->coeffs + 0) <= 0)
67     {
68         printf("FAIL\n");
69         flint_printf("Check gcd has positive lc\n"
70                                          "i = %wd, j = %wd, %s\n", i, j, name);
71         flint_abort();
72     }
73 
74     if ((i + j % 11) == 0)
75     {
76         fmpz_mpoly_set(cg, b, ctx);
77         fmpz_mpoly_gcd(cg, cg, a, ctx);
78         if (!fmpz_mpoly_equal(cg, g, ctx))
79         {
80             printf("FAIL\n");
81             flint_printf("Check aliasing 1\n"
82                                          "i = %wd, j = %wd, %s\n", i, j, name);
83             flint_abort();
84         }
85     }
86 
87     if ((i + j % 9) == 0)
88     {
89         fmpz_mpoly_set(cg, b, ctx);
90         fmpz_mpoly_gcd(cg, a, cg, ctx);
91         if (!fmpz_mpoly_equal(cg, g, ctx))
92         {
93             printf("FAIL\n");
94             flint_printf("Check aliasing 2\n"
95                                          "i = %wd, j = %wd, %s\n", i, j, name);
96             flint_abort();
97         }
98     }
99 
100     res = 1;
101     res = res && fmpz_mpoly_divides(ca, a, g, ctx);
102     res = res && fmpz_mpoly_divides(cb, b, g, ctx);
103     if (!res)
104     {
105         printf("FAIL\n");
106         flint_printf("Check divisibility\n"
107                                          "i = %wd, j = %wd, %s\n", i, j, name);
108         flint_abort();
109     }
110 
111     res = fmpz_mpoly_gcd(cg, ca, cb, ctx);
112     fmpz_mpoly_assert_canonical(cg, ctx);
113 
114     if (!res)
115     {
116         printf("FAIL\n");
117         flint_printf("Check gcd of cofactors can be computed\n"
118                                          "i = %wd, j = %wd, %s\n", i, j, name);
119         flint_abort();
120     }
121 
122     if (!fmpz_mpoly_is_one(cg, ctx))
123     {
124         printf("FAIL\n");
125         flint_printf("Check gcd of cofactors is one\n"
126                                          "i = %wd, j = %wd, %s\n", i, j, name);
127         flint_abort();
128     }
129 
130 cleanup:
131 
132     fmpz_mpoly_clear(ca, ctx);
133     fmpz_mpoly_clear(cb, ctx);
134     fmpz_mpoly_clear(cg, ctx);
135 }
136 
137 int
main(void)138 main(void)
139 {
140     const slong max_threads = 5;
141     slong i, j, k, tmul = 10;
142     FLINT_TEST_INIT(state);
143 #ifdef _WIN32
144     tmul = 1;
145 #endif
146 
147     flint_printf("gcd....");
148     fflush(stdout);
149 
150     for (i = 3; i <= 8; i++)
151     {
152         fmpz_mpoly_ctx_t ctx;
153         fmpz_mpoly_t g, a, b, t;
154 
155         fmpz_mpoly_ctx_init(ctx, i, ORD_DEGREVLEX);
156         fmpz_mpoly_init(a, ctx);
157         fmpz_mpoly_init(b, ctx);
158         fmpz_mpoly_init(g, ctx);
159         fmpz_mpoly_init(t, ctx);
160 
161         fmpz_mpoly_one(g, ctx);
162         fmpz_mpoly_one(a, ctx);
163         fmpz_mpoly_one(b, ctx);
164         for (j = 0; j < i; j++)
165         {
166             fmpz_mpoly_gen(t, j, ctx);
167             fmpz_mpoly_add_ui(t, t, 1, ctx);
168             fmpz_mpoly_mul(g, g, t, ctx);
169             fmpz_mpoly_gen(t, j, ctx);
170             fmpz_mpoly_sub_ui(t, t, 2, ctx);
171             fmpz_mpoly_mul(a, a, t, ctx);
172             fmpz_mpoly_gen(t, j, ctx);
173             fmpz_mpoly_add_ui(t, t, 2, ctx);
174             fmpz_mpoly_mul(b, b, t, ctx);
175         }
176         fmpz_mpoly_sub_ui(g, g, 2, ctx);
177         fmpz_mpoly_add_ui(a, a, 2, ctx);
178         fmpz_mpoly_sub_ui(b, b, 2, ctx);
179 
180         fmpz_mpoly_mul(a, a, g, ctx);
181         fmpz_mpoly_mul(b, b, g, ctx);
182         fmpz_mpoly_set(t, g, ctx);
183 
184         gcd_check(g, a, b, t, ctx, i, 0, "dense examples");
185 
186         flint_set_num_threads(n_randint(state, max_threads) + 1);
187 
188         fmpz_mpoly_clear(a, ctx);
189         fmpz_mpoly_clear(b, ctx);
190         fmpz_mpoly_clear(g, ctx);
191         fmpz_mpoly_clear(t, ctx);
192         fmpz_mpoly_ctx_clear(ctx);
193     }
194 
195     {
196         fmpz_mpoly_ctx_t ctx;
197         fmpz_mpoly_t g, a, b, t;
198         const char * vars[] = {"t" ,"x", "y", "z"};
199 
200         fmpz_mpoly_ctx_init(ctx, 4, ORD_LEX);
201         fmpz_mpoly_init(a, ctx);
202         fmpz_mpoly_init(b, ctx);
203         fmpz_mpoly_init(g, ctx);
204         fmpz_mpoly_init(t, ctx);
205 
206         fmpz_mpoly_set_str_pretty(t, "39 - t*x + 39*x^100 - t*x^101 + 39*x^3*y - t*x^4*y - 7*x^2*y^3*z^11 - 7*x^102*y^3*z^11 - 7*x^5*y^4*z^11 + 78*t^15*x^78*y^3*z^13 - 2*t^16*x^79*y^3*z^13 + x^1000*y^3*z^20 + x^1100*y^3*z^20 + x^1003*y^4*z^20 - 14*t^15*x^80*y^6*z^24 + 2*t^15*x^1078*y^6*z^33", vars, ctx);
207         fmpz_mpoly_set_str_pretty(a, "39 - t*x - 7*x^2*y^3*z^11 + x^1000*y^3*z^20", vars, ctx);
208         fmpz_mpoly_set_str_pretty(b, "1 + x^100 + x^3*y + 2*t^15*x^78*y^3*z^13", vars, ctx);
209         fmpz_mpoly_mul(a, a, t, ctx);
210         fmpz_mpoly_mul(b, b, t, ctx);
211 
212         gcd_check(g, a, b, t, ctx, 0, 0, "example");
213 
214         fmpz_mpoly_clear(a, ctx);
215         fmpz_mpoly_clear(b, ctx);
216         fmpz_mpoly_clear(g, ctx);
217         fmpz_mpoly_clear(t, ctx);
218         fmpz_mpoly_ctx_clear(ctx);
219     }
220 
221     {
222         int success;
223         fmpz_mpoly_ctx_t ctx;
224         fmpz_mpoly_t g, a, b;
225         const char * vars[] = {"x" ,"y", "z", "t"};
226 
227         fmpz_mpoly_ctx_init(ctx, 4, ORD_LEX);
228         fmpz_mpoly_init(a, ctx);
229         fmpz_mpoly_init(b, ctx);
230         fmpz_mpoly_init(g, ctx);
231 
232         fmpz_mpoly_set_str_pretty(a, "x^3 + 1", vars, ctx);
233         fmpz_mpoly_set_str_pretty(b, "x^9999999999999999999999 + x^3333333333333333333333 + x", vars, ctx);
234 
235         flint_set_num_threads(n_randint(state, max_threads) + 1);
236         success = fmpz_mpoly_gcd(g, a, b, ctx);
237         if (success)
238         {
239             printf("FAIL\n");
240             flint_printf("Check non-example\n");
241             flint_abort();
242         }
243 
244         fmpz_mpoly_clear(a, ctx);
245         fmpz_mpoly_clear(b, ctx);
246         fmpz_mpoly_clear(g, ctx);
247         fmpz_mpoly_ctx_clear(ctx);
248     }
249 
250     {
251         fmpz_mpoly_ctx_t ctx;
252         fmpz_mpoly_t g, a, b, t;
253         const char * vars[] = {"x" ,"y", "z", "t"};
254 
255         fmpz_mpoly_ctx_init(ctx, 4, ORD_LEX);
256         fmpz_mpoly_init(a, ctx);
257         fmpz_mpoly_init(b, ctx);
258         fmpz_mpoly_init(g, ctx);
259         fmpz_mpoly_init(t, ctx);
260 
261         fmpz_mpoly_set_str_pretty(a, "(1 + x)^1*(2 + y)^1*(1 + z)^2", vars, ctx);
262         fmpz_mpoly_set_str_pretty(b, "(2 + x)^1*(1 + y)^1*(1 - z)^2", vars, ctx);
263         fmpz_mpoly_set_str_pretty(t, "(1 - x)^1*(2 - y)^1*(1 - z)^2", vars, ctx);
264         fmpz_mpoly_mul(a, a, t, ctx);
265         fmpz_mpoly_mul(b, b, t, ctx);
266 
267         flint_set_num_threads(n_randint(state, max_threads) + 1);
268         gcd_check(g, a, b, t, ctx, 0, 0, "total dense example");
269 
270         fmpz_mpoly_clear(a, ctx);
271         fmpz_mpoly_clear(b, ctx);
272         fmpz_mpoly_clear(g, ctx);
273         fmpz_mpoly_clear(t, ctx);
274         fmpz_mpoly_ctx_clear(ctx);
275     }
276 
277     /* The gcd should always work when one input is a monomial */
278     for (i = 0; i < tmul * flint_test_multiplier(); i++)
279     {
280         fmpz_mpoly_ctx_t ctx;
281         fmpz_mpoly_t a, b, g, t;
282         slong len, len1, len2;
283         flint_bitcnt_t coeff_bits, exp_bits, exp_bits1, exp_bits2;
284 
285         fmpz_mpoly_ctx_init_rand(ctx, state, 10);
286 
287         fmpz_mpoly_init(g, ctx);
288         fmpz_mpoly_init(a, ctx);
289         fmpz_mpoly_init(b, ctx);
290         fmpz_mpoly_init(t, ctx);
291 
292         len = n_randint(state, 25);
293         len1 = n_randint(state, 50);
294         len2 = n_randint(state, 50);
295         if (n_randlimb(state) & UWORD(1))
296             len1 = FLINT_MIN(len1, WORD(1));
297         else
298             len2 = FLINT_MIN(len2, WORD(1));
299 
300         exp_bits = n_randint(state, 70) + 2;
301         exp_bits1 = n_randint(state, 100) + 2;
302         exp_bits2 = n_randint(state, 100) + 2;
303 
304         coeff_bits = n_randint(state, 200);
305 
306         for (j = 0; j < 4; j++)
307         {
308             do {
309                 fmpz_mpoly_randtest_bits(t, state, 1, coeff_bits + 1, exp_bits, ctx);
310             } while (t->length != 1);
311             fmpz_mpoly_randtest_bits(a, state, len1, coeff_bits, exp_bits1, ctx);
312             fmpz_mpoly_randtest_bits(b, state, len2, coeff_bits, exp_bits2, ctx);
313             fmpz_mpoly_mul(a, a, t, ctx);
314             fmpz_mpoly_mul(b, b, t, ctx);
315 
316             fmpz_mpoly_randtest_bits(g, state, len, coeff_bits, exp_bits, ctx);
317 
318             flint_set_num_threads(n_randint(state, max_threads) + 1);
319             gcd_check(g, a, b, t, ctx, i, j, "monomial");
320         }
321 
322         fmpz_mpoly_clear(g, ctx);
323         fmpz_mpoly_clear(a, ctx);
324         fmpz_mpoly_clear(b, ctx);
325         fmpz_mpoly_clear(t, ctx);
326         fmpz_mpoly_ctx_clear(ctx);
327     }
328 
329     /* The gcd should always work when both cofactors are monomials */
330     for (i = 0; i < tmul * flint_test_multiplier(); i++)
331     {
332         fmpz_mpoly_ctx_t ctx;
333         fmpz_mpoly_t a, b, g, t1, t2;
334         slong len, len1;
335         flint_bitcnt_t coeff_bits, exp_bits, exp_bits1, exp_bits2;
336 
337         fmpz_mpoly_ctx_init_rand(ctx, state, 10);
338 
339         fmpz_mpoly_init(g, ctx);
340         fmpz_mpoly_init(a, ctx);
341         fmpz_mpoly_init(b, ctx);
342         fmpz_mpoly_init(t1, ctx);
343         fmpz_mpoly_init(t2, ctx);
344 
345         len = n_randint(state, 25);
346         len1 = n_randint(state, 25);
347 
348         exp_bits = n_randint(state, 70) + 2;
349         exp_bits1 = n_randint(state, 100) + 2;
350         exp_bits2 = n_randint(state, 100) + 2;
351 
352         coeff_bits = n_randint(state, 200);
353 
354         for (j = 0; j < 4; j++)
355         {
356             do {
357                 fmpz_mpoly_randtest_bits(t1, state, 1, coeff_bits + 1, exp_bits1, ctx);
358             } while (t1->length != 1);
359             do {
360                 fmpz_mpoly_randtest_bits(t2, state, 1, coeff_bits + 1, exp_bits2, ctx);
361             } while (t2->length != 1);
362             fmpz_mpoly_randtest_bits(a, state, len1, coeff_bits, exp_bits, ctx);
363             fmpz_mpoly_mul(b, a, t1, ctx);
364             fmpz_mpoly_mul(t2, a, t2, ctx);
365 
366             fmpz_mpoly_randtest_bits(g, state, len, coeff_bits, exp_bits, ctx);
367 
368             flint_set_num_threads(n_randint(state, max_threads) + 1);
369             gcd_check(g, t2, b, a, ctx, i, j, "monomial cofactors");
370         }
371 
372         fmpz_mpoly_clear(g, ctx);
373         fmpz_mpoly_clear(a, ctx);
374         fmpz_mpoly_clear(b, ctx);
375         fmpz_mpoly_clear(t1, ctx);
376         fmpz_mpoly_clear(t2, ctx);
377         fmpz_mpoly_ctx_clear(ctx);
378     }
379 
380     /* one input divides the other */
381     for (i = 0; i < tmul * flint_test_multiplier(); i++)
382     {
383         fmpz_t c;
384         fmpz_mpoly_ctx_t ctx;
385         fmpz_mpoly_t a, b, g, t1, t2;
386         slong len, len1, len2;
387         mp_limb_t exp_bound, exp_bound1, exp_bound2;
388         flint_bitcnt_t coeff_bits;
389 
390         fmpz_mpoly_ctx_init_rand(ctx, state, 10);
391 
392         fmpz_init(c);
393         fmpz_mpoly_init(g, ctx);
394         fmpz_mpoly_init(a, ctx);
395         fmpz_mpoly_init(b, ctx);
396         fmpz_mpoly_init(t1, ctx);
397         fmpz_mpoly_init(t2, ctx);
398 
399         len = n_randint(state, 5);
400         len1 = n_randint(state, 5);
401         len2 = n_randint(state, 5);
402 
403         exp_bound = n_randint(state, 100) + 2;
404         exp_bound1 = n_randint(state, 100) + 2;
405         exp_bound2 = n_randint(state, 100) + 2;
406 
407         coeff_bits = n_randint(state, 200);
408 
409         for (j = 0; j < 4; j++)
410         {
411             fmpz_mpoly_randtest_bound(t1, state, len1, coeff_bits + 1, exp_bound1, ctx);
412             fmpz_mpoly_randtest_bound(t2, state, len2, coeff_bits + 1, exp_bound2, ctx);
413             fmpz_mpoly_mul(b, t1, t2, ctx);
414             fmpz_randtest(c, state, coeff_bits + 1);
415             fmpz_mpoly_scalar_mul_fmpz(a, t2, c, ctx);
416             fmpz_randtest(c, state, coeff_bits + 1);
417             fmpz_mpoly_scalar_mul_fmpz(b, b, c, ctx);
418 
419             fmpz_mpoly_randtest_bound(g, state, len, coeff_bits, exp_bound, ctx);
420 
421             flint_set_num_threads(n_randint(state, max_threads) + 1);
422 
423             if ((j%2) == 0)
424                 fmpz_mpoly_swap(a, b, ctx);
425 
426             gcd_check(g, a, b, t2, ctx, i, j, "one input divides the other");
427         }
428 
429         fmpz_clear(c);
430         fmpz_mpoly_clear(g, ctx);
431         fmpz_mpoly_clear(a, ctx);
432         fmpz_mpoly_clear(b, ctx);
433         fmpz_mpoly_clear(t1, ctx);
434         fmpz_mpoly_clear(t2, ctx);
435         fmpz_mpoly_ctx_clear(ctx);
436     }
437 
438     /* sparse inputs */
439     for (i = 0; i < tmul * flint_test_multiplier(); i++)
440     {
441         fmpz_mpoly_ctx_t ctx;
442         fmpz_mpoly_t a, b, g, t;
443         flint_bitcnt_t coeff_bits;
444         slong len, len1, len2;
445         slong degbound;
446 
447         fmpz_mpoly_ctx_init_rand(ctx, state, 5);
448 
449         fmpz_mpoly_init(g, ctx);
450         fmpz_mpoly_init(a, ctx);
451         fmpz_mpoly_init(b, ctx);
452         fmpz_mpoly_init(t, ctx);
453 
454         len = n_randint(state, 20) + 1;
455         len1 = n_randint(state, 30);
456         len2 = n_randint(state, 30);
457 
458         degbound = 30/(2*ctx->minfo->nvars - 1);
459 
460         coeff_bits = n_randint(state, 200);
461 
462         for (j = 0; j < 4; j++)
463         {
464             do {
465                 fmpz_mpoly_randtest_bound(t, state, len, coeff_bits + 1, degbound, ctx);
466             } while (t->length == 0);
467             fmpz_mpoly_randtest_bound(a, state, len1, coeff_bits, degbound, ctx);
468             fmpz_mpoly_randtest_bound(b, state, len2, coeff_bits, degbound, ctx);
469             fmpz_mpoly_mul(a, a, t, ctx);
470             fmpz_mpoly_mul(b, b, t, ctx);
471 
472             fmpz_mpoly_randtest_bits(g, state, len, coeff_bits, FLINT_BITS, ctx);
473 
474             flint_set_num_threads(n_randint(state, max_threads) + 1);
475             gcd_check(g, a, b, t, ctx, i, j, "sparse inputs");
476         }
477 
478         fmpz_mpoly_clear(g, ctx);
479         fmpz_mpoly_clear(a, ctx);
480         fmpz_mpoly_clear(b, ctx);
481         fmpz_mpoly_clear(t, ctx);
482         fmpz_mpoly_ctx_clear(ctx);
483     }
484 
485     /* sparse inputs with random repackings */
486     for (i = 0; i < tmul * flint_test_multiplier(); i++)
487     {
488         fmpz_mpoly_ctx_t ctx;
489         fmpz_mpoly_t a, b, g, t;
490         mp_limb_t rlimb;
491         flint_bitcnt_t coeff_bits, newbits;
492         slong len, len1, len2;
493         slong degbound;
494 
495         fmpz_mpoly_ctx_init_rand(ctx, state, 5);
496 
497         fmpz_mpoly_init(g, ctx);
498         fmpz_mpoly_init(a, ctx);
499         fmpz_mpoly_init(b, ctx);
500         fmpz_mpoly_init(t, ctx);
501 
502         len = n_randint(state, 20) + 1;
503         len1 = n_randint(state, 30);
504         len2 = n_randint(state, 30);
505 
506         degbound = 30/(2*ctx->minfo->nvars - 1);
507 
508         coeff_bits = n_randint(state, 200);
509 
510         for (j = 0; j < 4; j++)
511         {
512             do {
513                 fmpz_mpoly_randtest_bound(t, state, len, coeff_bits + 1, degbound, ctx);
514             } while (t->length == 0);
515             fmpz_mpoly_randtest_bound(a, state, len1, coeff_bits, degbound, ctx);
516             fmpz_mpoly_randtest_bound(b, state, len2, coeff_bits, degbound, ctx);
517             fmpz_mpoly_mul(a, a, t, ctx);
518             fmpz_mpoly_mul(b, b, t, ctx);
519 
520             rlimb = n_randlimb(state);
521 
522             if (rlimb & UWORD(3))
523             {
524                 newbits = a->bits + n_randint(state, 2*FLINT_BITS);
525                 newbits = mpoly_fix_bits(newbits, ctx->minfo);
526                 fmpz_mpoly_repack_bits(a, a, newbits, ctx);
527             }
528 
529             if (rlimb & UWORD(12))
530             {
531                 newbits = b->bits + n_randint(state, 2*FLINT_BITS);
532                 newbits = mpoly_fix_bits(newbits, ctx->minfo);
533                 fmpz_mpoly_repack_bits(b, b, newbits, ctx);
534             }
535 
536             fmpz_mpoly_randtest_bits(g, state, len, coeff_bits, FLINT_BITS, ctx);
537 
538             flint_set_num_threads(n_randint(state, max_threads) + 1);
539             gcd_check(g, a, b, t, ctx, i, j, "sparse input with repacking");
540         }
541 
542         fmpz_mpoly_clear(g, ctx);
543         fmpz_mpoly_clear(a, ctx);
544         fmpz_mpoly_clear(b, ctx);
545         fmpz_mpoly_clear(t, ctx);
546         fmpz_mpoly_ctx_clear(ctx);
547     }
548 
549     /* sparse inputs with random inflations */
550     for (i = 0; i < tmul * flint_test_multiplier(); i++)
551     {
552         fmpz_mpoly_ctx_t ctx;
553         fmpz_mpoly_t a, b, g, t;
554         flint_bitcnt_t coeff_bits;
555         fmpz * shifts1, * shifts2, * strides;
556         flint_bitcnt_t stride_bits, shift_bits;
557         slong len, len1, len2;
558         slong degbound;
559 
560         fmpz_mpoly_ctx_init_rand(ctx, state, 5);
561 
562         fmpz_mpoly_init(g, ctx);
563         fmpz_mpoly_init(a, ctx);
564         fmpz_mpoly_init(b, ctx);
565         fmpz_mpoly_init(t, ctx);
566 
567         len = n_randint(state, 20) + 1;
568         len1 = n_randint(state, 30);
569         len2 = n_randint(state, 30);
570 
571         degbound = 30/(2*ctx->minfo->nvars - 1);
572 
573         coeff_bits = n_randint(state, 200);
574 
575         stride_bits = n_randint(state, 100) + 2;
576         shift_bits = n_randint(state, 100) + 2;
577 
578         shifts1 = flint_malloc(ctx->minfo->nvars*sizeof(fmpz));
579         shifts2 = flint_malloc(ctx->minfo->nvars*sizeof(fmpz));
580         strides = flint_malloc(ctx->minfo->nvars*sizeof(fmpz));
581         for (k = 0; k < ctx->minfo->nvars; k++)
582         {
583             fmpz_init(shifts1 + k);
584             fmpz_init(shifts2 + k);
585             fmpz_init(strides + k);
586         }
587 
588         for (j = 0; j < 4; j++)
589         {
590             do {
591                 fmpz_mpoly_randtest_bound(t, state, len, coeff_bits + 1, degbound, ctx);
592             } while (t->length == 0);
593             fmpz_mpoly_randtest_bound(a, state, len1, coeff_bits, degbound, ctx);
594             fmpz_mpoly_randtest_bound(b, state, len2, coeff_bits, degbound, ctx);
595             fmpz_mpoly_mul(a, a, t, ctx);
596             fmpz_mpoly_mul(b, b, t, ctx);
597 
598             fmpz_mpoly_randtest_bits(g, state, len, coeff_bits, FLINT_BITS, ctx);
599 
600             for (k = 0; k < ctx->minfo->nvars; k++)
601             {
602                 fmpz_randtest_unsigned(shifts1 + k, state, shift_bits);
603                 fmpz_randtest_unsigned(shifts2 + k, state, shift_bits);
604                 fmpz_randtest_unsigned(strides + k, state, stride_bits);
605             }
606             fmpz_mpoly_inflate(a, a, shifts1, strides, ctx);
607             fmpz_mpoly_inflate(b, b, shifts2, strides, ctx);
608 
609             for (k = 0; k < ctx->minfo->nvars; k++)
610             {
611                 if (fmpz_cmp(shifts1 + k, shifts2 + k) > 0)
612                     fmpz_set(shifts1 + k, shifts2 + k);
613             }
614             fmpz_mpoly_inflate(t, t, shifts1, strides, ctx);
615 
616             flint_set_num_threads(n_randint(state, max_threads) + 1);
617             gcd_check(g, a, b, t, ctx, i, j, "sparse input with inflation");
618         }
619 
620         for (k = 0; k < ctx->minfo->nvars; k++)
621         {
622             fmpz_clear(shifts1 + k);
623             fmpz_clear(shifts2 + k);
624             fmpz_clear(strides + k);
625         }
626         flint_free(shifts1);
627         flint_free(shifts2);
628         flint_free(strides);
629 
630         fmpz_mpoly_clear(g, ctx);
631         fmpz_mpoly_clear(a, ctx);
632         fmpz_mpoly_clear(b, ctx);
633         fmpz_mpoly_clear(t, ctx);
634         fmpz_mpoly_ctx_clear(ctx);
635     }
636 
637     /* dense inputs */
638     for (i = 0; i < tmul * flint_test_multiplier(); i++)
639     {
640         fmpz_mpoly_ctx_t ctx;
641         fmpz_mpoly_t a, b, g, t;
642         flint_bitcnt_t coeff_bits1, coeff_bits2, coeff_bits3, coeff_bits4;
643         slong len1, len2, len3, len4;
644         ulong degbounds1[4];
645         ulong degbounds2[4];
646         ulong degbounds3[4];
647         flint_bitcnt_t bits4;
648 
649         fmpz_mpoly_ctx_init_rand(ctx, state, 4);
650 
651         fmpz_mpoly_init(g, ctx);
652         fmpz_mpoly_init(a, ctx);
653         fmpz_mpoly_init(b, ctx);
654         fmpz_mpoly_init(t, ctx);
655 
656         len1 = n_randint(state, 300) + 1;
657         len2 = n_randint(state, 300);
658         len3 = n_randint(state, 300);
659         len4 = n_randint(state, 300);
660 
661         for (j = 0; j < ctx->minfo->nvars; j++)
662         {
663             degbounds1[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
664             degbounds2[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
665             degbounds3[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
666         }
667 
668         bits4 = n_randint(state, 200);
669         coeff_bits1 = n_randint(state, 200);
670         coeff_bits2 = n_randint(state, 200);
671         coeff_bits3 = n_randint(state, 200);
672         coeff_bits4 = n_randint(state, 200);
673 
674         for (j = 0; j < 4; j++)
675         {
676             do {
677                 fmpz_mpoly_randtest_bounds(t, state, len1, coeff_bits1 + 1, degbounds1, ctx);
678             } while (t->length == 0);
679             fmpz_mpoly_randtest_bounds(a, state, len2, coeff_bits2, degbounds2, ctx);
680             fmpz_mpoly_randtest_bounds(b, state, len3, coeff_bits3, degbounds3, ctx);
681             fmpz_mpoly_mul(a, a, t, ctx);
682             fmpz_mpoly_mul(b, b, t, ctx);
683 
684             fmpz_mpoly_randtest_bits(g, state, len4, coeff_bits4, bits4, ctx);
685 
686             flint_set_num_threads(n_randint(state, max_threads) + 1);
687             gcd_check(g, a, b, t, ctx, i, j, "dense input");
688         }
689 
690         fmpz_mpoly_clear(g, ctx);
691         fmpz_mpoly_clear(a, ctx);
692         fmpz_mpoly_clear(b, ctx);
693         fmpz_mpoly_clear(t, ctx);
694         fmpz_mpoly_ctx_clear(ctx);
695     }
696 
697     /* dense inputs with repacking */
698     for (i = 0; i < tmul * flint_test_multiplier(); i++)
699     {
700         fmpz_mpoly_ctx_t ctx;
701         fmpz_mpoly_t a, b, g, t;
702         mp_limb_t rlimb;
703         flint_bitcnt_t newbits;
704         flint_bitcnt_t coeff_bits1, coeff_bits2, coeff_bits3, coeff_bits4;
705         slong len1, len2, len3, len4;
706         ulong degbounds1[4];
707         ulong degbounds2[4];
708         ulong degbounds3[4];
709         flint_bitcnt_t bits4;
710 
711         fmpz_mpoly_ctx_init_rand(ctx, state, 4);
712 
713         fmpz_mpoly_init(g, ctx);
714         fmpz_mpoly_init(a, ctx);
715         fmpz_mpoly_init(b, ctx);
716         fmpz_mpoly_init(t, ctx);
717 
718         len1 = n_randint(state, 300) + 1;
719         len2 = n_randint(state, 300);
720         len3 = n_randint(state, 300);
721         len4 = n_randint(state, 300);
722 
723         for (j = 0; j < ctx->minfo->nvars; j++)
724         {
725             degbounds1[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
726             degbounds2[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
727             degbounds3[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
728         }
729 
730         bits4 = n_randint(state, 200);
731         coeff_bits1 = n_randint(state, 200);
732         coeff_bits2 = n_randint(state, 200);
733         coeff_bits3 = n_randint(state, 200);
734         coeff_bits4 = n_randint(state, 200);
735 
736         for (j = 0; j < 4; j++)
737         {
738             do {
739                 fmpz_mpoly_randtest_bounds(t, state, len1, coeff_bits1 + 1, degbounds1, ctx);
740             } while (t->length == 0);
741             fmpz_mpoly_randtest_bounds(a, state, len2, coeff_bits2, degbounds2, ctx);
742             fmpz_mpoly_randtest_bounds(b, state, len3, coeff_bits3, degbounds3, ctx);
743             fmpz_mpoly_mul(a, a, t, ctx);
744             fmpz_mpoly_mul(b, b, t, ctx);
745 
746             rlimb = n_randlimb(state);
747 
748             if (rlimb & UWORD(3))
749             {
750                 newbits = a->bits + n_randint(state, 2*FLINT_BITS);
751                 newbits = mpoly_fix_bits(newbits, ctx->minfo);
752                 fmpz_mpoly_repack_bits(a, a, newbits, ctx);
753             }
754 
755             if (rlimb & UWORD(12))
756             {
757                 newbits = b->bits + n_randint(state, 2*FLINT_BITS);
758                 newbits = mpoly_fix_bits(newbits, ctx->minfo);
759                 fmpz_mpoly_repack_bits(b, b, newbits, ctx);
760             }
761 
762             fmpz_mpoly_randtest_bits(g, state, len4, coeff_bits4, bits4, ctx);
763 
764             flint_set_num_threads(n_randint(state, max_threads) + 1);
765             gcd_check(g, a, b, t, ctx, i, j, "dense input with repacking");
766         }
767 
768         fmpz_mpoly_clear(g, ctx);
769         fmpz_mpoly_clear(a, ctx);
770         fmpz_mpoly_clear(b, ctx);
771         fmpz_mpoly_clear(t, ctx);
772         fmpz_mpoly_ctx_clear(ctx);
773     }
774 
775     /* dense inputs with random inflations */
776     for (i = 0; i < tmul * flint_test_multiplier(); i++)
777     {
778         fmpz_mpoly_ctx_t ctx;
779         fmpz_mpoly_t a, b, g, t;
780         fmpz * shifts1, * shifts2, * strides;
781         flint_bitcnt_t stride_bits, shift_bits;
782         flint_bitcnt_t coeff_bits1, coeff_bits2, coeff_bits3, coeff_bits4;
783         slong len1, len2, len3, len4;
784         ulong degbounds1[4];
785         ulong degbounds2[4];
786         ulong degbounds3[4];
787         flint_bitcnt_t bits4;
788 
789         fmpz_mpoly_ctx_init_rand(ctx, state, 4);
790 
791         fmpz_mpoly_init(g, ctx);
792         fmpz_mpoly_init(a, ctx);
793         fmpz_mpoly_init(b, ctx);
794         fmpz_mpoly_init(t, ctx);
795 
796         len1 = n_randint(state, 300) + 1;
797         len2 = n_randint(state, 300);
798         len3 = n_randint(state, 300);
799         len4 = n_randint(state, 300);
800 
801         for (j = 0; j < ctx->minfo->nvars; j++)
802         {
803             degbounds1[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
804             degbounds2[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
805             degbounds3[j] = 1 + n_randint(state, 15/ctx->minfo->nvars);
806         }
807 
808         bits4 = n_randint(state, 200);
809         coeff_bits1 = n_randint(state, 200);
810         coeff_bits2 = n_randint(state, 200);
811         coeff_bits3 = n_randint(state, 200);
812         coeff_bits4 = n_randint(state, 200);
813 
814         stride_bits = n_randint(state, 100) + 2;
815         shift_bits = n_randint(state, 100) + 2;
816 
817         shifts1 = flint_malloc(ctx->minfo->nvars*sizeof(fmpz));
818         shifts2 = flint_malloc(ctx->minfo->nvars*sizeof(fmpz));
819         strides = flint_malloc(ctx->minfo->nvars*sizeof(fmpz));
820         for (k = 0; k < ctx->minfo->nvars; k++)
821         {
822             fmpz_init(shifts1 + k);
823             fmpz_init(shifts2 + k);
824             fmpz_init(strides + k);
825         }
826 
827         for (j = 0; j < 4; j++)
828         {
829             do {
830                 fmpz_mpoly_randtest_bounds(t, state, len1, coeff_bits1 + 1, degbounds1, ctx);
831             } while (t->length == 0);
832             fmpz_mpoly_randtest_bounds(a, state, len2, coeff_bits2, degbounds2, ctx);
833             fmpz_mpoly_randtest_bounds(b, state, len3, coeff_bits3, degbounds3, ctx);
834             fmpz_mpoly_mul(a, a, t, ctx);
835             fmpz_mpoly_mul(b, b, t, ctx);
836 
837             fmpz_mpoly_randtest_bits(g, state, len4, coeff_bits4, bits4, ctx);
838 
839             for (k = 0; k < ctx->minfo->nvars; k++)
840             {
841                 fmpz_randtest_unsigned(shifts1 + k, state, shift_bits);
842                 fmpz_randtest_unsigned(shifts2 + k, state, shift_bits);
843                 fmpz_randtest_unsigned(strides + k, state, stride_bits);
844             }
845             fmpz_mpoly_inflate(a, a, shifts1, strides, ctx);
846             fmpz_mpoly_inflate(b, b, shifts2, strides, ctx);
847 
848             for (k = 0; k < ctx->minfo->nvars; k++)
849             {
850                 if (fmpz_cmp(shifts1 + k, shifts2 + k) > 0)
851                     fmpz_set(shifts1 + k, shifts2 + k);
852             }
853             fmpz_mpoly_inflate(t, t, shifts1, strides, ctx);
854 
855             flint_set_num_threads(n_randint(state, max_threads) + 1);
856             gcd_check(g, a, b, t, ctx, i, j, "dense input with inflation");
857         }
858 
859         for (k = 0; k < ctx->minfo->nvars; k++)
860         {
861             fmpz_clear(shifts1 + k);
862             fmpz_clear(shifts2 + k);
863             fmpz_clear(strides + k);
864         }
865         flint_free(shifts1);
866         flint_free(shifts2);
867         flint_free(strides);
868 
869         fmpz_mpoly_clear(g, ctx);
870         fmpz_mpoly_clear(a, ctx);
871         fmpz_mpoly_clear(b, ctx);
872         fmpz_mpoly_clear(t, ctx);
873         fmpz_mpoly_ctx_clear(ctx);
874     }
875 
876     FLINT_TEST_CLEANUP(state);
877 
878     printf("PASS\n");
879     return 0;
880 }
881 
882