1 /*
2 Copyright (C) 2018-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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include "mpoly.h"
13
14
mpoly_gcd_info_init(mpoly_gcd_info_t I,slong nvars)15 void mpoly_gcd_info_init(mpoly_gcd_info_t I, slong nvars)
16 {
17 char * d;
18
19 FLINT_ASSERT(nvars > 0);
20
21 d = (char *) flint_malloc(nvars*(10*sizeof(ulong) + 12*sizeof(slong)));
22
23 I->data = d;
24
25 I->Amax_exp = (ulong *) d; d += nvars*sizeof(ulong);
26 I->Amin_exp = (ulong *) d; d += nvars*sizeof(ulong);
27 I->Astride = (ulong *) d; d += nvars*sizeof(ulong);
28 I->Adeflate_deg = (slong *) d; d += nvars*sizeof(slong);
29 I->Alead_count = (slong *) d; d += nvars*sizeof(slong);
30 I->Atail_count = (slong *) d; d += nvars*sizeof(slong);
31
32 I->Bmax_exp = (ulong *) d; d += nvars*sizeof(ulong);
33 I->Bmin_exp = (ulong *) d; d += nvars*sizeof(ulong);
34 I->Bstride = (ulong *) d; d += nvars*sizeof(ulong);
35 I->Bdeflate_deg = (slong *) d; d += nvars*sizeof(slong);
36 I->Blead_count = (slong *) d; d += nvars*sizeof(slong);
37 I->Btail_count = (slong *) d; d += nvars*sizeof(slong);
38
39 I->Gmin_exp = (ulong *) d; d += nvars*sizeof(ulong);
40 I->Abarmin_exp = (ulong *) d; d += nvars*sizeof(ulong);
41 I->Bbarmin_exp = (ulong *) d; d += nvars*sizeof(ulong);
42 I->Gstride = (ulong *) d; d += nvars*sizeof(ulong);
43 I->Gdeflate_deg_bound = (slong *) d; d += nvars*sizeof(slong);
44 I->Gterm_count_est = (slong *) d; d += nvars*sizeof(slong);
45
46 I->hensel_perm = (slong *) d; d += nvars*sizeof(slong);
47 I->brown_perm = (slong *) d; d += nvars*sizeof(slong);
48 I->zippel_perm = (slong *) d; d += nvars*sizeof(slong);
49 I->zippel2_perm = (slong *) d; d += nvars*sizeof(slong);
50 }
51
mpoly_gcd_info_clear(mpoly_gcd_info_t I)52 void mpoly_gcd_info_clear(mpoly_gcd_info_t I)
53 {
54 flint_free(I->data);
55 }
56
57 /*
58 Scan A and fill in the min and max exponents of each variable along
59 with the count of terms attached to each.
60 */
mpoly_gcd_info_limits(ulong * Amax_exp,ulong * Amin_exp,slong * Amax_exp_count,slong * Amin_exp_count,const ulong * Aexps,flint_bitcnt_t Abits,slong Alength,const mpoly_ctx_t mctx)61 void mpoly_gcd_info_limits(ulong * Amax_exp, ulong * Amin_exp,
62 slong * Amax_exp_count, slong * Amin_exp_count,
63 const ulong * Aexps, flint_bitcnt_t Abits, slong Alength,
64 const mpoly_ctx_t mctx)
65 {
66 ulong * exps;
67 slong i, j, N;
68 slong nvars = mctx->nvars;
69 TMP_INIT;
70
71 FLINT_ASSERT(Alength > 0);
72 FLINT_ASSERT(Abits <= FLINT_BITS);
73
74 TMP_START;
75
76 exps = (ulong *) TMP_ALLOC(nvars*sizeof(ulong));
77
78 N = mpoly_words_per_exp(Abits, mctx);
79
80 i = 0;
81 mpoly_get_monomial_ui(exps, Aexps + N*i, Abits, mctx);
82 for (j = 0; j < nvars; j++)
83 {
84 Amin_exp[j] = exps[j];
85 Amax_exp[j] = exps[j];
86 Amin_exp_count[j] = 1;
87 Amax_exp_count[j] = 1;
88 }
89 for (i = 1; i < Alength; i++)
90 {
91 mpoly_get_monomial_ui(exps, Aexps + N*i, Abits, mctx);
92
93 for (j = 0; j < nvars; j++)
94 {
95 if (Amin_exp[j] > exps[j])
96 {
97 Amin_exp[j] = exps[j];
98 Amin_exp_count[j] = 1;
99 }
100 else if (Amin_exp[j] == exps[j])
101 {
102 Amin_exp_count[j] += 1;
103 }
104
105 if (Amax_exp[j] < exps[j])
106 {
107 Amax_exp[j] = exps[j];
108 Amax_exp_count[j] = 1;
109 }
110 else if (Amax_exp[j] == exps[j])
111 {
112 Amax_exp_count[j] += 1;
113 }
114 }
115 }
116
117 TMP_END;
118 }
119
120
121 /*
122 For each variable v, let SA[v] be the set of exponents of variable v in A.
123 Ditto for SB[v]. The function computes
124 strides[v] = GCD(SA[v] - min(SA[v]), SB[v] - min(SB[v]))
125 It is assumed that {A|B}{max|min}_exp are correct.
126 */
mpoly_gcd_info_stride(ulong * strides,const ulong * Aexps,flint_bitcnt_t Abits,slong Alength,const ulong * Amax_exp,const ulong * Amin_exp,const ulong * Bexps,flint_bitcnt_t Bbits,slong Blength,const ulong * Bmax_exp,const ulong * Bmin_exp,const mpoly_ctx_t mctx)127 void mpoly_gcd_info_stride(ulong * strides,
128 const ulong * Aexps, flint_bitcnt_t Abits, slong Alength,
129 const ulong * Amax_exp, const ulong * Amin_exp,
130 const ulong * Bexps, flint_bitcnt_t Bbits, slong Blength,
131 const ulong * Bmax_exp, const ulong * Bmin_exp,
132 const mpoly_ctx_t mctx)
133 {
134 slong i, j, NA, NB;
135 slong nvars = mctx->nvars;
136 ulong mask;
137 ulong * exps;
138 TMP_INIT;
139
140 FLINT_ASSERT(Abits <= FLINT_BITS);
141 FLINT_ASSERT(Bbits <= FLINT_BITS);
142
143 for (j = 0; j < nvars; j++)
144 {
145 strides[j] = n_gcd(Amax_exp[j] - Amin_exp[j],
146 Bmax_exp[j] - Bmin_exp[j]);
147 }
148
149 TMP_START;
150 exps = (ulong *) TMP_ALLOC(nvars*sizeof(ulong));
151
152 NA = mpoly_words_per_exp(Abits, mctx);
153
154 for (i = 0; i < Alength; i++)
155 {
156 mpoly_get_monomial_ui(exps, Aexps + NA*i, Abits, mctx);
157 mask = 0;
158 for (j = 0; j < nvars; j++)
159 {
160 strides[j] = n_gcd(strides[j], exps[j] - Amin_exp[j]);
161 mask |= strides[j];
162 }
163 if (mask < UWORD(2))
164 {
165 goto cleanup;
166 }
167 }
168
169 NB = mpoly_words_per_exp(Bbits, mctx);
170
171 for (i = 0; i < Blength; i++)
172 {
173 mpoly_get_monomial_ui(exps, Bexps + NB*i, Bbits, mctx);
174 mask = 0;
175 for (j = 0; j < nvars; j++)
176 {
177 strides[j] = n_gcd(strides[j], exps[j] - Bmin_exp[j]);
178 mask |= strides[j];
179 }
180 if (mask < UWORD(2))
181 {
182 goto cleanup;
183 }
184 }
185
186 cleanup:
187
188 TMP_END;
189
190 return;
191 }
192
mpoly_gcd_info_set_perm(mpoly_gcd_info_t I,slong Alength,slong Blength,const mpoly_ctx_t mctx)193 void mpoly_gcd_info_set_perm(
194 mpoly_gcd_info_t I,
195 slong Alength,
196 slong Blength,
197 const mpoly_ctx_t mctx)
198 {
199 slong j, m;
200
201 I->Adensity = Alength;
202 I->Bdensity = Blength;
203
204 m = 0;
205 for (j = 0; j < mctx->nvars; j++)
206 {
207 if (I->Amax_exp[j] > I->Amin_exp[j])
208 {
209 FLINT_ASSERT(I->Gstride[j] != UWORD(0));
210 FLINT_ASSERT((I->Amax_exp[j] - I->Amin_exp[j]) % I->Gstride[j] == 0);
211 FLINT_ASSERT((I->Bmax_exp[j] - I->Bmin_exp[j]) % I->Gstride[j] == 0);
212
213 I->Adensity /= UWORD(1) + (ulong)(I->Adeflate_deg[j]);
214 I->Bdensity /= UWORD(1) + (ulong)(I->Bdeflate_deg[j]);
215
216 I->hensel_perm[m] = j;
217 I->brown_perm[m] = j;
218 I->zippel_perm[m] = j;
219 I->zippel2_perm[m] = j;
220 m++;
221 }
222 else
223 {
224 FLINT_ASSERT(I->Amax_exp[j] == I->Amin_exp[j]);
225 FLINT_ASSERT(I->Bmax_exp[j] == I->Bmin_exp[j]);
226 }
227 }
228
229 I->mvars = m;
230
231 I->can_use = 0;
232 }
233
234
mpoly_gcd_info_measure_hensel(mpoly_gcd_info_t I,slong Alength,slong Blength,const mpoly_ctx_t mctx)235 void mpoly_gcd_info_measure_hensel(
236 mpoly_gcd_info_t I,
237 slong Alength,
238 slong Blength,
239 const mpoly_ctx_t mctx)
240 {
241 slong i, k;
242 slong m = I->mvars;
243 slong * perm = I->hensel_perm;
244 flint_bitcnt_t abits, bbits;
245 double te, tg, ta, tb;
246 double stgab, mtgab, iblend, eblend;
247
248 /* need at least 2 variables */
249 if (m < 2)
250 return;
251
252 abits = FLINT_BIT_COUNT(Alength);
253 bbits = FLINT_BIT_COUNT(Blength);
254
255 te = tg = ta = tb = 1;
256 for (i = 0; i < m; i++)
257 {
258 double x;
259
260 k = perm[i];
261 if (abits + FLINT_BIT_COUNT(I->Adeflate_deg[k]) > FLINT_BITS ||
262 bbits + FLINT_BIT_COUNT(I->Bdeflate_deg[k]) > FLINT_BITS)
263 {
264 /* each variable is eventually converted to dense storage */
265 return;
266 }
267
268 te *= 1 + FLINT_MAX(I->Adeflate_deg[k], I->Bdeflate_deg[k]);
269
270 x = I->Gdeflate_deg_bound[k];
271 tg *= 1 + x + 0.005*x*x;
272
273 x = FLINT_MAX(0, I->Adeflate_deg[k] - I->Gdeflate_deg_bound[k]);
274 ta *= 1 + x + 0.005*x*x;
275
276 x = FLINT_MAX(0, I->Bdeflate_deg[k] - I->Gdeflate_deg_bound[k]);
277 tb *= 1 + x + 0.005*x*x;
278 }
279
280 iblend = 1;
281 eblend = 1;
282
283 stgab = tg + ta + tb;
284 mtgab = FLINT_MIN(tg, ta);
285 mtgab = FLINT_MIN(mtgab, tb);
286
287 I->can_use |= MPOLY_GCD_USE_HENSEL;
288 I->hensel_time = 0.005*te*(I->Adensity + I->Bdensity)*eblend +
289 0.004*(iblend*stgab + (1 - iblend)*mtgab);
290 }
291
292 /*
293 limit past which we should not test divisibility
294 */
mpoly_gcd_info_get_brown_upper_limit(const mpoly_gcd_info_t I,slong var,slong bound)295 slong mpoly_gcd_info_get_brown_upper_limit(
296 const mpoly_gcd_info_t I,
297 slong var,
298 slong bound)
299 {
300 if (I == NULL || !I->Gdeflate_deg_bounds_are_nice)
301 {
302 return 0;
303 }
304 else
305 {
306 slong k, max;
307 slong density;
308 slong limit;
309
310 FLINT_ASSERT(var < I->mvars);
311 k = I->brown_perm[var];
312 max = FLINT_MAX(I->Adeflate_deg[k], I->Bdeflate_deg[k]);
313 bound = FLINT_MAX(bound, 1 + max);
314 density = 0.5*(I->Adensity + I->Bdensity);
315 limit = bound*((1.125 - density)*(1.125 - density))*0.375;
316 return FLINT_MIN(limit, bound/2);
317 }
318 }
319
mpoly_gcd_info_measure_brown(mpoly_gcd_info_t I,slong Alength,slong Blength,const mpoly_ctx_t mctx)320 void mpoly_gcd_info_measure_brown(
321 mpoly_gcd_info_t I,
322 slong Alength,
323 slong Blength,
324 const mpoly_ctx_t mctx)
325 {
326 slong i, k;
327 slong m = I->mvars;
328 slong * perm = I->brown_perm;
329 flint_bitcnt_t abits, bbits;
330 double te, tg, ta, tb;
331 double stgab, mtgab, iblend, eblend;
332
333 /* need at least 2 variables */
334 if (m < 2)
335 return;
336
337 abits = FLINT_BIT_COUNT(Alength);
338 bbits = FLINT_BIT_COUNT(Blength);
339
340 te = tg = ta = tb = 1;
341 for (i = 0; i < m; i++)
342 {
343 double x;
344
345 k = perm[i];
346 if (abits + FLINT_BIT_COUNT(I->Adeflate_deg[k]) > FLINT_BITS ||
347 bbits + FLINT_BIT_COUNT(I->Bdeflate_deg[k]) > FLINT_BITS)
348 {
349 /* each variable is eventually converted to dense storage */
350 return;
351 }
352
353 te *= 1 + FLINT_MAX(I->Adeflate_deg[k], I->Bdeflate_deg[k]);
354
355 x = I->Gdeflate_deg_bound[k];
356 tg *= 1 + x + 0.005*x*x;
357
358 x = FLINT_MAX(0, I->Adeflate_deg[k] - I->Gdeflate_deg_bound[k]);
359 ta *= 1 + x + 0.005*x*x;
360
361 x = FLINT_MAX(0, I->Bdeflate_deg[k] - I->Gdeflate_deg_bound[k]);
362 tb *= 1 + x + 0.005*x*x;
363 }
364
365 iblend = 1;
366 eblend = 1;
367 if (I->Gdeflate_deg_bounds_are_nice)
368 {
369 slong k = perm[m - 1];
370 slong limit = mpoly_gcd_info_get_brown_upper_limit(I, m - 1, 0);
371 slong expected_stab;
372
373 expected_stab = FLINT_MIN(I->Adeflate_deg[k], I->Bdeflate_deg[k]);
374 expected_stab = expected_stab - I->Gdeflate_deg_bound[k];
375 expected_stab = FLINT_MIN(expected_stab, I->Gdeflate_deg_bound[k]);
376
377 if (expected_stab < limit)
378 {
379 slong bound = 1 + FLINT_MAX(I->Adeflate_deg[k], I->Bdeflate_deg[k]);
380 iblend = (I->Adensity + I->Bdensity);
381 iblend = FLINT_MIN(iblend, 1);
382 iblend = FLINT_MAX(iblend, 0.01);
383 eblend = 0.25 + 0.75*(double)(expected_stab)/(double)(bound);
384 }
385 }
386
387 stgab = tg + ta + tb;
388 mtgab = FLINT_MIN(tg, ta);
389 mtgab = FLINT_MIN(mtgab, tb);
390
391 I->can_use |= MPOLY_GCD_USE_BROWN;
392 I->brown_time = 0.005*te*(I->Adensity + I->Bdensity)*eblend +
393 0.004*(iblend*stgab + (1 - iblend)*mtgab);
394 }
395
mpoly_gcd_info_measure_bma(mpoly_gcd_info_t I,slong Alength,slong Blength,const mpoly_ctx_t mctx)396 void mpoly_gcd_info_measure_bma(
397 mpoly_gcd_info_t I,
398 slong Alength,
399 slong Blength,
400 const mpoly_ctx_t mctx)
401 {
402 slong i, j, k;
403 slong m = I->mvars;
404 slong * perm = I->zippel2_perm;
405 slong max_main_degree;
406 double Glength, Glead_count_X, Gtail_count_X, Glead_count_Y, Gtail_count_Y;
407 double evals, bivar, reconstruct;
408
409 /* need at least 3 variables */
410 if (m < 3)
411 return;
412
413 /* figure out the two main variables y_0, y_1 */
414 for (k = 0; k < 2; k++)
415 {
416 slong main_var;
417 ulong count, deg, new_count, new_deg;
418
419 main_var = k;
420 j = perm[main_var];
421 count = FLINT_MIN(I->Alead_count[j], I->Blead_count[j]);
422 deg = FLINT_MAX(I->Adeflate_deg[j], I->Bdeflate_deg[j]);
423 for (i = k + 1; i < m; i++)
424 {
425 j = perm[i];
426 new_count = FLINT_MIN(I->Alead_count[j], I->Blead_count[j]);
427 new_deg = FLINT_MAX(I->Adeflate_deg[j], I->Bdeflate_deg[j]);
428 if (new_deg + new_count/256 < deg + count/256)
429 {
430 count = new_count;
431 deg = new_deg;
432 main_var = i;
433 }
434 }
435
436 if (main_var != k)
437 {
438 slong t = perm[main_var];
439 perm[main_var] = perm[k];
440 perm[k] = t;
441 }
442 }
443
444 max_main_degree = 0;
445 for (i = 0; i < 2; i++)
446 {
447 k = perm[i];
448 max_main_degree = FLINT_MAX(max_main_degree, I->Adeflate_deg[k]);
449 max_main_degree = FLINT_MAX(max_main_degree, I->Bdeflate_deg[k]);
450 }
451
452 /* two main variables must be packed into bits = FLINT_BITS/2 */
453 if (FLINT_BIT_COUNT(max_main_degree) >= FLINT_BITS/2)
454 return;
455
456 /* estimate length of gcd */
457 Glength = 0.5*(I->Adensity + I->Bdensity);
458 for (i = 0; i < m; i++)
459 {
460 k = perm[i];
461 Glength *= UWORD(1) + (ulong)(I->Gdeflate_deg_bound[k]);
462 }
463
464 /* estimate number of lead/tail terms of G wrt X,Y */
465 {
466 double a, b;
467 double Alead_density_X, Atail_density_X;
468 double Blead_density_X, Btail_density_X;
469 double Alead_density_Y, Atail_density_Y;
470 double Blead_density_Y, Btail_density_Y;
471
472 k = perm[0];
473 a = I->Adensity*(UWORD(1) + (ulong)(I->Adeflate_deg[k]))/Alength;
474 b = I->Bdensity*(UWORD(1) + (ulong)(I->Bdeflate_deg[k]))/Blength;
475 Alead_density_X = a*I->Alead_count[k];
476 Atail_density_X = a*I->Atail_count[k];
477 Blead_density_X = b*I->Blead_count[k];
478 Btail_density_X = b*I->Btail_count[k];
479
480 k = perm[1];
481 a = I->Adensity*(UWORD(1) + (ulong)(I->Adeflate_deg[k]))/Alength;
482 b = I->Bdensity*(UWORD(1) + (ulong)(I->Bdeflate_deg[k]))/Blength;
483 Alead_density_Y = a*I->Alead_count[k];
484 Atail_density_Y = a*I->Atail_count[k];
485 Blead_density_Y = b*I->Blead_count[k];
486 Btail_density_Y = b*I->Btail_count[k];
487
488 Glead_count_X = 0.5*(Alead_density_X + Blead_density_X);
489 Gtail_count_X = 0.5*(Atail_density_X + Btail_density_X);
490 Glead_count_Y = 0.5*(Alead_density_Y + Blead_density_Y);
491 Gtail_count_Y = 0.5*(Atail_density_Y + Btail_density_Y);
492 for (i = 0; i < m; i++)
493 {
494 k = perm[i];
495 if (i != 0)
496 {
497 Glead_count_X *= UWORD(1) + (ulong)(I->Gdeflate_deg_bound[k]);
498 Gtail_count_X *= UWORD(1) + (ulong)(I->Gdeflate_deg_bound[k]);
499 }
500
501 if (i != 1)
502 {
503 Glead_count_Y *= UWORD(1) + (ulong)(I->Gdeflate_deg_bound[k]);
504 Gtail_count_Y *= UWORD(1) + (ulong)(I->Gdeflate_deg_bound[k]);
505 }
506 }
507 }
508
509 /* evaluations needed is the max length of the coefficients of G wrt X,Y */
510 {
511 double Gmax_terms_X, Gmax_terms_Y;
512
513 k = perm[0];
514 Gmax_terms_X = Glength/(UWORD(1) + (ulong)(I->Gterm_count_est[k]));
515 Gmax_terms_X = FLINT_MAX(Gmax_terms_X, Glead_count_X);
516 Gmax_terms_X = FLINT_MAX(Gmax_terms_X, Gtail_count_X);
517 Gmax_terms_X = FLINT_MAX(Gmax_terms_X, 1);
518
519 k = perm[1];
520 Gmax_terms_Y = Glength/(UWORD(1) + (ulong)(I->Gterm_count_est[k]));
521 Gmax_terms_Y = FLINT_MAX(Gmax_terms_Y, Glead_count_Y);
522 Gmax_terms_Y = FLINT_MAX(Gmax_terms_Y, Gtail_count_Y);
523 Gmax_terms_Y = FLINT_MAX(Gmax_terms_Y, 1);
524
525 evals = Gmax_terms_X*Gmax_terms_Y/(1 + Glength);
526 }
527
528 /* time for bivar gcd */
529 {
530 double te, tg, ta, tb;
531 te = tg = ta = tb = 1;
532 for (i = 0; i < 2; i++)
533 {
534 k = perm[i];
535 /* already checked reasonable degrees with max_main_degree above */
536 te *= 1 + FLINT_MAX(I->Adeflate_deg[k], I->Bdeflate_deg[k]);
537 tg *= 1 + I->Gdeflate_deg_bound[k];
538 ta *= 1 + FLINT_MAX(0, I->Adeflate_deg[k] - I->Gdeflate_deg_bound[k]);
539 tb *= 1 + FLINT_MAX(0, I->Bdeflate_deg[k] - I->Gdeflate_deg_bound[k]);
540 }
541 bivar = te + (tg + ta + tb)*0.1;
542 }
543
544 /* time for reconstruction */
545 {
546 reconstruct = I->Gterm_count_est[perm[0]];
547 reconstruct += I->Gterm_count_est[perm[1]];
548 reconstruct = Glength*Glength/(1 + reconstruct);
549 }
550
551 I->can_use |= MPOLY_GCD_USE_ZIPPEL2;
552 I->zippel2_time = 0.00000002*bivar*evals*(Alength + Blength) +
553 0.0003*reconstruct;
554 }
555
mpoly_gcd_info_measure_zippel(mpoly_gcd_info_t I,slong Alength,slong Blength,const mpoly_ctx_t mctx)556 void mpoly_gcd_info_measure_zippel(
557 mpoly_gcd_info_t I,
558 slong Alength,
559 slong Blength,
560 const mpoly_ctx_t mctx)
561 {
562 slong i, j, k;
563 slong m = I->mvars;
564 slong * perm = I->zippel_perm;
565
566 /* need at least two variables */
567 if (m < 2)
568 return;
569
570 /* figure out a main variable y_0 */
571 {
572 slong main_var;
573 ulong count, deg, new_count, new_deg;
574
575 main_var = 0;
576 j = I->zippel_perm[main_var];
577 count = FLINT_MIN(I->Atail_count[j], I->Alead_count[j]);
578 count = FLINT_MIN(count, I->Btail_count[j]);
579 count = FLINT_MIN(count, I->Blead_count[j]);
580 deg = FLINT_MAX(I->Adeflate_deg[j], I->Bdeflate_deg[j]);
581 for (i = 1; i < m; i++)
582 {
583 j = perm[i];
584 new_count = FLINT_MIN(I->Atail_count[j], I->Alead_count[j]);
585 new_count = FLINT_MIN(new_count, I->Btail_count[j]);
586 new_count = FLINT_MIN(new_count, I->Blead_count[j]);
587 new_deg = FLINT_MAX(I->Adeflate_deg[j], I->Bdeflate_deg[j]);
588 if (new_count < count || (new_count == count && new_deg < deg))
589 {
590 count = new_count;
591 deg = new_deg;
592 main_var = i;
593 }
594 }
595
596 if (main_var != 0)
597 {
598 slong t = perm[main_var];
599 perm[main_var] = perm[0];
600 perm[0] = t;
601 }
602 }
603
604 /* sort with hope that ddeg(G,y_1) >= ddeg(G,y_2) ... >= ddeg(G,y_m) */
605 for (k = 1; k + 1 < m; k++)
606 {
607 slong var;
608 ulong deg, new_deg;
609
610 var = k;
611 j = perm[var];
612 deg = FLINT_MIN(I->Adeflate_deg[j], I->Bdeflate_deg[j]);
613 for (i = k + 1; i < m; i++)
614 {
615 j = perm[i];
616 new_deg = FLINT_MIN(I->Adeflate_deg[j], I->Bdeflate_deg[j]);
617 if (new_deg > deg)
618 {
619 deg = new_deg;
620 var = i;
621 }
622 }
623 if (var != k)
624 {
625 slong t = I->zippel_perm[var];
626 perm[var] = perm[k];
627 perm[k] = t;
628 }
629 }
630
631 I->can_use |= MPOLY_GCD_USE_ZIPPEL;
632 I->zippel_time = 0.3456;
633 }
634
mpoly_gcd_info_measure_zippel2(mpoly_gcd_info_t I,slong Alength,slong Blength,const mpoly_ctx_t mctx)635 void mpoly_gcd_info_measure_zippel2(
636 mpoly_gcd_info_t I,
637 slong Alength,
638 slong Blength,
639 const mpoly_ctx_t mctx)
640 {
641 slong i, j, k;
642 slong m = I->mvars;
643 slong * perm = I->zippel2_perm;
644 slong max_main_degree;
645
646 /* need at least 3 variables */
647 if (m < 3)
648 return;
649
650 /* figure out the two main variables y_0, y_1 */
651
652 #define NEEDS_SWAP \
653 FLINT_MIN(I->Adeflate_deg[perm[j]], I->Bdeflate_deg[perm[j]]) < \
654 FLINT_MIN(I->Adeflate_deg[perm[j-1]], I->Bdeflate_deg[perm[j-1]]) \
655
656 for (i = 1; i < m; i++)
657 for (j = i; j > 0 && NEEDS_SWAP; j--)
658 SLONG_SWAP(perm[j], perm[j - 1]);
659
660
661 #define NEEDS_SWAP2 \
662 FLINT_MIN(I->Adeflate_deg[perm[j]], I->Bdeflate_deg[perm[j]]) > \
663 FLINT_MIN(I->Adeflate_deg[perm[j-1]], I->Bdeflate_deg[perm[j-1]]) \
664
665 for (i = 3; i < m; i++)
666 for (j = i; j > 2 && NEEDS_SWAP; j--)
667 SLONG_SWAP(perm[j], perm[j - 1]);
668
669
670 max_main_degree = 0;
671 for (i = 0; i < 2; i++)
672 {
673 k = perm[i];
674 max_main_degree = FLINT_MAX(max_main_degree, I->Adeflate_deg[k]);
675 max_main_degree = FLINT_MAX(max_main_degree, I->Bdeflate_deg[k]);
676 }
677
678 /* two main variables must be packed into bits = FLINT_BITS/2 */
679 if (FLINT_BIT_COUNT(max_main_degree) >= FLINT_BITS/2)
680 return;
681
682 I->can_use |= MPOLY_GCD_USE_ZIPPEL2;
683 I->zippel2_time = 0.243;
684 }
685