1 /*
2 Copyright (C) 2016-2017 William Hart
3 Copyright (C) 2017-2020 Daniel Schultz
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <https://www.gnu.org/licenses/>.
11 */
12
13 #ifndef MPOLY_H
14 #define MPOLY_H
15
16 #ifdef MPOLY_INLINES_C
17 #define MPOLY_INLINE FLINT_DLL
18 #else
19 #define MPOLY_INLINE static __inline__
20 #endif
21
22 #undef ulong
23 #define ulong ulongxx /* interferes with system includes */
24 #include <stdio.h>
25 #undef ulong
26
27 #include <gmp.h>
28 #define ulong mp_limb_t
29
30 #include "string.h"
31 #include "flint.h"
32 #include "fmpz.h"
33 #include "fmpz_mat.h"
34 #include "fmpq.h"
35 #include "ulong_extras.h"
36 #include "thread_pool.h"
37
38 #ifdef __cplusplus
39 extern "C" {
40 #endif
41
42 #define MPOLY_MIN_BITS (UWORD(8)) /* minimum number of bits to pack into */
43
44 /* choose m so that (m + 1)/(n - m) ~= la/lb, i.e. m = (n*la - lb)/(la + lb) */
mpoly_divide_threads(slong n,double la,double lb)45 MPOLY_INLINE slong mpoly_divide_threads(slong n, double la, double lb)
46 {
47 double m_double = (n*la - lb)/(la + lb);
48 slong m = m_double + (2*m_double > n ? -0.5 : 0.5);
49
50 /* input must satisfy */
51 FLINT_ASSERT(n > 0);
52
53 if (m <= 0)
54 m = 0;
55
56 if (m >= n - 1)
57 m = n - 1;
58
59 /* output must satisfy */
60 FLINT_ASSERT(m >= 0);
61 FLINT_ASSERT(m < n);
62 return m;
63 }
64
65
66 typedef enum {
67 ORD_LEX,
68 ORD_DEGLEX,
69 ORD_DEGREVLEX
70 } ordering_t;
71
72 #define MPOLY_NUM_ORDERINGS 3
73
74 /* context *******************************************************************/
75
76 typedef struct
77 {
78 slong nvars; /* number of variables */
79 slong nfields; /* number of fields in exponent vector */
80 ordering_t ord; /* monomial ordering */
81 int deg; /* is ord a degree ordering? */
82 int rev; /* is ord a reversed ordering? */
83 slong lut_words_per_exp[FLINT_BITS];
84 unsigned char lut_fix_bits[FLINT_BITS]; /* FLINT_BITS < 256 */
85 } mpoly_ctx_struct;
86
87 typedef mpoly_ctx_struct mpoly_ctx_t[1];
88
89 FLINT_DLL void mpoly_ctx_init(mpoly_ctx_t ctx, slong nvars, const ordering_t ord);
90
91 FLINT_DLL void mpoly_ctx_init_rand(mpoly_ctx_t mctx, flint_rand_t state, slong max_nvars);
92
93 FLINT_DLL void mpoly_monomial_randbits_fmpz(fmpz * exp, flint_rand_t state, flint_bitcnt_t exp_bits, const mpoly_ctx_t mctx);
94
95 FLINT_DLL void mpoly_ctx_clear(mpoly_ctx_t mctx);
96
97 /*
98 number of words used by an exponent vector packed into "bits" bits:
99 we must have either
100 (mp) bits > FLINT_BITS and bits % FLINT_BITS == 0, or
101 (sp) MPOLY_MIN_BITS <= bits <= FLINT_BITS
102 */
103 MPOLY_INLINE
mpoly_words_per_exp_sp(flint_bitcnt_t bits,const mpoly_ctx_t mctx)104 slong mpoly_words_per_exp_sp(flint_bitcnt_t bits, const mpoly_ctx_t mctx)
105 {
106 FLINT_ASSERT(0 < bits);
107 FLINT_ASSERT(bits <= FLINT_BITS);
108 FLINT_ASSERT(mctx->lut_words_per_exp[bits - 1]
109 == (mctx->nfields - 1)/(FLINT_BITS/bits) + 1);
110 return mctx->lut_words_per_exp[bits - 1];
111 }
112
113 MPOLY_INLINE
mpoly_words_per_exp_mp(flint_bitcnt_t bits,const mpoly_ctx_t mctx)114 slong mpoly_words_per_exp_mp(flint_bitcnt_t bits, const mpoly_ctx_t mctx)
115 {
116 FLINT_ASSERT(bits % FLINT_BITS == 0);
117 return bits/FLINT_BITS*mctx->nfields;
118 }
119
120 MPOLY_INLINE
mpoly_words_per_exp(flint_bitcnt_t bits,const mpoly_ctx_t mctx)121 slong mpoly_words_per_exp(flint_bitcnt_t bits, const mpoly_ctx_t mctx)
122 {
123 if (bits <= FLINT_BITS)
124 return mpoly_words_per_exp_sp(bits, mctx);
125 else
126 return mpoly_words_per_exp_mp(bits, mctx);
127 }
128
129 /*
130 If "bits" is simply the number of bits needed to pack an exponent vector,
131 possibly upgrade it so that it is either
132 (mp) a multiple of FLINT_BITS in the mp case, or
133 (sp) as big as possible without increasing words_per_exp in the sp case
134 The upgrade in (mp) is manditory, while the upgrade in (sp) is simply nice.
135 */
136 MPOLY_INLINE
mpoly_fix_bits(flint_bitcnt_t bits,const mpoly_ctx_t mctx)137 flint_bitcnt_t mpoly_fix_bits(flint_bitcnt_t bits, const mpoly_ctx_t mctx)
138 {
139 FLINT_ASSERT(bits > 0);
140 if (bits <= FLINT_BITS)
141 return mctx->lut_fix_bits[bits - 1];
142 else
143 return (bits + FLINT_BITS - 1)/FLINT_BITS*FLINT_BITS;
144 }
145
146
147 /* heaps *********************************************************************/
148 typedef struct mpoly_heap_t
149 {
150 ulong i;
151 ulong j;
152 struct mpoly_heap_t * next;
153 } mpoly_heap_t;
154
155 typedef struct mpoly_nheap_t
156 {
157 ulong i;
158 ulong j;
159 struct mpoly_nheap_t * next;
160 slong p;
161 } mpoly_nheap_t;
162
163 typedef struct mpoly_heap1_s
164 {
165 ulong exp;
166 void * next;
167 } mpoly_heap1_s;
168
169 typedef struct mpoly_heap_s
170 {
171 ulong * exp;
172 void * next;
173 } mpoly_heap_s;
174
175 /* trees *********************************************************************/
176
177 /* red-black with ui keys */
178 typedef struct {
179 ulong key;
180 slong up;
181 slong left;
182 slong right;
183 int color;
184 } mpoly_rbnode_ui_struct;
185
186 typedef struct {
187 slong length;
188 mpoly_rbnode_ui_struct * nodes;
189 slong node_alloc;
190 char * data;
191 slong data_alloc;
192 slong data_size;
193 } mpoly_rbtree_ui_struct;
194
195 typedef mpoly_rbtree_ui_struct mpoly_rbtree_ui_t[1];
196
197 FLINT_DLL void mpoly_rbtree_ui_init(mpoly_rbtree_ui_t T, slong data_size);
198
199 FLINT_DLL void mpoly_rbtree_ui_clear(mpoly_rbtree_ui_t T);
200
201 FLINT_DLL void * mpoly_rbtree_ui_lookup(mpoly_rbtree_ui_t T, int * its_new,
202 ulong key);
203
mpoly_rbtree_ui_head(const mpoly_rbtree_ui_t T)204 MPOLY_INLINE slong mpoly_rbtree_ui_head(const mpoly_rbtree_ui_t T)
205 {
206 FLINT_ASSERT(T->nodes[1].left >= 0 || T->length < 1);
207 return T->nodes[1].left;
208 }
209
210 /* red-black with fmpz keys */
211 typedef struct {
212 fmpz_t key;
213 slong up;
214 slong left;
215 slong right;
216 int color;
217 } mpoly_rbnode_fmpz_struct;
218
219 typedef struct {
220 slong length;
221 mpoly_rbnode_fmpz_struct * nodes;
222 slong node_alloc;
223 char * data;
224 slong data_alloc;
225 slong data_size;
226 } mpoly_rbtree_fmpz_struct;
227
228 typedef mpoly_rbtree_fmpz_struct mpoly_rbtree_fmpz_t[1];
229
230 FLINT_DLL void mpoly_rbtree_fmpz_init(mpoly_rbtree_fmpz_t T, slong data_size);
231
232 FLINT_DLL void mpoly_rbtree_fmpz_clear(mpoly_rbtree_fmpz_t T);
233
234 FLINT_DLL void * mpoly_rbtree_fmpz_lookup(mpoly_rbtree_fmpz_t T, int * its_new,
235 const fmpz_t key);
236
mpoly_rbtree_fmpz_head(const mpoly_rbtree_fmpz_t T)237 MPOLY_INLINE slong mpoly_rbtree_fmpz_head(const mpoly_rbtree_fmpz_t T)
238 {
239 FLINT_ASSERT(T->nodes[1].left >= 0 || T->length < 1);
240 return T->nodes[1].left;
241 }
242
243 /* Orderings *****************************************************************/
244
245 MPOLY_INLINE
mpoly_ordering_randtest(flint_rand_t state)246 ordering_t mpoly_ordering_randtest(flint_rand_t state)
247 {
248 return (ordering_t) n_randint(state, MPOLY_NUM_ORDERINGS);
249 }
250
251 MPOLY_INLINE
mpoly_ordering_isdeg(const mpoly_ctx_t mctx)252 int mpoly_ordering_isdeg(const mpoly_ctx_t mctx)
253 {
254 return mctx->ord == ORD_DEGLEX || mctx->ord == ORD_DEGREVLEX;
255 }
256
257 MPOLY_INLINE
mpoly_ordering_isrev(const mpoly_ctx_t mctx)258 int mpoly_ordering_isrev(const mpoly_ctx_t mctx)
259 {
260 return mctx->ord == ORD_DEGREVLEX;
261 }
262
263 MPOLY_INLINE
mpoly_ordering_print(ordering_t ord)264 void mpoly_ordering_print(ordering_t ord)
265 {
266 switch (ord)
267 {
268 case ORD_LEX:
269 printf("lex");
270 break;
271 case ORD_DEGLEX:
272 printf("deglex");
273 break;
274 case ORD_DEGREVLEX:
275 printf("degrevlex");
276 break;
277 default:
278 printf("Unknown ordering in mpoly_ordering_print.");
279 }
280 }
281
282 /* Monomials ****************************************************************/
283
284 MPOLY_INLINE
mpoly_monomial_zero(ulong * exp_ptr,slong N)285 void mpoly_monomial_zero(ulong * exp_ptr, slong N)
286 {
287 slong i;
288 for (i = 0; i < N; i++)
289 exp_ptr[i] = 0;
290 }
291
292 MPOLY_INLINE
mpoly_monomial_add(ulong * exp_ptr,const ulong * exp2,const ulong * exp3,slong N)293 void mpoly_monomial_add(ulong * exp_ptr, const ulong * exp2,
294 const ulong * exp3, slong N)
295 {
296 slong i;
297 for (i = 0; i < N; i++)
298 exp_ptr[i] = exp2[i] + exp3[i];
299 }
300
301 MPOLY_INLINE
mpoly_monomial_add_mp(ulong * exp_ptr,const ulong * exp2,const ulong * exp3,slong N)302 void mpoly_monomial_add_mp(ulong * exp_ptr, const ulong * exp2,
303 const ulong * exp3, slong N)
304 {
305 mpn_add_n(exp_ptr, exp2, exp3, N);
306 }
307
308 MPOLY_INLINE
mpoly_monomial_sub(ulong * exp_ptr,const ulong * exp2,const ulong * exp3,slong N)309 void mpoly_monomial_sub(ulong * exp_ptr, const ulong * exp2,
310 const ulong * exp3, slong N)
311 {
312 slong i;
313 for (i = 0; i < N; i++)
314 exp_ptr[i] = exp2[i] - exp3[i];
315 }
316
317 MPOLY_INLINE
mpoly_monomial_sub_mp(ulong * exp_ptr,const ulong * exp2,const ulong * exp3,slong N)318 void mpoly_monomial_sub_mp(ulong * exp_ptr, const ulong * exp2,
319 const ulong * exp3, slong N)
320 {
321 mpn_sub_n(exp_ptr, exp2, exp3, N);
322 }
323
324 MPOLY_INLINE
mpoly_monomial_madd(ulong * exp1,const ulong * exp2,ulong scalar,const ulong * exp3,slong N)325 void mpoly_monomial_madd(ulong * exp1, const ulong * exp2, ulong scalar,
326 const ulong * exp3, slong N)
327 {
328 slong i;
329 for (i = 0; i < N; i++)
330 exp1[i] = exp2[i] + scalar*exp3[i];
331 }
332
333 MPOLY_INLINE
mpoly_monomial_madd_mp(ulong * exp1,const ulong * exp2,ulong scalar,const ulong * exp3,slong N)334 void mpoly_monomial_madd_mp(ulong * exp1, const ulong * exp2, ulong scalar,
335 const ulong * exp3, slong N)
336 {
337 slong i;
338 for (i = 0; i < N; i++)
339 exp1[i] = exp2[i];
340 mpn_addmul_1(exp1, exp3, N, scalar);
341 }
342
343 MPOLY_INLINE
mpoly_monomial_madd_inplace_mp(ulong * exp12,ulong scalar,const ulong * exp3,slong N)344 void mpoly_monomial_madd_inplace_mp(ulong * exp12, ulong scalar,
345 const ulong * exp3, slong N)
346 {
347 mpn_addmul_1(exp12, exp3, N, scalar);
348 }
349
350 MPOLY_INLINE
mpoly_monomial_msub(ulong * exp1,const ulong * exp2,ulong scalar,const ulong * exp3,slong N)351 void mpoly_monomial_msub(ulong * exp1, const ulong * exp2, ulong scalar,
352 const ulong * exp3, slong N)
353 {
354 slong i;
355 for (i = 0; i < N; i++)
356 exp1[i] = exp2[i] - scalar*exp3[i];
357 }
358
359 MPOLY_INLINE
mpoly_monomial_msub_mp(ulong * exp1,const ulong * exp2,ulong scalar,const ulong * exp3,slong N)360 void mpoly_monomial_msub_mp(ulong * exp1, const ulong * exp2, ulong scalar,
361 const ulong * exp3, slong N)
362 {
363 slong i;
364 for (i = 0; i < N; i++)
365 exp1[i] = exp2[i];
366 FLINT_ASSERT(N > 0);
367 mpn_submul_1(exp1, exp3, N, scalar);
368 }
369
370 MPOLY_INLINE
mpoly_monomial_msub_ui_array(ulong * exp1,const ulong * exp2,const ulong * scalar,slong scalar_limbs,const ulong * exp3,slong N)371 void mpoly_monomial_msub_ui_array(ulong * exp1, const ulong * exp2,
372 const ulong * scalar, slong scalar_limbs,
373 const ulong * exp3, slong N)
374 {
375 slong i;
376 for (i = 0; i < N; i++)
377 exp1[i] = exp2[i];
378 FLINT_ASSERT(scalar_limbs <= N);
379 for (i = 0; i < scalar_limbs; i++)
380 {
381 FLINT_ASSERT(N > i);
382 mpn_submul_1(exp1 + i, exp3, N - i, scalar[i]);
383 }
384 }
385
386 MPOLY_INLINE
mpoly_monomial_madd_ui_array(ulong * exp1,const ulong * exp2,const ulong * scalar,slong scalar_limbs,const ulong * exp3,slong N)387 void mpoly_monomial_madd_ui_array(ulong * exp1, const ulong * exp2,
388 const ulong * scalar, slong scalar_limbs,
389 const ulong * exp3, slong N)
390 {
391 slong i;
392 for (i = 0; i < N; i++)
393 exp1[i] = exp2[i];
394 FLINT_ASSERT(scalar_limbs <= N);
395 for (i = 0; i < scalar_limbs; i++)
396 mpn_addmul_1(exp1 + i, exp3, N - i, scalar[i]);
397 }
398
399 MPOLY_INLINE
mpoly_monomial_madd_fmpz(ulong * exp1,const ulong * exp2,const fmpz_t scalar,const ulong * exp3,slong N)400 void mpoly_monomial_madd_fmpz(ulong * exp1, const ulong * exp2,
401 const fmpz_t scalar, const ulong * exp3, slong N)
402 {
403 if (COEFF_IS_MPZ(*scalar))
404 {
405 __mpz_struct * mpz = COEFF_TO_PTR(*scalar);
406 mpoly_monomial_madd_ui_array(exp1, exp2,
407 mpz->_mp_d, mpz->_mp_size, exp3, N);
408 }
409 else
410 {
411 mpoly_monomial_madd_mp(exp1, exp2, *scalar, exp3, N);
412 }
413 }
414
415 MPOLY_INLINE
mpoly_overflow_mask_sp(flint_bitcnt_t bits)416 ulong mpoly_overflow_mask_sp(flint_bitcnt_t bits)
417 {
418 slong i;
419 ulong mask = 0;
420
421 FLINT_ASSERT(bits <= FLINT_BITS);
422
423 for (i = 0; i < FLINT_BITS/bits; i++)
424 mask = (mask << bits) + (UWORD(1) << (bits - 1));
425
426 return mask;
427 }
428
429 MPOLY_INLINE
mpoly_monomial_max1(ulong exp2,ulong exp3,flint_bitcnt_t bits,ulong mask)430 ulong mpoly_monomial_max1(ulong exp2, ulong exp3,
431 flint_bitcnt_t bits, ulong mask)
432 {
433 ulong s, m, exp1;
434 s = mask + exp2 - exp3;
435 m = mask & s;
436 m = m - (m >> (bits - 1));
437 exp1 = exp3 + (s & m);
438 return exp1;
439 }
440
441 MPOLY_INLINE
mpoly_monomial_max(ulong * exp1,const ulong * exp2,const ulong * exp3,flint_bitcnt_t bits,slong N,ulong mask)442 void mpoly_monomial_max(ulong * exp1, const ulong * exp2, const ulong * exp3,
443 flint_bitcnt_t bits, slong N, ulong mask)
444 {
445 ulong i, s, m;
446 for (i = 0; i < N; i++)
447 {
448 s = mask + exp2[i] - exp3[i];
449 m = mask & s;
450 m = m - (m >> (bits - 1));
451 exp1[i] = exp3[i] + (s & m);
452 }
453 }
454
455 MPOLY_INLINE
mpoly_monomial_min1(ulong exp2,ulong exp3,flint_bitcnt_t bits,ulong mask)456 ulong mpoly_monomial_min1(ulong exp2, ulong exp3,
457 flint_bitcnt_t bits, ulong mask)
458 {
459 ulong s, m, exp1;
460 s = mask + exp2 - exp3;
461 m = mask & s;
462 m = m - (m >> (bits - 1));
463 exp1 = exp2 - (s & m);
464 return exp1;
465 }
466
467 MPOLY_INLINE
mpoly_monomial_min(ulong * exp1,const ulong * exp2,const ulong * exp3,flint_bitcnt_t bits,slong N,ulong mask)468 void mpoly_monomial_min(ulong * exp1, const ulong * exp2, const ulong * exp3,
469 flint_bitcnt_t bits, slong N, ulong mask)
470 {
471 ulong i, s, m;
472 for (i = 0; i < N; i++)
473 {
474 s = mask + exp2[i] - exp3[i];
475 m = mask & s;
476 m = m - (m >> (bits - 1));
477 exp1[i] = exp2[i] - (s & m);
478 }
479 }
480
481 MPOLY_INLINE
mpoly_monomial_max_mp(ulong * exp1,const ulong * exp2,const ulong * exp3,flint_bitcnt_t bits,slong N)482 void mpoly_monomial_max_mp(ulong * exp1, const ulong * exp2, const ulong * exp3,
483 flint_bitcnt_t bits, slong N)
484 {
485 slong i, j;
486 for (i = 0; i < N; i += bits/FLINT_BITS)
487 {
488 const ulong * t = exp2;
489 for (j = bits/FLINT_BITS - 1; j >= 0; j--)
490 {
491 if (exp3[i + j] != exp2[i + j])
492 {
493 if (exp3[i + j] > exp2[i + j])
494 t = exp3;
495 break;
496 }
497 }
498 for (j = 0; j < bits/FLINT_BITS; j++)
499 {
500 exp1[i + j] = t[i + j];
501 }
502 }
503 }
504
505 MPOLY_INLINE
mpoly_monomial_min_mp(ulong * exp1,const ulong * exp2,const ulong * exp3,flint_bitcnt_t bits,slong N)506 void mpoly_monomial_min_mp(ulong * exp1, const ulong * exp2, const ulong * exp3,
507 flint_bitcnt_t bits, slong N)
508 {
509 slong i, j;
510 for (i = 0; i < N; i += bits/FLINT_BITS)
511 {
512 const ulong * t = exp2;
513 for (j = bits/FLINT_BITS - 1; j >= 0; j--)
514 {
515 if (exp3[i + j] != exp2[i + j])
516 {
517 if (exp3[i + j] < exp2[i + j])
518 t = exp3;
519 break;
520 }
521 }
522 for (j = 0; j < bits/FLINT_BITS; j++)
523 {
524 exp1[i + j] = t[i + j];
525 }
526 }
527 }
528
529 MPOLY_INLINE
mpoly_monomial_overflows(ulong * exp2,slong N,ulong mask)530 int mpoly_monomial_overflows(ulong * exp2, slong N, ulong mask)
531 {
532 slong i;
533 for (i = 0; i < N; i++)
534 {
535 if ((exp2[i] & mask) != 0)
536 return 1;
537 }
538 return 0;
539 }
540
541 MPOLY_INLINE
mpoly_monomial_overflows_mp(ulong * exp_ptr,slong N,flint_bitcnt_t bits)542 int mpoly_monomial_overflows_mp(ulong * exp_ptr, slong N, flint_bitcnt_t bits)
543 {
544 slong i = bits/FLINT_BITS - 1;
545 do {
546 if ((slong)(exp_ptr[i]) < 0)
547 return 1;
548 i += bits/FLINT_BITS;
549 } while (i < N);
550
551 return 0;
552 }
553
554 MPOLY_INLINE
mpoly_monomial_overflows1(ulong exp,ulong mask)555 int mpoly_monomial_overflows1(ulong exp, ulong mask)
556 {
557 return (exp & mask) != 0;
558 }
559
560 MPOLY_INLINE
mpoly_monomial_divides(ulong * exp_ptr,const ulong * exp2,const ulong * exp3,slong N,ulong mask)561 int mpoly_monomial_divides(ulong * exp_ptr, const ulong * exp2,
562 const ulong * exp3, slong N, ulong mask)
563 {
564 slong i;
565 for (i = 0; i < N; i++)
566 {
567 exp_ptr[i] = exp2[i] - exp3[i];
568
569 if ((exp_ptr[i] & mask) != 0)
570 return 0;
571 }
572
573 return 1;
574 }
575
576 MPOLY_INLINE
mpoly_monomial_halves(ulong * exp_ptr,const ulong * exp2,slong N,ulong mask)577 int mpoly_monomial_halves(ulong * exp_ptr, const ulong * exp2,
578 slong N, ulong mask)
579 {
580 slong i;
581 ulong bw;
582
583 bw = mpn_rshift(exp_ptr, exp2, N, 1);
584
585 if (bw != 0)
586 return 0;
587
588 for (i = 0; i < N; i++)
589 {
590 if ((exp_ptr[i] & mask) != 0)
591 return 0;
592 }
593
594 return 1;
595 }
596
597 MPOLY_INLINE
mpoly_monomial_divides_mp(ulong * exp_ptr,const ulong * exp2,const ulong * exp3,slong N,flint_bitcnt_t bits)598 int mpoly_monomial_divides_mp(ulong * exp_ptr, const ulong * exp2,
599 const ulong * exp3, slong N, flint_bitcnt_t bits)
600 {
601 slong i;
602
603 mpn_sub_n(exp_ptr, exp2, exp3, N);
604
605 i = bits/FLINT_BITS - 1;
606 do {
607 if ((slong)(exp_ptr[i]) < 0)
608 return 0;
609 i += bits/FLINT_BITS;
610 } while (i < N);
611
612 return 1;
613 }
614
615 MPOLY_INLINE
mpoly_monomial_halves_mp(ulong * exp_ptr,const ulong * exp2,slong N,flint_bitcnt_t bits)616 int mpoly_monomial_halves_mp(ulong * exp_ptr, const ulong * exp2,
617 slong N, flint_bitcnt_t bits)
618 {
619 slong i;
620 ulong bw;
621
622 bw = mpn_rshift(exp_ptr, exp2, N, 1);
623
624 if (bw != 0)
625 return 0;
626
627 i = bits/FLINT_BITS - 1;
628 do {
629 if ((slong)(exp_ptr[i]) < 0)
630 return 0;
631 i += bits/FLINT_BITS;
632 } while (i < N);
633
634 return 1;
635 }
636
637 MPOLY_INLINE
mpoly_monomial_divides_test(const ulong * exp2,const ulong * exp3,slong N,ulong mask)638 int mpoly_monomial_divides_test(const ulong * exp2,
639 const ulong * exp3, slong N, ulong mask)
640 {
641 slong i;
642 for (i = 0; i < N; i++)
643 if (((exp2[i] - exp3[i]) & mask) != 0)
644 return 0;
645
646 return 1;
647 }
648
649 MPOLY_INLINE
mpoly_monomial_divides_mp_test(const ulong * exp2,const ulong * exp3,slong N,flint_bitcnt_t bits)650 int mpoly_monomial_divides_mp_test(const ulong * exp2,
651 const ulong * exp3, slong N, flint_bitcnt_t bits)
652 {
653 slong i, j;
654
655 i = 0;
656 do {
657 for (j = bits/FLINT_BITS - 1; j >= 0; j--)
658 {
659 if (exp2[i + j] > exp3[i + j])
660 break;
661 if (exp2[i + j] < exp3[i + j])
662 return 0;
663 }
664 i += bits/FLINT_BITS;
665 } while (i < N);
666
667 return 1;
668 }
669
670 MPOLY_INLINE
mpoly_monomial_divides1(ulong * exp_ptr,const ulong exp2,const ulong exp3,ulong mask)671 int mpoly_monomial_divides1(ulong * exp_ptr, const ulong exp2,
672 const ulong exp3, ulong mask)
673 {
674 (*exp_ptr) = exp2 - exp3;
675
676 if (((exp2 - exp3) & mask) != 0)
677 return 0;
678
679 return 1;
680 }
681
682 MPOLY_INLINE
mpoly_monomial_halves1(ulong * exp_ptr,const ulong exp2,ulong mask)683 int mpoly_monomial_halves1(ulong * exp_ptr, const ulong exp2, ulong mask)
684 {
685 if (exp2 & 1)
686 return 0;
687
688 (*exp_ptr) = exp2 >> 1;
689
690 if (((exp2 >> 1) & mask) != 0)
691 return 0;
692
693 return 1;
694 }
695
696 MPOLY_INLINE
mpoly_monomial_set(ulong * exp2,const ulong * exp3,slong N)697 void mpoly_monomial_set(ulong * exp2, const ulong * exp3, slong N)
698 {
699 slong i;
700 for (i = 0; i < N; i++)
701 exp2[i] = exp3[i];
702 }
703
704 MPOLY_INLINE
mpoly_monomial_set_extra(ulong * exp2,const ulong * exp3,slong N,slong offset,ulong extra)705 void mpoly_monomial_set_extra(ulong * exp2, const ulong * exp3,
706 slong N, slong offset, ulong extra)
707 {
708 slong i;
709 for (i = 0; i < N; i++)
710 {
711 exp2[i] = exp3[i] + (i == offset ? extra : 0);
712 }
713 }
714
715 MPOLY_INLINE
mpoly_copy_monomials(ulong * exp1,const ulong * exp2,slong len,slong N)716 void mpoly_copy_monomials(ulong * exp1, const ulong * exp2, slong len, slong N)
717 {
718 memcpy(exp1, exp2, N*len*sizeof(ulong));
719 }
720
721 MPOLY_INLINE
mpoly_monomial_swap(ulong * exp2,ulong * exp3,slong N)722 void mpoly_monomial_swap(ulong * exp2, ulong * exp3, slong N)
723 {
724 slong i;
725 ulong t;
726 for (i = 0; i < N; i++)
727 {
728 t = exp2[i];
729 exp2[i] = exp3[i];
730 exp3[i] = t;
731 }
732 }
733
734 MPOLY_INLINE
mpoly_monomial_mul_ui(ulong * exp2,const ulong * exp3,slong N,ulong c)735 void mpoly_monomial_mul_ui(ulong * exp2, const ulong * exp3, slong N, ulong c)
736 {
737 slong i;
738 for (i = 0; i < N; i++)
739 exp2[i] = exp3[i]*c;
740 }
741
742 MPOLY_INLINE
mpoly_monomial_mul_ui_mp(ulong * exp2,const ulong * exp3,slong N,ulong c)743 void mpoly_monomial_mul_ui_mp(ulong * exp2, const ulong * exp3, slong N, ulong c)
744 {
745 FLINT_ASSERT(N > 0);
746 mpn_mul_1(exp2, exp3, N, c);
747 }
748
749 FLINT_DLL void mpoly_monomial_mul_fmpz(ulong * exp2, const ulong * exp3,
750 slong N, const fmpz_t c);
751
752 MPOLY_INLINE
mpoly_monomial_is_zero(const ulong * exp,slong N)753 int mpoly_monomial_is_zero(const ulong * exp, slong N)
754 {
755 slong i;
756 for (i = 0; i < N; i++)
757 {
758 if (exp[i] != 0)
759 return 0;
760 }
761
762 return 1;
763 }
764
765 MPOLY_INLINE
mpoly_monomial_equal(const ulong * exp2,const ulong * exp3,slong N)766 int mpoly_monomial_equal(const ulong * exp2, const ulong * exp3, slong N)
767 {
768 slong i;
769
770 for (i = 0; i < N; i++)
771 {
772 if (exp2[i] != exp3[i])
773 return 0;
774 }
775
776 return 1;
777 }
778
779 MPOLY_INLINE
mpoly_monomial_equal_extra(const ulong * exp2,const ulong * exp3,slong N,slong offset,ulong extra)780 int mpoly_monomial_equal_extra(const ulong * exp2, const ulong * exp3,
781 slong N, slong offset, ulong extra)
782 {
783 slong i;
784
785 for (i = 0; i < N; i++)
786 {
787 ulong e3 = exp3[i] + ((i == offset) ? extra : 0);
788 if (exp2[i] != e3)
789 return 0;
790 }
791
792 return 1;
793 }
794
795 MPOLY_INLINE
mpoly_monomial_cmp1(ulong a,ulong b,ulong cmpmask)796 int mpoly_monomial_cmp1(ulong a, ulong b, ulong cmpmask)
797 {
798 if ((a^cmpmask) != (b^cmpmask))
799 {
800 if ((a^cmpmask) > (b^cmpmask))
801 return 1;
802 else
803 return -1;
804 }
805 return 0;
806 }
807
808 MPOLY_INLINE
mpoly_monomial_gt1(ulong a,ulong b,ulong cmpmask)809 int mpoly_monomial_gt1(ulong a, ulong b, ulong cmpmask)
810 {
811 return (a^cmpmask) > (b^cmpmask);
812 }
813
814 MPOLY_INLINE
mpoly_monomial_ge1(ulong a,ulong b,ulong cmpmask)815 int mpoly_monomial_ge1(ulong a, ulong b, ulong cmpmask)
816 {
817 return (a^cmpmask) >= (b^cmpmask);
818 }
819
820 MPOLY_INLINE
mpoly_monomial_lt(const ulong * exp3,const ulong * exp2,slong N,const ulong * cmpmask)821 int mpoly_monomial_lt(const ulong * exp3, const ulong * exp2,
822 slong N, const ulong * cmpmask)
823 {
824 slong i = N - 1;
825 do {
826 if (exp2[i] != exp3[i])
827 {
828 return (exp3[i]^cmpmask[i]) < (exp2[i]^cmpmask[i]);
829 }
830 } while (--i >= 0);
831 return 0;
832 }
833
834 MPOLY_INLINE
mpoly_monomial_gt(const ulong * exp3,const ulong * exp2,slong N,const ulong * cmpmask)835 int mpoly_monomial_gt(const ulong * exp3, const ulong * exp2,
836 slong N, const ulong * cmpmask)
837 {
838 slong i = N - 1;
839 do {
840 if (exp2[i] != exp3[i])
841 {
842 return (exp3[i]^cmpmask[i]) > (exp2[i]^cmpmask[i]);
843 }
844 } while (--i >= 0);
845 return 0;
846 }
847
848 MPOLY_INLINE
mpoly_monomial_lt_nomask(const ulong * exp2,const ulong * exp3,slong N)849 int mpoly_monomial_lt_nomask(const ulong * exp2, const ulong * exp3, slong N)
850 {
851 slong i = N - 1;
852 do {
853 if (exp2[i] != exp3[i])
854 {
855 return exp2[i] < exp3[i];
856 }
857 } while (--i >= 0);
858 return 0;
859 }
860
861 MPOLY_INLINE
mpoly_monomial_gt_nomask(const ulong * exp2,const ulong * exp3,slong N)862 int mpoly_monomial_gt_nomask(const ulong * exp2, const ulong * exp3, slong N)
863 {
864 slong i = N - 1;
865 do {
866 if (exp2[i] != exp3[i])
867 {
868 return exp2[i] > exp3[i];
869 }
870 } while (--i >= 0);
871 return 0;
872 }
873
874 MPOLY_INLINE
mpoly_monomial_lt_nomask_extra(const ulong * exp2,const ulong * exp3,slong N,slong offset,ulong extra)875 int mpoly_monomial_lt_nomask_extra(const ulong * exp2, const ulong * exp3,
876 slong N, slong offset, ulong extra)
877 {
878 slong i = N - 1;
879 do {
880 ulong e3 = exp3[i] + ((i == offset) ? extra : 0);
881 if (exp2[i] != e3)
882 {
883 return exp2[i] < e3;
884 }
885 } while (--i >= 0);
886 return 0;
887 }
888
889 MPOLY_INLINE
mpoly_monomial_gt_nomask_extra(const ulong * exp2,const ulong * exp3,slong N,slong offset,ulong extra)890 int mpoly_monomial_gt_nomask_extra(const ulong * exp2, const ulong * exp3,
891 slong N, slong offset, ulong extra)
892 {
893 slong i = N - 1;
894 do {
895 ulong e3 = exp3[i] + ((i == offset) ? extra : 0);
896 if (exp2[i] != e3)
897 {
898 return exp2[i] > e3;
899 }
900 } while (--i >= 0);
901 return 0;
902 }
903
904 MPOLY_INLINE
mpoly_monomial_cmp(const ulong * exp2,const ulong * exp3,slong N,const ulong * cmpmask)905 int mpoly_monomial_cmp(const ulong * exp2, const ulong * exp3,
906 slong N, const ulong * cmpmask)
907 {
908 slong i = N - 1;
909 do {
910 if (exp2[i] != exp3[i])
911 {
912 if ((exp2[i]^cmpmask[i]) > (exp3[i]^cmpmask[i]))
913 return 1;
914 else
915 return -1;
916 }
917 } while (--i >= 0);
918 return 0;
919 }
920
921 MPOLY_INLINE
mpoly_monomial_cmp_nomask(const ulong * exp2,const ulong * exp3,slong N)922 int mpoly_monomial_cmp_nomask(const ulong * exp2, const ulong * exp3, slong N)
923 {
924 slong i = N - 1;
925 do {
926 if (exp2[i] != exp3[i])
927 {
928 if (exp2[i] > exp3[i])
929 return 1;
930 else
931 return -1;
932 }
933 } while (--i >= 0);
934 return 0;
935 }
936
937 MPOLY_INLINE
mpoly_monomial_cmp_nomask_extra(const ulong * exp2,const ulong * exp3,slong N,slong offset,ulong extra)938 int mpoly_monomial_cmp_nomask_extra(const ulong * exp2, const ulong * exp3,
939 slong N, slong offset, ulong extra)
940 {
941 slong i = N - 1;
942 do {
943 ulong e3 = exp3[i] + ((i == offset) ? extra : 0);
944 if (exp2[i] != e3)
945 {
946 if (exp2[i] > e3)
947 return 1;
948 else
949 return -1;
950 }
951 } while (--i >= 0);
952 return 0;
953 }
954
955 MPOLY_INLINE
mpoly_monomial_divides_tight(slong e1,slong e2,slong * prods,slong num)956 int mpoly_monomial_divides_tight(slong e1, slong e2, slong * prods, slong num)
957 {
958 slong j;
959
960 for (j = 0; j < num; j++)
961 {
962 slong d1 = (e1 % prods[j + 1])/prods[j];
963 slong d2 = (e2 % prods[j + 1])/prods[j];
964
965 if (d1 < d2)
966 return 0;
967 }
968
969 return 1;
970 }
971
972 MPOLY_INLINE
mpoly_max_degrees_tight(slong * max_exp,ulong * exps,slong len,slong * prods,slong num)973 void mpoly_max_degrees_tight(slong * max_exp,
974 ulong * exps, slong len, slong * prods, slong num)
975 {
976 slong i, j;
977
978 for (j = 0; j < num; j++)
979 max_exp[j] = 0;
980
981 for (i = 0; i < len; i++)
982 {
983 for (j = 0; j < num; j++)
984 {
985 slong d1 = (exps[i] % prods[j + 1])/prods[j];
986
987 if (d1 > max_exp[j])
988 max_exp[j] = d1;
989 }
990 }
991 }
992
993
994 /* ceiling(log_4(x)) - 1 */
mpoly_geobucket_clog4(slong x)995 MPOLY_INLINE slong mpoly_geobucket_clog4(slong x)
996 {
997 if (x <= 4)
998 return 0;
999 /*
1000 FLINT_BIT_COUNT returns unsigned int.
1001 Signed division is not defined.
1002 Do the calculation with unsigned ints and then convert to slong.
1003 */
1004 return (slong)((FLINT_BIT_COUNT(x - 1) - UWORD(1))/(UWORD(2)));
1005 }
1006
1007 /* single-limb packings ******************************************************/
1008
1009 MPOLY_INLINE
pack_exp2(ulong e0,ulong e1)1010 ulong pack_exp2(ulong e0, ulong e1)
1011 {
1012 return (e0 << (1*(FLINT_BITS/2))) +
1013 (e1 << (0*(FLINT_BITS/2)));
1014 }
1015
1016 MPOLY_INLINE
pack_exp3(ulong e0,ulong e1,ulong e2)1017 ulong pack_exp3(ulong e0, ulong e1, ulong e2)
1018 {
1019 return (e0 << (2*(FLINT_BITS/3))) +
1020 (e1 << (1*(FLINT_BITS/3))) +
1021 (e2 << (0*(FLINT_BITS/3)));
1022 }
1023
1024 MPOLY_INLINE
extract_exp(ulong e,int idx,int nvars)1025 ulong extract_exp(ulong e, int idx, int nvars)
1026 {
1027 return (e >> (idx*(FLINT_BITS/nvars))) &
1028 ((-UWORD(1)) >> (FLINT_BITS - FLINT_BITS/nvars));
1029 }
1030
1031 FLINT_DLL ulong _mpoly_bidegree(const ulong * Aexps, flint_bitcnt_t Abits,
1032 const mpoly_ctx_t mctx);
1033
1034 /* generators ****************************************************************/
1035
1036 FLINT_DLL void mpoly_gen_fields_ui(ulong * exp, slong var,
1037 const mpoly_ctx_t mctx);
1038
1039 FLINT_DLL void mpoly_gen_fields_fmpz(fmpz * exp, slong var,
1040 const mpoly_ctx_t mctx);
1041
1042 FLINT_DLL flint_bitcnt_t mpoly_gen_bits_required(slong var, const mpoly_ctx_t mctx);
1043
1044 /* return the index in the fields where the generator of index v is stored */
mpoly_gen_index(slong v,const mpoly_ctx_t mctx)1045 MPOLY_INLINE slong mpoly_gen_index(slong v, const mpoly_ctx_t mctx)
1046 {
1047 return mctx->rev ? v : mctx->nvars - 1 - v;
1048 }
1049
1050 FLINT_DLL void mpoly_gen_offset_shift_sp(slong * offset, slong * shift,
1051 slong var, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1052
1053 FLINT_DLL void mpoly_gen_monomial_offset_shift_sp(ulong * mexp, slong * offset,
1054 slong * shift, slong var, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1055
1056 FLINT_DLL void mpoly_gen_monomial_sp(ulong * oneexp, slong var,
1057 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1058
1059 FLINT_DLL slong mpoly_gen_offset_mp(slong var,
1060 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1061
1062 FLINT_DLL slong mpoly_gen_monomial_offset_mp(ulong * mexp, slong var,
1063 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1064
1065 FLINT_DLL void fmpz_mat_mul_vec(fmpz * v, const fmpz_mat_t M, fmpz * u);
1066
1067 FLINT_DLL void mpoly_compose_mat_gen(fmpz_mat_t M, const slong * c,
1068 const mpoly_ctx_t mctxB, const mpoly_ctx_t mctxAC);
1069
1070 FLINT_DLL void mpoly_compose_mat_fill_column(fmpz_mat_t M, const ulong * Cexp,
1071 flint_bitcnt_t Cbits, slong Bvar, const mpoly_ctx_t mctxB,
1072 const mpoly_ctx_t mctxAC);
1073
1074 /* Monomial arrays ***********************************************************/
1075
1076 FLINT_DLL void mpoly_get_cmpmask(ulong * cmpmask, slong N, flint_bitcnt_t bits,
1077 const mpoly_ctx_t mctx);
1078
1079 FLINT_DLL void mpoly_get_ovfmask(ulong * ovfmask, slong N, flint_bitcnt_t bits,
1080 const mpoly_ctx_t mctx);
1081
1082 FLINT_DLL int mpoly_monomials_cmp(const ulong * Aexps, flint_bitcnt_t Abits,
1083 const ulong * Bexps, flint_bitcnt_t Bbits,
1084 slong length, const mpoly_ctx_t mctx);
1085
1086 FLINT_DLL flint_bitcnt_t mpoly_exp_bits_required_ui(const ulong * user_exp,
1087 const mpoly_ctx_t mctx);
1088 FLINT_DLL flint_bitcnt_t mpoly_exp_bits_required_ffmpz(const fmpz * user_exp,
1089 const mpoly_ctx_t mctx);
1090 FLINT_DLL flint_bitcnt_t mpoly_exp_bits_required_pfmpz(fmpz * const * user_exp,
1091 const mpoly_ctx_t mctx);
1092
1093 MPOLY_INLINE
mpoly_gen_pow_exp_bits_required(slong v,ulong e,const mpoly_ctx_t mctx)1094 flint_bitcnt_t mpoly_gen_pow_exp_bits_required(slong v, ulong e,
1095 const mpoly_ctx_t mctx)
1096 {
1097 return 1 + FLINT_BIT_COUNT(e); /* only lex and deg supported */
1098 }
1099
1100 FLINT_DLL int mpoly_is_poly(const ulong * Aexps, slong Alen,
1101 flint_bitcnt_t Abits, slong var, const mpoly_ctx_t mctx);
1102
1103 FLINT_DLL void mpoly_pack_vec_ui(ulong * exp1, const ulong * exp2,
1104 flint_bitcnt_t bits, slong nfields, slong len);
1105 FLINT_DLL void mpoly_pack_vec_fmpz(ulong * exp1, const fmpz * exp2,
1106 flint_bitcnt_t bits, slong nfields, slong len);
1107
1108 FLINT_DLL void mpoly_unpack_vec_ui(ulong * exp1, const ulong * exp2,
1109 flint_bitcnt_t bits, slong nfields, slong len);
1110
1111 FLINT_DLL void mpoly_unpack_vec_fmpz(fmpz * exp1, const ulong * exp2,
1112 flint_bitcnt_t bits, slong nfields, slong len);
1113
1114 FLINT_DLL void mpoly_get_monomial_ui_unpacked_ffmpz(ulong * user_exps,
1115 const fmpz * poly_exps, const mpoly_ctx_t mctx);
1116
1117 FLINT_DLL void mpoly_get_monomial_ffmpz_unpacked_ffmpz(fmpz * user_exps,
1118 const fmpz * poly_exps, const mpoly_ctx_t mctx);
1119
1120 FLINT_DLL void mpoly_get_monomial_pfmpz_unpacked_ffmpz(fmpz ** user_exps,
1121 const fmpz * poly_exps, const mpoly_ctx_t mctx);
1122
1123 FLINT_DLL void mpoly_get_monomial_ui_unpacked_ui(ulong * user_exps,
1124 const ulong * poly_exps, const mpoly_ctx_t mctx);
1125
1126 FLINT_DLL void mpoly_get_monomial_ui_sp(ulong * user_exps,
1127 const ulong * poly_exps, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1128
1129 FLINT_DLL void mpoly_get_monomial_ui_mp(ulong * user_exps,
1130 const ulong * poly_exps, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1131
1132 FLINT_DLL void mpoly_get_monomial_si_mp(slong * user_exps,
1133 const ulong * poly_exps, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1134
mpoly_get_monomial_ui(ulong * user_exps,const ulong * poly_exps,flint_bitcnt_t bits,const mpoly_ctx_t mctx)1135 MPOLY_INLINE void mpoly_get_monomial_ui(ulong * user_exps,
1136 const ulong * poly_exps, flint_bitcnt_t bits, const mpoly_ctx_t mctx)
1137 {
1138 if (bits <= FLINT_BITS)
1139 mpoly_get_monomial_ui_sp(user_exps, poly_exps, bits, mctx);
1140 else
1141 mpoly_get_monomial_ui_mp(user_exps, poly_exps, bits, mctx);
1142 }
1143
mpoly_get_monomial_si(slong * user_exps,const ulong * poly_exps,flint_bitcnt_t bits,const mpoly_ctx_t mctx)1144 MPOLY_INLINE void mpoly_get_monomial_si(slong * user_exps,
1145 const ulong * poly_exps, flint_bitcnt_t bits, const mpoly_ctx_t mctx)
1146 {
1147 /* if bits <= FLINT_BITS and poly_exps is canonical, everything should be ok */
1148 if (bits <= FLINT_BITS)
1149 mpoly_get_monomial_ui_sp((ulong *) user_exps, poly_exps, bits, mctx);
1150 else
1151 mpoly_get_monomial_si_mp(user_exps, poly_exps, bits, mctx);
1152 }
1153
1154 FLINT_DLL ulong mpoly_get_monomial_var_exp_ui_sp(const ulong * poly_exps,
1155 slong var, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1156
1157 FLINT_DLL ulong mpoly_get_monomial_var_exp_ui_mp(const ulong * poly_exps,
1158 slong var, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1159
1160 FLINT_DLL slong mpoly_get_monomial_var_exp_si_mp(const ulong * poly_exps,
1161 slong var, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1162
mpoly_get_monomial_var_exp_ui(const ulong * poly_exps,slong var,flint_bitcnt_t bits,const mpoly_ctx_t mctx)1163 MPOLY_INLINE ulong mpoly_get_monomial_var_exp_ui(const ulong * poly_exps,
1164 slong var, flint_bitcnt_t bits, const mpoly_ctx_t mctx)
1165 {
1166 if (bits <= FLINT_BITS)
1167 return mpoly_get_monomial_var_exp_ui_sp(poly_exps, var, bits, mctx);
1168 else
1169 return mpoly_get_monomial_var_exp_ui_mp(poly_exps, var, bits, mctx);
1170 }
1171
mpoly_get_monomial_var_exp_si(const ulong * poly_exps,slong var,flint_bitcnt_t bits,const mpoly_ctx_t mctx)1172 MPOLY_INLINE slong mpoly_get_monomial_var_exp_si(const ulong * poly_exps,
1173 slong var, flint_bitcnt_t bits, const mpoly_ctx_t mctx)
1174 {
1175 if (bits <= FLINT_BITS)
1176 return (slong) mpoly_get_monomial_var_exp_ui_sp(poly_exps, var, bits, mctx);
1177 else
1178 return mpoly_get_monomial_var_exp_si_mp(poly_exps, var, bits, mctx);
1179 }
1180
1181 FLINT_DLL void mpoly_get_monomial_ffmpz(fmpz * exps, const ulong * poly_exps,
1182 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1183
1184 FLINT_DLL void mpoly_get_monomial_pfmpz(fmpz ** exps, const ulong * poly_exps,
1185 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1186
1187 FLINT_DLL void mpoly_set_monomial_ui(ulong * exp1, const ulong * exp2,
1188 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1189 FLINT_DLL void mpoly_set_monomial_ffmpz(ulong * exp1, const fmpz * exp2,
1190 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1191 FLINT_DLL void mpoly_set_monomial_pfmpz(ulong * exp1, fmpz * const * exp2,
1192 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1193
1194 FLINT_DLL int mpoly_repack_monomials(ulong * exps1, flint_bitcnt_t bits1,
1195 const ulong * exps2, flint_bitcnt_t bits2, slong len,
1196 const mpoly_ctx_t mctx);
1197
1198 FLINT_DLL void mpoly_pack_monomials_tight(ulong * exp1,
1199 const ulong * exp2, slong len, const slong * mults,
1200 slong num, slong bits);
1201
1202 FLINT_DLL void mpoly_unpack_monomials_tight(ulong * e1, ulong * e2, slong len,
1203 slong * mults, slong num, slong bits);
1204
1205 FLINT_DLL int mpoly_monomial_exists(slong * index, const ulong * poly_exps,
1206 const ulong * exp, slong len, slong N, const ulong * cmpmask);
1207
1208 FLINT_DLL slong mpoly_monomial_index1_nomask(ulong * Aexps, slong Alen, ulong e);
1209
1210 FLINT_DLL slong mpoly_monomial_index_ui(const ulong * Aexp, flint_bitcnt_t Abits,
1211 slong Alength, const ulong * exp, const mpoly_ctx_t mctx);
1212
1213 FLINT_DLL slong mpoly_monomial_index_pfmpz(const ulong * Aexp, flint_bitcnt_t Abits,
1214 slong Alength, fmpz * const * exp, const mpoly_ctx_t mctx);
1215
1216 FLINT_DLL slong mpoly_monomial_index_monomial(const ulong * Aexp,
1217 flint_bitcnt_t Abits, slong Alength, const ulong * Mexp,
1218 flint_bitcnt_t Mbits, const mpoly_ctx_t mctx);
1219
1220 FLINT_DLL void mpoly_min_fields_ui_sp(ulong * min_fields, const ulong * poly_exps,
1221 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1222
1223 FLINT_DLL void mpoly_min_fields_fmpz(fmpz * min_fields, const ulong * poly_exps,
1224 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1225
1226 FLINT_DLL void mpoly_max_fields_ui_sp(ulong * max_fields, const ulong * poly_exps,
1227 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1228
1229 FLINT_DLL void mpoly_max_fields_fmpz(fmpz * max_fields, const ulong * poly_exps,
1230 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1231
1232 FLINT_DLL int mpoly_degrees_fit_si(const ulong * poly_exps,
1233 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1234
1235 FLINT_DLL void mpoly_degrees_si(slong * user_degs, const ulong * poly_exps,
1236 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1237
1238 FLINT_DLL void mpoly_degrees_si_threaded(slong * user_degs, const ulong * poly_exps,
1239 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx,
1240 const thread_pool_handle * handles, slong num_handles);
1241
1242 FLINT_DLL void mpoly_degrees_ffmpz(fmpz * user_degs, const ulong * poly_exps,
1243 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1244
1245 FLINT_DLL void mpoly_degrees_pfmpz(fmpz ** user_degs, const ulong * poly_exps,
1246 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1247
1248 FLINT_DLL slong mpoly_degree_si(const ulong * poly_exps,
1249 slong len, flint_bitcnt_t bits, slong var, const mpoly_ctx_t mctx);
1250
1251 FLINT_DLL void mpoly_degree_fmpz(fmpz_t deg, const ulong * poly_exps,
1252 slong len, flint_bitcnt_t bits, slong var, const mpoly_ctx_t mctx);
1253
1254 FLINT_DLL int mpoly_total_degree_fits_si(const ulong * exps,
1255 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1256
1257 FLINT_DLL slong mpoly_total_degree_si(const ulong * exps,
1258 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1259
1260 FLINT_DLL void mpoly_total_degree_fmpz(fmpz_t totdeg, const ulong * exps,
1261 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1262
1263 FLINT_DLL void mpoly_total_degree_fmpz_ref(fmpz_t totdeg, const ulong * exps,
1264 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1265
1266 FLINT_DLL void mpoly_used_vars_or(int * used, const ulong * exps,
1267 slong len, flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1268
1269 FLINT_DLL int mpoly_monomial_cmp_general(ulong * Aexp, flint_bitcnt_t Abits,
1270 ulong * Bexp, flint_bitcnt_t Bbits, const mpoly_ctx_t mctx);
1271
1272 FLINT_DLL void mpoly_search_monomials(
1273 slong ** e_ind, ulong * e, slong * e_score,
1274 slong * t1, slong * t2, slong *t3,
1275 slong lower, slong upper,
1276 const ulong * a, slong a_len, const ulong * b, slong b_len,
1277 slong N, const ulong * cmpmask);
1278
1279 FLINT_DLL void mpoly_main_variable_split_LEX(slong * ind, ulong * pexp,
1280 const ulong * Aexp, slong l1, slong Alen,
1281 const ulong * mults, slong num, slong Abits);
1282
1283 FLINT_DLL void mpoly_main_variable_split_DEG(slong * ind, ulong * pexp,
1284 const ulong * Aexp, slong l1, slong Alen,
1285 ulong deg, slong num, slong Abits);
1286
1287 FLINT_DLL int mpoly_term_exp_fits_si(ulong * exps, flint_bitcnt_t bits,
1288 slong n, const mpoly_ctx_t mctx);
1289
1290 FLINT_DLL int mpoly_term_exp_fits_ui(ulong * exps, flint_bitcnt_t bits,
1291 slong n, const mpoly_ctx_t mctx);
1292
1293 FLINT_DLL int mpoly_is_gen(ulong * exps, slong var,
1294 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1295
1296 FLINT_DLL int mpoly_monomials_valid_test(ulong * exps, slong len,
1297 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1298
1299 FLINT_DLL int mpoly_monomials_overflow_test(ulong * exps, slong len,
1300 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1301
1302 FLINT_DLL int mpoly_monomials_inorder_test(ulong * exps, slong len,
1303 flint_bitcnt_t bits, const mpoly_ctx_t mctx);
1304
1305 FLINT_DLL void mpoly_reverse(ulong * Aexp, const ulong * Bexp, slong len, slong N);
1306
1307 FLINT_DLL void mpoly_monomials_deflation(fmpz * shift, fmpz * stride,
1308 const ulong * Aexps, flint_bitcnt_t Abits, slong Alength,
1309 const mpoly_ctx_t mctx);
1310
1311 FLINT_DLL void mpoly_monomials_deflate(ulong * Aexps, flint_bitcnt_t Abits,
1312 const ulong * Bexps, flint_bitcnt_t Bbits, slong Blength,
1313 const fmpz * shift, const fmpz * stride, const mpoly_ctx_t mctx);
1314
1315 FLINT_DLL void mpoly_monomials_inflate(ulong * Aexps, flint_bitcnt_t Abits,
1316 const ulong * Bexps, flint_bitcnt_t Bbits, slong Blength,
1317 const fmpz * shift, const fmpz * stride, const mpoly_ctx_t mctx);
1318
1319 FLINT_DLL void _mpoly_gen_shift_right(ulong * Aexp, flint_bitcnt_t Abits,
1320 slong Alength, slong var, ulong amount, const mpoly_ctx_t mctx);
1321
1322 FLINT_DLL void _mpoly_gen_shift_right_fmpz(ulong * Aexp, flint_bitcnt_t Abits,
1323 slong Alength, slong var, const fmpz_t amount, const mpoly_ctx_t mctx);
1324
1325 FLINT_DLL void _mpoly_gen_shift_left(ulong * Aexp, flint_bitcnt_t Abits,
1326 slong Alength, slong var, ulong amount, const mpoly_ctx_t mctx);
1327
1328 FLINT_DLL void mpoly_monomials_shift_right_ui(ulong * Aexps, flint_bitcnt_t Abits,
1329 slong Alength, const ulong * user_exps, const mpoly_ctx_t mctx);
1330
1331 FLINT_DLL void mpoly_monomials_shift_right_ffmpz(ulong * Aexps, flint_bitcnt_t Abits,
1332 slong Alength, const fmpz * user_exps, const mpoly_ctx_t mctx);
1333
1334 FLINT_DLL void mpoly1_fill_marks(ulong ** Dcoeffs, slong * Dlen, slong * Dalloc,
1335 const ulong * Aexps, slong Alen, flint_bitcnt_t Abits,
1336 const mpoly_ctx_t mctx);
1337
1338 FLINT_DLL void mpoly2_fill_marks(ulong ** Dcoeffs, slong * Dlen, slong * Dalloc,
1339 const ulong * Aexps, slong Alen, flint_bitcnt_t Abits,
1340 const mpoly_ctx_t mctx);
1341
1342 FLINT_DLL void mpoly_to_mpolyl_perm_deflate(
1343 ulong * Aexps,
1344 flint_bitcnt_t Abits,
1345 const mpoly_ctx_t Actx,
1346 ulong * Bexps,
1347 flint_bitcnt_t Bbits,
1348 const mpoly_ctx_t Bctx,
1349 slong length,
1350 const slong * perm,
1351 const ulong * shift,
1352 const ulong * stride);
1353
1354 FLINT_DLL void mpoly_from_mpolyl_perm_inflate(
1355 ulong * Bexps,
1356 flint_bitcnt_t Bbits,
1357 const mpoly_ctx_t Bctx,
1358 ulong * Aexps,
1359 flint_bitcnt_t Abits,
1360 const mpoly_ctx_t Actx,
1361 slong length,
1362 const slong * perm,
1363 const ulong * shift,
1364 const ulong * stride);
1365
1366 /* gcd ***********************************************************************/
1367
1368 #define MPOLY_GCD_USE_HENSEL 1
1369 #define MPOLY_GCD_USE_BROWN 2
1370 #define MPOLY_GCD_USE_ZIPPEL 4
1371 #define MPOLY_GCD_USE_ZIPPEL2 8
1372 #define MPOLY_GCD_USE_PRS 16
1373 #define MPOLY_GCD_USE_ALL 31
1374
1375 typedef struct
1376 {
1377 ulong * Amax_exp;
1378 ulong * Amin_exp;
1379 ulong * Astride;
1380 slong * Adeflate_deg;
1381 slong * Alead_count;
1382 slong * Atail_count;
1383
1384 ulong * Bmax_exp;
1385 ulong * Bmin_exp;
1386 ulong * Bstride;
1387 slong * Bdeflate_deg;
1388 slong * Blead_count;
1389 slong * Btail_count;
1390
1391 ulong * Gmin_exp;
1392 ulong * Abarmin_exp;
1393 ulong * Bbarmin_exp;
1394 ulong * Gstride;
1395 slong * Gterm_count_est;
1396 slong * Gdeflate_deg_bound;
1397
1398 flint_bitcnt_t Gbits, Abarbits, Bbarbits;
1399
1400 slong mvars;
1401 slong Adeflate_tdeg;
1402 slong Bdeflate_tdeg;
1403
1404 double Adensity;
1405 double Bdensity;
1406
1407 double hensel_time, brown_time, zippel_time, zippel2_time;
1408 slong * hensel_perm, * brown_perm, * zippel_perm, * zippel2_perm;
1409 unsigned int can_use;
1410 int Gdeflate_deg_bounds_are_nice; /* all of Gdeflate_deg_bound came from real gcd computations */
1411
1412 char * data;
1413 } mpoly_gcd_info_struct;
1414
1415 typedef mpoly_gcd_info_struct mpoly_gcd_info_t[1];
1416
1417 FLINT_DLL void mpoly_gcd_info_init(mpoly_gcd_info_t I, slong nvars);
1418
1419 FLINT_DLL void mpoly_gcd_info_clear(mpoly_gcd_info_t I);
1420
1421 FLINT_DLL void mpoly_gcd_info_limits(ulong * Amax_exp, ulong * Amin_exp,
1422 slong * Amax_exp_count, slong * Amin_exp_count,
1423 const ulong * Aexps, flint_bitcnt_t Abits, slong Alength,
1424 const mpoly_ctx_t mctx);
1425 FLINT_DLL void mpoly_gcd_info_stride(ulong * strides,
1426 const ulong * Aexps, flint_bitcnt_t Abits, slong Alength,
1427 const ulong * Amax_exp, const ulong * Amin_exp,
1428 const ulong * Bexps, flint_bitcnt_t Bbits, slong Blength,
1429 const ulong * Bmax_exp, const ulong * Bmin_exp,
1430 const mpoly_ctx_t mctx);
1431
1432 FLINT_DLL void mpoly_gcd_info_set_perm(mpoly_gcd_info_t I,
1433 slong Alength, slong Blength, const mpoly_ctx_t mctx);
1434
1435 FLINT_DLL slong mpoly_gcd_info_get_brown_upper_limit(const mpoly_gcd_info_t I,
1436 slong var, slong bound);
1437
1438 FLINT_DLL void mpoly_gcd_info_measure_hensel(mpoly_gcd_info_t I,
1439 slong Alength, slong Blength, const mpoly_ctx_t mctx);
1440
1441 FLINT_DLL void mpoly_gcd_info_measure_brown(mpoly_gcd_info_t I,
1442 slong Alength, slong Blength, const mpoly_ctx_t mctx);
1443
1444 FLINT_DLL void mpoly_gcd_info_measure_bma(mpoly_gcd_info_t I,
1445 slong Alength, slong Blength, const mpoly_ctx_t mctx);
1446
1447 FLINT_DLL void mpoly_gcd_info_measure_zippel(mpoly_gcd_info_t I,
1448 slong Alength, slong Blength, const mpoly_ctx_t mctx);
1449
1450 FLINT_DLL void mpoly_gcd_info_measure_zippel2(mpoly_gcd_info_t I,
1451 slong Alength, slong Blength, const mpoly_ctx_t mctx);
1452
1453 FLINT_DLL int mpoly_monomial_cofactors(fmpz * Abarexps, fmpz * Bbarexps,
1454 const ulong * Aexps, flint_bitcnt_t Abits,
1455 const ulong * Bexps, flint_bitcnt_t Bbits,
1456 slong length, const mpoly_ctx_t mctx);
1457
1458 /* factoring ****************************************************************/
1459
1460 #define MPOLY_FACTOR_USE_ZAS 1
1461 #define MPOLY_FACTOR_USE_WANG 2
1462 #define MPOLY_FACTOR_USE_ZIP 4
1463 #define MPOLY_FACTOR_USE_ALL 7
1464
1465 FLINT_DLL int mpoly_is_proved_not_square(const ulong * Aexps,
1466 slong Alen, flint_bitcnt_t Abits, slong N, ulong * t);
1467
1468 FLINT_DLL void mpoly_remove_var_powers(fmpz * var_powers, ulong * Aexps,
1469 flint_bitcnt_t Abits, slong Alen, const mpoly_ctx_t mctx);
1470
1471 FLINT_DLL slong _mpoly_compress_exps(slong * V, slong * D, slong * deg,
1472 slong * S, slong n, slong l);
1473
1474 FLINT_DLL int mpoly_test_irreducible(ulong * Aexps, flint_bitcnt_t Abits,
1475 slong Alen, const mpoly_ctx_t ctx);
1476
1477 FLINT_DLL int _mpoly_test_irreducible(slong * Aexps, slong stride, slong Alen,
1478 slong nvars, flint_rand_t state, slong tries_left);
1479
1480 typedef struct {
1481 slong mvars;
1482 slong nvars;
1483 slong * exps;
1484 slong exps_alloc;
1485 slong * rest;
1486 slong rest_alloc;
1487 slong * umat;
1488 slong * deltas;
1489 slong * degs;
1490 int is_trivial;
1491 int is_perm;
1492 int is_irred;
1493 } mpoly_compression_struct;
1494
1495 typedef mpoly_compression_struct mpoly_compression_t[1];
1496
1497 FLINT_DLL void mpoly_compression_init(mpoly_compression_t M);
1498
1499 FLINT_DLL void mpoly_compression_clear(mpoly_compression_t M);
1500
1501 FLINT_DLL void mpoly_compression_set(mpoly_compression_t M, const ulong * Aexps,
1502 flint_bitcnt_t Abits, slong Alen, const mpoly_ctx_t mctx);
1503
1504 FLINT_DLL void mpoly_bivar_cld_bounds(slong * l, slong n);
1505
1506 MPOLY_INLINE
_slong_array_fit_length(slong ** array,slong * alloc,slong len)1507 void _slong_array_fit_length(slong ** array, slong * alloc, slong len)
1508 {
1509 if (len <= *alloc)
1510 return;
1511 len = FLINT_MAX(len, *alloc + *alloc/4 + 1);
1512 *alloc = len;
1513 *array = (slong *) flint_realloc(*array, len*sizeof(slong));
1514 }
1515
1516 /* Heap **********************************************************************/
1517
1518 #define HEAP_LEFT(i) (2*(i))
1519 #define HEAP_RIGHT(i) (2*(i) + 1)
1520 #define HEAP_PARENT(i) ((i)/2)
1521
1522 #define HEAP_ASSIGN(h, c1, c2) \
1523 do { \
1524 (h).exp = (c1); \
1525 (h).next = (c2); \
1526 } while (0)
1527
1528 MPOLY_INLINE
_mpoly_heap_pop1(mpoly_heap1_s * heap,slong * heap_len,ulong maskhi)1529 void * _mpoly_heap_pop1(mpoly_heap1_s * heap, slong * heap_len, ulong maskhi)
1530 {
1531 ulong exp;
1532 slong i, j, s = --(*heap_len);
1533 void * x = heap[1].next;
1534
1535 i = 1;
1536 j = 2;
1537
1538 while (j < s)
1539 {
1540 if ((heap[j].exp^maskhi) <= (heap[j + 1].exp^maskhi))
1541 j++;
1542 heap[i] = heap[j];
1543 i = j;
1544 j = HEAP_LEFT(j);
1545 }
1546
1547 /* insert last element into heap[i] */
1548 exp = heap[s].exp;
1549 j = HEAP_PARENT(i);
1550
1551 while (i > 1 && (exp^maskhi) > (heap[j].exp^maskhi))
1552 {
1553 heap[i] = heap[j];
1554 i = j;
1555 j = HEAP_PARENT(j);
1556 }
1557
1558 heap[i] = heap[s];
1559
1560 return x;
1561 }
1562
1563 MPOLY_INLINE
_mpoly_heap_insert1(mpoly_heap1_s * heap,ulong exp,void * x,slong * next_loc,slong * heap_len,ulong maskhi)1564 void _mpoly_heap_insert1(mpoly_heap1_s * heap, ulong exp, void * x,
1565 slong * next_loc, slong * heap_len, ulong maskhi)
1566 {
1567 slong i = *heap_len, j, n = *heap_len;
1568
1569 if (i != 1 && exp == heap[1].exp)
1570 {
1571 ((mpoly_heap_t *) x)->next = (mpoly_heap_t *) heap[1].next;
1572 heap[1].next = x;
1573
1574 return;
1575 }
1576
1577 if (*next_loc < *heap_len)
1578 {
1579 if (exp == heap[*next_loc].exp)
1580 {
1581 ((mpoly_heap_t *) x)->next = (mpoly_heap_t *) heap[*next_loc].next;
1582 heap[*next_loc].next = x;
1583 return;
1584 }
1585 }
1586
1587
1588
1589 while ((j = HEAP_PARENT(i)) >= 1)
1590 {
1591 if (exp == heap[j].exp)
1592 {
1593 ((mpoly_heap_t *) x)->next = (mpoly_heap_t *) heap[j].next;
1594 heap[j].next = x;
1595
1596 *next_loc = j;
1597
1598
1599 return;
1600 } else if ((exp^maskhi) > (heap[j].exp^maskhi))
1601 i = j;
1602 else
1603 break;
1604 }
1605
1606 (*heap_len)++;
1607
1608 while (n > i)
1609 {
1610 heap[n] = heap[HEAP_PARENT(n)];
1611 n = HEAP_PARENT(n);
1612 }
1613
1614 HEAP_ASSIGN(heap[i], exp, x);
1615 }
1616
1617 MPOLY_INLINE
_mpoly_heap_pop(mpoly_heap_s * heap,slong * heap_len,slong N,const ulong * cmpmask)1618 void * _mpoly_heap_pop(mpoly_heap_s * heap, slong * heap_len, slong N,
1619 const ulong * cmpmask)
1620 {
1621 ulong * exp;
1622 slong i, j, s = --(*heap_len);
1623 mpoly_heap_t * x = (mpoly_heap_t *) heap[1].next;
1624
1625 i = 1;
1626 j = 2;
1627
1628 while (j < s)
1629 {
1630 if (!mpoly_monomial_gt(heap[j].exp, heap[j + 1].exp, N, cmpmask))
1631 j++;
1632 heap[i] = heap[j];
1633 i = j;
1634 j = HEAP_LEFT(j);
1635 }
1636
1637 /* insert last element into heap[i] */
1638 exp = heap[s].exp;
1639 j = HEAP_PARENT(i);
1640
1641 while (i > 1 && mpoly_monomial_gt(exp, heap[j].exp, N, cmpmask))
1642 {
1643 heap[i] = heap[j];
1644 i = j;
1645 j = HEAP_PARENT(j);
1646 }
1647
1648 heap[i] = heap[s];
1649
1650 return x;
1651 }
1652
1653 MPOLY_INLINE
_mpoly_heap_insert(mpoly_heap_s * heap,ulong * exp,void * x,slong * next_loc,slong * heap_len,slong N,const ulong * cmpmask)1654 int _mpoly_heap_insert(mpoly_heap_s * heap, ulong * exp, void * x,
1655 slong * next_loc, slong * heap_len, slong N, const ulong * cmpmask)
1656 {
1657 slong i = *heap_len, j, n = *heap_len;
1658
1659 if (i != 1 && mpoly_monomial_equal(exp, heap[1].exp, N))
1660 {
1661 ((mpoly_heap_t *) x)->next = (mpoly_heap_t *) heap[1].next;
1662 heap[1].next = x;
1663
1664 return 0;
1665 }
1666
1667 if (*next_loc < *heap_len)
1668 {
1669 if (mpoly_monomial_equal(exp, heap[*next_loc].exp, N))
1670 {
1671 ((mpoly_heap_t *) x)->next = (mpoly_heap_t *) heap[*next_loc].next;
1672 heap[*next_loc].next = x;
1673 return 0;
1674 }
1675 }
1676
1677 while ((j = HEAP_PARENT(i)) >= 1)
1678 {
1679 if (!mpoly_monomial_gt(exp, heap[j].exp, N, cmpmask))
1680 break;
1681
1682 i = j;
1683 }
1684
1685 if (j >= 1 && mpoly_monomial_equal(exp, heap[j].exp, N))
1686 {
1687 ((mpoly_heap_t *) x)->next = (mpoly_heap_t *) heap[j].next;
1688 heap[j].next = x;
1689 *next_loc = j;
1690
1691 return 0;
1692 }
1693
1694 (*heap_len)++;
1695
1696 while (n > i)
1697 {
1698 heap[n] = heap[HEAP_PARENT(n)];
1699 n = HEAP_PARENT(n);
1700 }
1701
1702 HEAP_ASSIGN(heap[i], exp, x);
1703
1704 return 1;
1705 }
1706
1707 /* generic parsing ***********************************************************/
1708
1709 typedef struct {
1710 char * coeffs;
1711 fmpz * exps;
1712 slong length;
1713 slong alloc;
1714 } mpoly_univar_struct;
1715
1716 typedef mpoly_univar_struct mpoly_univar_t[1];
1717
1718 typedef struct {
1719 slong elem_size;
1720 const void * ctx;
1721 void (*init)(void *, const void *);
1722 void (*clear)(void *, const void *);
1723 int (*is_zero)(const void *, const void *);
1724 void (*zero)(void *, const void *);
1725 void (*one)(void *, const void *);
1726 void (*set_fmpz)(void *, const fmpz_t, const void *);
1727 void (*set)(void *, const void *, const void *);
1728 void (*swap)(void *, void *, const void *);
1729 void (*neg)(void *, const void *, const void *);
1730 void (*add)(void *, const void *, const void *, const void *);
1731 void (*sub)(void *, const void *, const void *, const void *);
1732 void (*mul_fmpz)(void *, const void *, const fmpz_t, const void *);
1733 void (*mul)(void *, const void *, const void *, const void *);
1734 void (*divexact)(void *, const void *, const void *, const void *);
1735 int (*divides)(void *, const void *, const void *, const void *);
1736 int (*pow_fmpz)(void *, const void *, const fmpz_t, const void *);
1737 slong (*length)(const void *, const void *);
1738 } mpoly_void_ring_t[1];
1739
1740 FLINT_DLL void * mpoly_void_ring_elem_init(mpoly_void_ring_t R);
1741
1742 FLINT_DLL void mpoly_void_ring_elem_clear(void * a, mpoly_void_ring_t R);
1743
1744 FLINT_DLL void mpoly_univar_init(mpoly_univar_t A, mpoly_void_ring_t R);
1745
1746 FLINT_DLL void mpoly_univar_clear(mpoly_univar_t A, mpoly_void_ring_t R);
1747
1748 FLINT_DLL void mpoly_univar_swap(mpoly_univar_t A, mpoly_univar_t B);
1749
1750 FLINT_DLL void mpoly_univar_fit_length(mpoly_univar_t A, slong len,
1751 mpoly_void_ring_t R);
1752
1753 FLINT_DLL void mpoly_univar_init2(mpoly_univar_t A, slong len,
1754 mpoly_void_ring_t R);
1755
1756 FLINT_DLL int mpoly_univar_pseudo_gcd_ducos(mpoly_univar_t G,
1757 mpoly_univar_t B, mpoly_univar_t A, mpoly_void_ring_t R);
1758
1759 FLINT_DLL int mpoly_univar_resultant(void * r, mpoly_univar_t fx,
1760 mpoly_univar_t gx, mpoly_void_ring_t R);
1761
1762 FLINT_DLL int mpoly_univar_discriminant(void * d, mpoly_univar_t fx,
1763 mpoly_void_ring_t R);
1764
1765 typedef struct {
1766 char * str;
1767 slong str_len;
1768 } string_with_length_struct;
1769
1770 typedef struct {
1771 mpoly_void_ring_t R;
1772 slong * stack;
1773 slong stack_len;
1774 slong stack_alloc;
1775 char * estore;
1776 slong estore_len;
1777 slong estore_alloc;
1778 void * tmp;
1779 string_with_length_struct * terminal_strings;
1780 char * terminal_values;
1781 slong terminals_alloc;
1782 slong terminals_len;
1783 } mpoly_parse_struct;
1784
1785 typedef mpoly_parse_struct mpoly_parse_t[1];
1786
1787 FLINT_DLL void mpoly_parse_init(mpoly_parse_t E);
1788
1789 FLINT_DLL void mpoly_parse_clear(mpoly_parse_t E);
1790
1791 FLINT_DLL void mpoly_parse_add_terminal(mpoly_parse_t E,
1792 const char * s, const void * v);
1793
1794 FLINT_DLL int mpoly_parse_parse(mpoly_parse_t E, void * res,
1795 const char * s, slong len);
1796
1797 /* chunking */
1798
1799 /*
1800 Set i1[i] to the index of the i-th "coefficient" in variable k of num
1801 variables, each taking the given number of bits in the exponent. This
1802 assumes there are l1 "coefficients" in a list of len1 exponents.
1803 Note this doesn't currently mask the relevant bits.
1804 */
1805 MPOLY_INLINE
mpoly_main_variable_terms1(slong * i1,slong * n1,const ulong * exp1,slong l1,slong len1,slong k,slong num,slong bits)1806 void mpoly_main_variable_terms1(slong * i1, slong * n1, const ulong * exp1,
1807 slong l1, slong len1, slong k, slong num, slong bits)
1808 {
1809 slong i, j = 0;
1810 slong shift = bits*(k - 1);
1811
1812 i1[0] = 0;
1813 for (i = 0; i < l1 - 1; i++)
1814 {
1815 while (j < len1 && (l1 - i - 1) == (slong) (exp1[j] >> shift))
1816 j++;
1817
1818 i1[i + 1] = j;
1819 n1[i] = j - i1[i];
1820 }
1821 n1[l1 - 1] = len1 - j;
1822 }
1823
1824 #ifdef __cplusplus
1825 }
1826 #endif
1827
1828 #endif
1829
1830