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