My Project
Loading...
Searching...
No Matches
facMul.cc
Go to the documentation of this file.
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facMul.cc
5 *
6 * This file implements functions for fast multiplication and division with
7 * remainder.
8 *
9 * Nomenclature rules: kronSub* -> plain Kronecker substitution
10 * reverseSubst* -> reverse Kronecker substitution
11 * kronSubRecipro* -> reciprocal Kronecker substitution as
12 * described in D. Harvey "Faster
13 * polynomial multiplication via
14 * multipoint Kronecker substitution"
15 *
16 * @author Martin Lee
17 *
18 **/
19/*****************************************************************************/
20
21#include "debug.h"
22
23#include "config.h"
24
25#include <math.h>
26
27#include "canonicalform.h"
28#include "facMul.h"
29#include "cf_util.h"
30#include "cf_iter.h"
31#include "cf_algorithm.h"
33
34#ifdef HAVE_NTL
35#include <NTL/lzz_pEX.h>
36#include "NTLconvert.h"
37#endif
38
39#ifdef HAVE_FLINT
40#include "FLINTconvert.h"
41#include "flint/fq_nmod_vec.h"
42#if __FLINT_RELEASE >= 20503
43#include "flint/mpn_extras.h"
44#include "flint/fmpz_mod.h"
45#endif
46#endif
47
48// univariate polys
49#if defined(HAVE_NTL) || defined(HAVE_FLINT)
50
51#ifdef HAVE_FLINT
52void kronSubQa (fmpz_poly_t result, const CanonicalForm& A, int d)
53{
54 int degAy= degree (A);
55 fmpz_poly_init2 (result, d*(degAy + 1));
56 _fmpz_poly_set_length (result, d*(degAy + 1));
58 for (CFIterator i= A; i.hasTerms(); i++)
59 {
60 if (i.coeff().inBaseDomain())
61 convertCF2initFmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
62 else
63 for (j= i.coeff(); j.hasTerms(); j++)
64 convertCF2initFmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
65 j.coeff());
66 }
67 _fmpz_poly_normalise(result);
68}
69
70
72reverseSubstQa (const fmpz_poly_t F, int d, const Variable& x,
73 const Variable& alpha, const CanonicalForm& den)
74{
76 int i= 0;
77 int degf= fmpz_poly_degree (F);
78 int k= 0;
79 int degfSubK;
80 int repLength;
81 fmpq_poly_t buf;
82 fmpq_poly_t mipo;
84 while (degf >= k)
85 {
86 degfSubK= degf - k;
87 if (degfSubK >= d)
88 repLength= d;
89 else
90 repLength= degfSubK + 1;
91
92 fmpq_poly_init2 (buf, repLength);
93 _fmpq_poly_set_length (buf, repLength);
94 _fmpz_vec_set (buf->coeffs, F->coeffs + k, repLength);
95 _fmpq_poly_normalise (buf);
96 fmpq_poly_rem (buf, buf, mipo);
97
99 fmpq_poly_clear (buf);
100 i++;
101 k= d*i;
102 }
103 fmpq_poly_clear (mipo);
104 result /= den;
105 return result;
106}
107
110 const Variable& alpha)
111{
112 CanonicalForm A= F;
114
115 CanonicalForm denA= bCommonDen (A);
116 CanonicalForm denB= bCommonDen (B);
117
118 A *= denA;
119 B *= denB;
120 int degAa= degree (A, alpha);
121 int degBa= degree (B, alpha);
122 int d= degAa + 1 + degBa;
123
124 fmpz_poly_t FLINTA,FLINTB;
125 kronSubQa (FLINTA, A, d);
126 kronSubQa (FLINTB, B, d);
127
128 fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
129
130 denA *= denB;
131 A= reverseSubstQa (FLINTA, d, F.mvar(), alpha, denA);
132
133 fmpz_poly_clear (FLINTA);
134 fmpz_poly_clear (FLINTB);
135 return A;
136}
137
140{
141 CanonicalForm A= F;
143
144 CanonicalForm denA= bCommonDen (A);
145 CanonicalForm denB= bCommonDen (B);
146
147 A *= denA;
148 B *= denB;
149 fmpz_poly_t FLINTA,FLINTB;
150 convertFacCF2Fmpz_poly_t (FLINTA, A);
151 convertFacCF2Fmpz_poly_t (FLINTB, B);
152 fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
153 denA *= denB;
154 A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
155 A /= denA;
156 fmpz_poly_clear (FLINTA);
157 fmpz_poly_clear (FLINTB);
158
159 return A;
160}
161
162/*CanonicalForm
163mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
164{
165 CanonicalForm A= F;
166 CanonicalForm B= G;
167
168 fmpq_poly_t FLINTA,FLINTB;
169 convertFacCF2Fmpq_poly_t (FLINTA, A);
170 convertFacCF2Fmpq_poly_t (FLINTB, B);
171
172 fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
173 A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
174 fmpq_poly_clear (FLINTA);
175 fmpq_poly_clear (FLINTB);
176 return A;
177}*/
178
181{
182 CanonicalForm A= F;
184
185 fmpq_poly_t FLINTA,FLINTB;
186 convertFacCF2Fmpq_poly_t (FLINTA, A);
187 convertFacCF2Fmpq_poly_t (FLINTB, B);
188
189 fmpq_poly_div (FLINTA, FLINTA, FLINTB);
190 A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
191
192 fmpq_poly_clear (FLINTA);
193 fmpq_poly_clear (FLINTB);
194 return A;
195}
196
199{
200 CanonicalForm A= F;
202
203 fmpq_poly_t FLINTA,FLINTB;
204 convertFacCF2Fmpq_poly_t (FLINTA, A);
205 convertFacCF2Fmpq_poly_t (FLINTB, B);
206
207 fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
208 A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
209
210 fmpq_poly_clear (FLINTA);
211 fmpq_poly_clear (FLINTB);
212 return A;
213}
214
217 const Variable& alpha, int m)
218{
219 CanonicalForm A= F;
221
222 CanonicalForm denA= bCommonDen (A);
223 CanonicalForm denB= bCommonDen (B);
224
225 A *= denA;
226 B *= denB;
227
228 int degAa= degree (A, alpha);
229 int degBa= degree (B, alpha);
230 int d= degAa + 1 + degBa;
231
232 fmpz_poly_t FLINTA,FLINTB;
233 kronSubQa (FLINTA, A, d);
234 kronSubQa (FLINTB, B, d);
235
236 int k= d*m;
237 fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
238
239 denA *= denB;
240 A= reverseSubstQa (FLINTA, d, F.mvar(), alpha, denA);
241 fmpz_poly_clear (FLINTA);
242 fmpz_poly_clear (FLINTB);
243 return A;
244}
245
248{
249 if (F.inCoeffDomain() && G.inCoeffDomain())
250 return F*G;
251 if (F.inCoeffDomain())
252 return mod (F*G, power (G.mvar(), m));
253 if (G.inCoeffDomain())
254 return mod (F*G, power (F.mvar(), m));
257 return mulFLINTQaTrunc (F, G, alpha, m);
258
259 CanonicalForm A= F;
261
262 CanonicalForm denA= bCommonDen (A);
263 CanonicalForm denB= bCommonDen (B);
264
265 A *= denA;
266 B *= denB;
267 fmpz_poly_t FLINTA,FLINTB;
268 convertFacCF2Fmpz_poly_t (FLINTA, A);
269 convertFacCF2Fmpz_poly_t (FLINTB, B);
270 fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
271 denA *= denB;
272 A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
273 A /= denA;
274 fmpz_poly_clear (FLINTA);
275 fmpz_poly_clear (FLINTB);
276
277 return A;
278}
279
281{
282 if (d == 0)
283 return F;
284 if (F.inCoeffDomain())
285 return F*power (x,d);
287 CFIterator i= F;
288 while (d - i.exp() < 0)
289 i++;
290
291 for (; i.hasTerms() && (d - i.exp() >= 0); i++)
292 result += i.coeff()*power (x, d - i.exp());
293 return result;
294}
295
297newtonInverse (const CanonicalForm& F, const int n, const Variable& x)
298{
299 int l= ilog2(n);
300
302 if (F.inCoeffDomain())
303 g= F;
304 else
305 g= F [0];
306
307 if (!F.inCoeffDomain())
308 ASSERT (F.mvar() == x, "main variable of F and x differ");
309 ASSERT (!g.isZero(), "expected a unit");
310
311 if (!g.isOne())
312 g = 1/g;
314 int exp= 0;
315 if (n & 1)
316 {
317 result= g;
318 exp= 1;
319 }
321
322 for (int i= 1; i <= l; i++)
323 {
324 h= mulNTL (g, mod (F, power (x, (1 << i))));
325 h= mod (h, power (x, (1 << i)) - 1);
326 h= div (h, power (x, (1 << (i - 1))));
327 g -= power (x, (1 << (i - 1)))*
328 mulFLINTQTrunc (g, h, 1 << (i-1));
329
330 if (n & (1 << i))
331 {
332 if (exp)
333 {
334 h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
335 h= mod (h, power (x, exp + (1 << i)) - 1);
336 h= div (h, power (x, exp));
337 result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
338 exp += (1 << i);
339 }
340 else
341 {
342 exp= (1 << i);
343 result= g;
344 }
345 }
346 }
347
348 return result;
349}
350
351void
354{
355 ASSERT (F.level() == G.level(), "F and G have different level");
356 CanonicalForm A= F;
358 Variable x= A.mvar();
359 int degA= degree (A);
360 int degB= degree (B);
361 int m= degA - degB;
362
363 if (m < 0)
364 {
365 R= A;
366 Q= 0;
367 return;
368 }
369
370 if (degB <= 1)
371 divrem (A, B, Q, R);
372 else
373 {
374 R= uniReverse (A, degA, x);
375
376 CanonicalForm revB= uniReverse (B, degB, x);
377 revB= newtonInverse (revB, m + 1, x);
378 Q= mulFLINTQTrunc (R, revB, m + 1);
379 Q= uniReverse (Q, m, x);
380
381 R= A - mulNTL (Q, B);
382 }
383}
384
385void
387{
388 ASSERT (F.level() == G.level(), "F and G have different level");
389 CanonicalForm A= F;
391 Variable x= A.mvar();
392 int degA= degree (A);
393 int degB= degree (B);
394 int m= degA - degB;
395
396 if (m < 0)
397 {
398 Q= 0;
399 return;
400 }
401
402 if (degB <= 1)
403 Q= div (A, B);
404 else
405 {
406 CanonicalForm R= uniReverse (A, degA, x);
407 CanonicalForm revB= uniReverse (B, degB, x);
408 revB= newtonInverse (revB, m + 1, x);
409 Q= mulFLINTQTrunc (R, revB, m + 1);
410 Q= uniReverse (Q, m, x);
411 }
412}
413
414#endif
415
417mulNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
418{
420 return F*G;
421 if (getCharacteristic() == 0)
422 {
424 if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
426 {
427 if (b.getp() != 0)
428 {
430 bool is_rat= isOn (SW_RATIONAL);
431 if (!is_rat)
432 On (SW_RATIONAL);
433 mipo *=bCommonDen (mipo);
434 if (!is_rat)
436#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
437 fmpz_t FLINTp;
438 fmpz_mod_poly_t FLINTmipo;
439 fq_ctx_t fq_con;
440 fq_poly_t FLINTF, FLINTG;
441
442 fmpz_init (FLINTp);
443
444 convertCF2initFmpz (FLINTp, b.getpk());
445
446 convertFacCF2Fmpz_mod_poly_t (FLINTmipo, mipo, FLINTp);
447
448 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
449 fmpz_mod_ctx_t fmpz_ctx;
450 fmpz_mod_ctx_init(fmpz_ctx,FLINTp);
451 fq_ctx_init_modulus (fq_con, FLINTmipo, fmpz_ctx, "Z");
452 #else
453 fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
454 #endif
455
456 convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
458
459 fq_poly_mul (FLINTF, FLINTF, FLINTG, fq_con);
460
462 alpha, fq_con);
463
464 fmpz_clear (FLINTp);
465 fq_poly_clear (FLINTF, fq_con);
466 fq_poly_clear (FLINTG, fq_con);
467 fq_ctx_clear (fq_con);
468 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
469 fmpz_mod_poly_clear (FLINTmipo,fmpz_ctx);
470 fmpz_mod_ctx_clear(fmpz_ctx);
471 #else
472 fmpz_mod_poly_clear(FLINTmipo);
473 #endif
474 return b (result);
475#endif
476#ifdef HAVE_NTL
477 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
478 ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (mipo));
479 ZZ_pE::init (NTLmipo);
480 ZZ_pEX NTLg= convertFacCF2NTLZZ_pEX (G, NTLmipo);
481 ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
482 mul (NTLf, NTLf, NTLg);
483
484 return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
485#endif
486 }
487#ifdef HAVE_FLINT
489 return result;
490#else
491 return F*G;
492#endif
493 }
494 else if (!F.inCoeffDomain() && !G.inCoeffDomain())
495 {
496#ifdef HAVE_FLINT
497 if (b.getp() != 0)
498 {
499 fmpz_t FLINTpk;
500 fmpz_init (FLINTpk);
501 convertCF2initFmpz (FLINTpk, b.getpk());
502 fmpz_mod_poly_t FLINTF, FLINTG;
503 convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
504 convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
505 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
506 fmpz_mod_ctx_t fmpz_ctx;
507 fmpz_mod_ctx_init(fmpz_ctx,FLINTpk);
508 fmpz_mod_poly_mul (FLINTF, FLINTF, FLINTG, fmpz_ctx);
509 #else
510 fmpz_mod_poly_mul (FLINTF, FLINTF, FLINTG);
511 #endif
513 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
514 fmpz_mod_poly_clear (FLINTG,fmpz_ctx);
515 fmpz_mod_poly_clear (FLINTF,fmpz_ctx);
516 fmpz_mod_ctx_clear(fmpz_ctx);
517 #else
518 fmpz_mod_poly_clear (FLINTG);
519 fmpz_mod_poly_clear (FLINTF);
520 #endif
521 fmpz_clear (FLINTpk);
522 return result;
523 }
524 return mulFLINTQ (F, G);
525#endif
526#ifdef HAVE_NTL
527 if (b.getp() != 0)
528 {
529 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
530 ZZX ZZf= convertFacCF2NTLZZX (F);
531 ZZX ZZg= convertFacCF2NTLZZX (G);
532 ZZ_pX NTLf= to_ZZ_pX (ZZf);
533 ZZ_pX NTLg= to_ZZ_pX (ZZg);
534 mul (NTLf, NTLf, NTLg);
535 return b (convertNTLZZX2CF (to_ZZX (NTLf), F.mvar()));
536 }
537 return F*G;
538#endif
539 }
540 if (b.getp() != 0)
541 {
542 if (!F.inBaseDomain() && !G.inBaseDomain())
543 {
545 {
546#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
547 fmpz_t FLINTp;
548 fmpz_mod_poly_t FLINTmipo;
549 fq_ctx_t fq_con;
550
551 fmpz_init (FLINTp);
552 convertCF2initFmpz (FLINTp, b.getpk());
553
555 bool rat=isOn(SW_RATIONAL);
558 mipo *= cd;
559 if (!rat) Off(SW_RATIONAL);
560 convertFacCF2Fmpz_mod_poly_t (FLINTmipo, mipo, FLINTp);
561
562 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
563 fmpz_mod_ctx_t fmpz_ctx;
564 fmpz_mod_ctx_init(fmpz_ctx,FLINTp);
565 fq_ctx_init_modulus (fq_con, FLINTmipo, fmpz_ctx, "Z");
566 #else
567 fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
568 #endif
569
571
572 if (F.inCoeffDomain() && !G.inCoeffDomain())
573 {
574 fq_poly_t FLINTG;
575 fmpz_poly_t FLINTF;
576 convertFacCF2Fmpz_poly_t (FLINTF, F);
578
579 fq_poly_scalar_mul_fq (FLINTG, FLINTG, FLINTF, fq_con);
580
581 result= convertFq_poly_t2FacCF (FLINTG, G.mvar(), alpha, fq_con);
582 fmpz_poly_clear (FLINTF);
583 fq_poly_clear (FLINTG, fq_con);
584 }
585 else if (!F.inCoeffDomain() && G.inCoeffDomain())
586 {
587 fq_poly_t FLINTF;
588 fmpz_poly_t FLINTG;
589
590 convertFacCF2Fmpz_poly_t (FLINTG, G);
591 convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
592
593 fq_poly_scalar_mul_fq (FLINTF, FLINTF, FLINTG, fq_con);
594
596 fmpz_poly_clear (FLINTG);
597 fq_poly_clear (FLINTF, fq_con);
598 }
599 else
600 {
601 fq_t FLINTF, FLINTG;
602
603 convertFacCF2Fq_t (FLINTF, F, fq_con);
604 convertFacCF2Fq_t (FLINTG, G, fq_con);
605
606 fq_mul (FLINTF, FLINTF, FLINTG, fq_con);
607
608 result= convertFq_t2FacCF (FLINTF, alpha);
609 fq_clear (FLINTF, fq_con);
610 fq_clear (FLINTG, fq_con);
611 }
612
613 fmpz_clear (FLINTp);
614 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
615 fmpz_mod_poly_clear (FLINTmipo,fmpz_ctx);
616 fmpz_mod_ctx_clear(fmpz_ctx);
617 #else
618 fmpz_mod_poly_clear (FLINTmipo);
619 #endif
620 fq_ctx_clear (fq_con);
621
622 return b (result);
623#endif
624#ifdef HAVE_NTL
625 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
626 ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
627 ZZ_pE::init (NTLmipo);
628
629 if (F.inCoeffDomain() && !G.inCoeffDomain())
630 {
631 ZZ_pEX NTLg= convertFacCF2NTLZZ_pEX (G, NTLmipo);
632 ZZ_pX NTLf= convertFacCF2NTLZZpX (F);
633 mul (NTLg, to_ZZ_pE (NTLf), NTLg);
634 return b (convertNTLZZ_pEX2CF (NTLg, G.mvar(), alpha));
635 }
636 else if (!F.inCoeffDomain() && G.inCoeffDomain())
637 {
638 ZZ_pX NTLg= convertFacCF2NTLZZpX (G);
639 ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
640 mul (NTLf, NTLf, to_ZZ_pE (NTLg));
641 return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
642 }
643 else
644 {
645 ZZ_pX NTLg= convertFacCF2NTLZZpX (G);
646 ZZ_pX NTLf= convertFacCF2NTLZZpX (F);
647 ZZ_pE result;
648 mul (result, to_ZZ_pE (NTLg), to_ZZ_pE (NTLf));
649 return b (convertNTLZZpX2CF (rep (result), alpha));
650 }
651#endif
652 }
653 }
654 return b (F*G);
655 }
656 return F*G;
657 }
658 else if (F.inCoeffDomain() || G.inCoeffDomain())
659 return F*G;
660 ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
661 ASSERT (F.level() == G.level(), "expected polys of same level");
662#ifdef HAVE_NTL
663#if (!defined(HAVE_FLINT) || __FLINT_RELEASE < 20400)
665 {
667 zz_p::init (getCharacteristic());
668 }
669#endif
670#endif
674 {
675 if (!getReduce (alpha))
676 {
677 result= 0;
678 for (CFIterator i= F; i.hasTerms(); i++)
679 result += i.coeff()*G*power (F.mvar(),i.exp());
680 return result;
681 }
682#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
683 nmod_poly_t FLINTmipo;
684 fq_nmod_ctx_t fq_con;
685
686 nmod_poly_init (FLINTmipo, getCharacteristic());
688
689 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
690
691 fq_nmod_poly_t FLINTF, FLINTG;
694
695 fq_nmod_poly_mul (FLINTF, FLINTF, FLINTG, fq_con);
696
698
699 fq_nmod_poly_clear (FLINTF, fq_con);
700 fq_nmod_poly_clear (FLINTG, fq_con);
701 nmod_poly_clear (FLINTmipo);
703 return result;
704#elif defined(AHVE_NTL)
705 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
706 zz_pE::init (NTLMipo);
707 zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
708 zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
709 mul (NTLF, NTLF, NTLG);
711 return result;
712#endif
713 }
714 else
715 {
716#ifdef HAVE_FLINT
717 nmod_poly_t FLINTF, FLINTG;
718 convertFacCF2nmod_poly_t (FLINTF, F);
719 convertFacCF2nmod_poly_t (FLINTG, G);
720 nmod_poly_mul (FLINTF, FLINTF, FLINTG);
721 result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
722 nmod_poly_clear (FLINTF);
723 nmod_poly_clear (FLINTG);
724 return result;
725#endif
726#ifdef HAVE_NTL
727 zz_pX NTLF= convertFacCF2NTLzzpX (F);
728 zz_pX NTLG= convertFacCF2NTLzzpX (G);
729 mul (NTLF, NTLF, NTLG);
730 return convertNTLzzpX2CF(NTLF, F.mvar());
731#endif
732 }
733 return F*G;
734}
735
737modNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
738{
740 return mod (F, G);
741 if (F.inCoeffDomain() && G.isUnivariate() && !G.inCoeffDomain())
742 {
743 if (b.getp() != 0)
744 return b(F);
745 return F;
746 }
747 else if (F.inCoeffDomain() && G.inCoeffDomain())
748 {
749 if (b.getp() != 0)
750 return b(F%G);
751 return mod (F, G);
752 }
753 else if (F.isUnivariate() && G.inCoeffDomain())
754 {
755 if (b.getp() != 0)
756 return b(F%G);
757 return mod (F,G);
758 }
759
760 if (getCharacteristic() == 0)
761 {
763 if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
764 {
765#ifdef HAVE_FLINT
766 if (b.getp() != 0)
767 {
768 fmpz_t FLINTpk;
769 fmpz_init (FLINTpk);
770 convertCF2initFmpz (FLINTpk, b.getpk());
771 fmpz_mod_poly_t FLINTF, FLINTG;
772 convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
773 convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
774 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
775 fmpz_mod_ctx_t fmpz_ctx;
776 fmpz_mod_ctx_init(fmpz_ctx,FLINTpk);
777 fmpz_mod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG, fmpz_ctx);
778 #else
779 fmpz_mod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
780 #endif
782 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
783 fmpz_mod_poly_clear (FLINTG, fmpz_ctx);
784 fmpz_mod_poly_clear (FLINTF, fmpz_ctx);
785 fmpz_mod_ctx_clear(fmpz_ctx);
786 #else
787 fmpz_mod_poly_clear (FLINTG);
788 fmpz_mod_poly_clear (FLINTF);
789 #endif
790 fmpz_clear (FLINTpk);
791 return result;
792 }
793 return modFLINTQ (F, G);
794#else
795 if (b.getp() != 0)
796 {
797 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
798 ZZX ZZf= convertFacCF2NTLZZX (F);
799 ZZX ZZg= convertFacCF2NTLZZX (G);
800 ZZ_pX NTLf= to_ZZ_pX (ZZf);
801 ZZ_pX NTLg= to_ZZ_pX (ZZg);
802 rem (NTLf, NTLf, NTLg);
803 return b (convertNTLZZX2CF (to_ZZX (NTLf), F.mvar()));
804 }
805 return mod (F, G);
806#endif
807 }
808 else
809 {
810 if (b.getp() != 0)
811 {
812#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
813 fmpz_t FLINTp;
814 fmpz_mod_poly_t FLINTmipo;
815 fq_ctx_t fq_con;
816 fq_poly_t FLINTF, FLINTG;
817
818 fmpz_init (FLINTp);
819
820 convertCF2initFmpz (FLINTp, b.getpk());
821
823 bool rat=isOn(SW_RATIONAL);
826 mipo *= cd;
827 if (!rat) Off(SW_RATIONAL);
828 convertFacCF2Fmpz_mod_poly_t (FLINTmipo, mipo, FLINTp);
829
830 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
831 fmpz_mod_ctx_t fmpz_ctx;
832 fmpz_mod_ctx_init(fmpz_ctx,FLINTp);
833 fq_ctx_init_modulus (fq_con, FLINTmipo, fmpz_ctx, "Z");
834 #else
835 fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
836 #endif
837
838 convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
840
841 fq_poly_rem (FLINTF, FLINTF, FLINTG, fq_con);
842
844 alpha, fq_con);
845
846 fmpz_clear (FLINTp);
847 fq_poly_clear (FLINTF, fq_con);
848 fq_poly_clear (FLINTG, fq_con);
849 fq_ctx_clear (fq_con);
850 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
851 fmpz_mod_poly_clear (FLINTmipo, fmpz_ctx);
852 fmpz_mod_ctx_clear(fmpz_ctx);
853 #else
854 fmpz_mod_poly_clear (FLINTmipo);
855 #endif
856
857 return b(result);
858#else
859 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
860 ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
861 ZZ_pE::init (NTLmipo);
862 ZZ_pEX NTLg= convertFacCF2NTLZZ_pEX (G, NTLmipo);
863 ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
864 rem (NTLf, NTLf, NTLg);
865 return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
866#endif
867 }
868#ifdef HAVE_FLINT
870 newtonDivrem (F, G, Q, R);
871 return R;
872#else
873 return mod (F,G);
874#endif
875 }
876 }
877
878 ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
879 ASSERT (F.level() == G.level(), "expected polys of same level");
880#if (!defined(HAVE_FLINT) || __FLINT_RELEASE < 20400)
882 {
884 zz_p::init (getCharacteristic());
885 }
886#endif
890 {
891#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
892 nmod_poly_t FLINTmipo;
893 fq_nmod_ctx_t fq_con;
894
895 nmod_poly_init (FLINTmipo, getCharacteristic());
897
898 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
899
900 fq_nmod_poly_t FLINTF, FLINTG;
903
904 fq_nmod_poly_rem (FLINTF, FLINTF, FLINTG, fq_con);
905
907
908 fq_nmod_poly_clear (FLINTF, fq_con);
909 fq_nmod_poly_clear (FLINTG, fq_con);
910 nmod_poly_clear (FLINTmipo);
912#else
913 zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
914 zz_pE::init (NTLMipo);
915 zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
916 zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
917 rem (NTLF, NTLF, NTLG);
919#endif
920 }
921 else
922 {
923#ifdef HAVE_FLINT
924 nmod_poly_t FLINTF, FLINTG;
925 convertFacCF2nmod_poly_t (FLINTF, F);
926 convertFacCF2nmod_poly_t (FLINTG, G);
927 nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
928 result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
929 nmod_poly_clear (FLINTF);
930 nmod_poly_clear (FLINTG);
931#else
932 zz_pX NTLF= convertFacCF2NTLzzpX (F);
933 zz_pX NTLG= convertFacCF2NTLzzpX (G);
934 rem (NTLF, NTLF, NTLG);
935 result= convertNTLzzpX2CF(NTLF, F.mvar());
936#endif
937 }
938 return result;
939}
940
942divNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
943{
945 return div (F, G);
946 if (F.inCoeffDomain() && G.isUnivariate() && !G.inCoeffDomain())
947 {
948 return 0;
949 }
950 else if (F.inCoeffDomain() && G.inCoeffDomain())
951 {
952 if (b.getp() != 0)
953 {
954 if (!F.inBaseDomain() || !G.inBaseDomain())
955 {
959#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
960 fmpz_t FLINTp;
961 fmpz_mod_poly_t FLINTmipo;
962 fq_ctx_t fq_con;
963 fq_t FLINTF, FLINTG;
964
965 fmpz_init (FLINTp);
966 convertCF2initFmpz (FLINTp, b.getpk());
967
968 convertFacCF2Fmpz_mod_poly_t (FLINTmipo, getMipo (alpha), FLINTp);
969
970 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
971 fmpz_mod_ctx_t fmpz_ctx;
972 fmpz_mod_ctx_init(fmpz_ctx,FLINTp);
973 fq_ctx_init_modulus (fq_con, FLINTmipo, fmpz_ctx, "Z");
974 #else
975 fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
976 #endif
977
978 convertFacCF2Fq_t (FLINTF, F, fq_con);
979 convertFacCF2Fq_t (FLINTG, G, fq_con);
980
981 fq_inv (FLINTG, FLINTG, fq_con);
982 fq_mul (FLINTF, FLINTF, FLINTG, fq_con);
983
985
986 fmpz_clear (FLINTp);
987 fq_clear (FLINTF, fq_con);
988 fq_clear (FLINTG, fq_con);
989 fq_ctx_clear (fq_con);
990 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
991 fmpz_mod_poly_clear (FLINTmipo, fmpz_ctx);
992 fmpz_mod_ctx_clear(fmpz_ctx);
993 #else
994 fmpz_mod_poly_clear (FLINTmipo);
995 #endif
996 return b (result);
997#else
998 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
999 ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
1000 ZZ_pE::init (NTLmipo);
1001 ZZ_pX NTLg= convertFacCF2NTLZZpX (G);
1002 ZZ_pX NTLf= convertFacCF2NTLZZpX (F);
1003 ZZ_pE result;
1004 div (result, to_ZZ_pE (NTLf), to_ZZ_pE (NTLg));
1005 return b (convertNTLZZpX2CF (rep (result), alpha));
1006#endif
1007 }
1008 return b(div (F,G));
1009 }
1010 return div (F, G);
1011 }
1012 else if (F.isUnivariate() && G.inCoeffDomain())
1013 {
1014 if (b.getp() != 0)
1015 {
1016 if (!G.inBaseDomain())
1017 {
1020#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
1021 fmpz_t FLINTp;
1022 fmpz_mod_poly_t FLINTmipo;
1023 fq_ctx_t fq_con;
1024 fq_poly_t FLINTF;
1025 fq_t FLINTG;
1026
1027 fmpz_init (FLINTp);
1028 convertCF2initFmpz (FLINTp, b.getpk());
1029
1030 convertFacCF2Fmpz_mod_poly_t (FLINTmipo, getMipo (alpha), FLINTp);
1031
1032 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
1033 fmpz_mod_ctx_t fmpz_ctx;
1034 fmpz_mod_ctx_init(fmpz_ctx,FLINTp);
1035 fq_ctx_init_modulus (fq_con, FLINTmipo, fmpz_ctx, "Z");
1036 #else
1037 fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
1038 #endif
1039
1040 convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
1041 convertFacCF2Fq_t (FLINTG, G, fq_con);
1042
1043 fq_inv (FLINTG, FLINTG, fq_con);
1044 fq_poly_scalar_mul_fq (FLINTF, FLINTF, FLINTG, fq_con);
1045
1047 alpha, fq_con);
1048
1049 fmpz_clear (FLINTp);
1050 fq_poly_clear (FLINTF, fq_con);
1051 fq_clear (FLINTG, fq_con);
1052 fq_ctx_clear (fq_con);
1053 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
1054 fmpz_mod_poly_clear (FLINTmipo, fmpz_ctx);
1055 fmpz_mod_ctx_clear(fmpz_ctx);
1056 #else
1057 fmpz_mod_poly_clear (FLINTmipo);
1058 #endif
1059 return b (result);
1060#else
1061 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
1062 ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
1063 ZZ_pE::init (NTLmipo);
1064 ZZ_pX NTLg= convertFacCF2NTLZZpX (G);
1065 ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
1066 div (NTLf, NTLf, to_ZZ_pE (NTLg));
1067 return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
1068#endif
1069 }
1070 return b(div (F,G));
1071 }
1072 return div (F, G);
1073 }
1074
1075 if (getCharacteristic() == 0)
1076 {
1077
1079 if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
1080 {
1081#ifdef HAVE_FLINT
1082 if (b.getp() != 0)
1083 {
1084 fmpz_t FLINTpk;
1085 fmpz_init (FLINTpk);
1086 convertCF2initFmpz (FLINTpk, b.getpk());
1087 fmpz_mod_poly_t FLINTF, FLINTG;
1088 convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
1089 convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
1090 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
1091 fmpz_mod_ctx_t fmpz_ctx;
1092 fmpz_mod_ctx_init(fmpz_ctx,FLINTpk);
1093 fmpz_mod_poly_divrem (FLINTF, FLINTG, FLINTF, FLINTG, fmpz_ctx);
1094 #else
1095 fmpz_mod_poly_divrem (FLINTF, FLINTG, FLINTF, FLINTG);
1096 #endif
1098 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
1099 fmpz_mod_poly_clear (FLINTG, fmpz_ctx);
1100 fmpz_mod_poly_clear (FLINTF, fmpz_ctx);
1101 fmpz_mod_ctx_clear(fmpz_ctx);
1102 #else
1103 fmpz_mod_poly_clear (FLINTG);
1104 fmpz_mod_poly_clear (FLINTF);
1105 #endif
1106 fmpz_clear (FLINTpk);
1107 return result;
1108 }
1109 return divFLINTQ (F,G);
1110#else
1111 if (b.getp() != 0)
1112 {
1113 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
1114 ZZX ZZf= convertFacCF2NTLZZX (F);
1115 ZZX ZZg= convertFacCF2NTLZZX (G);
1116 ZZ_pX NTLf= to_ZZ_pX (ZZf);
1117 ZZ_pX NTLg= to_ZZ_pX (ZZg);
1118 div (NTLf, NTLf, NTLg);
1119 return b (convertNTLZZX2CF (to_ZZX (NTLf), F.mvar()));
1120 }
1121 return div (F, G);
1122#endif
1123 }
1124 else
1125 {
1126 if (b.getp() != 0)
1127 {
1128#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
1129 fmpz_t FLINTp;
1130 fmpz_mod_poly_t FLINTmipo;
1131 fq_ctx_t fq_con;
1132 fq_poly_t FLINTF, FLINTG;
1133
1134 fmpz_init (FLINTp);
1135 convertCF2initFmpz (FLINTp, b.getpk());
1136
1137 convertFacCF2Fmpz_mod_poly_t (FLINTmipo, getMipo (alpha), FLINTp);
1138
1139 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
1140 fmpz_mod_ctx_t fmpz_ctx;
1141 fmpz_mod_ctx_init(fmpz_ctx,FLINTp);
1142 fq_ctx_init_modulus (fq_con, FLINTmipo, fmpz_ctx, "Z");
1143 #else
1144 fq_ctx_init_modulus (fq_con, FLINTmipo, "Z");
1145 #endif
1146
1147 convertFacCF2Fq_poly_t (FLINTF, F, fq_con);
1148 convertFacCF2Fq_poly_t (FLINTG, G, fq_con);
1149
1150 fq_poly_divrem (FLINTF, FLINTG, FLINTF, FLINTG, fq_con);
1151
1153 alpha, fq_con);
1154
1155 fmpz_clear (FLINTp);
1156 fq_poly_clear (FLINTF, fq_con);
1157 fq_poly_clear (FLINTG, fq_con);
1158 fq_ctx_clear (fq_con);
1159 #if (HAVE_FLINT && __FLINT_RELEASE >= 20700)
1160 fmpz_mod_poly_clear (FLINTmipo, fmpz_ctx);
1161 fmpz_mod_ctx_clear(fmpz_ctx);
1162 #else
1163 fmpz_mod_poly_clear (FLINTmipo);
1164 #endif
1165 return b (result);
1166#else
1167 ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
1168 ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (alpha)));
1169 ZZ_pE::init (NTLmipo);
1170 ZZ_pEX NTLg= convertFacCF2NTLZZ_pEX (G, NTLmipo);
1171 ZZ_pEX NTLf= convertFacCF2NTLZZ_pEX (F, NTLmipo);
1172 div (NTLf, NTLf, NTLg);
1173 return b (convertNTLZZ_pEX2CF (NTLf, F.mvar(), alpha));
1174#endif
1175 }
1176#ifdef HAVE_FLINT
1178 newtonDiv (F, G, Q);
1179 return Q;
1180#else
1181 return div (F,G);
1182#endif
1183 }
1184 }
1185
1186 ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
1187 ASSERT (F.level() == G.level(), "expected polys of same level");
1188#if (!defined(HAVE_FLINT) || __FLINT_RELEASE < 20400)
1190 {
1192 zz_p::init (getCharacteristic());
1193 }
1194#endif
1197 if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
1198 {
1199#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
1200 nmod_poly_t FLINTmipo;
1201 fq_nmod_ctx_t fq_con;
1202
1203 nmod_poly_init (FLINTmipo, getCharacteristic());
1204 convertFacCF2nmod_poly_t (FLINTmipo, getMipo (alpha));
1205
1206 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
1207
1208 fq_nmod_poly_t FLINTF, FLINTG;
1211
1212 fq_nmod_poly_divrem (FLINTF, FLINTG, FLINTF, FLINTG, fq_con);
1213
1215
1216 fq_nmod_poly_clear (FLINTF, fq_con);
1217 fq_nmod_poly_clear (FLINTG, fq_con);
1218 nmod_poly_clear (FLINTmipo);
1220#else
1221 zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
1222 zz_pE::init (NTLMipo);
1223 zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
1224 zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
1225 div (NTLF, NTLF, NTLG);
1226 result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
1227#endif
1228 }
1229 else
1230 {
1231#ifdef HAVE_FLINT
1232 nmod_poly_t FLINTF, FLINTG;
1233 convertFacCF2nmod_poly_t (FLINTF, F);
1234 convertFacCF2nmod_poly_t (FLINTG, G);
1235 nmod_poly_div (FLINTF, FLINTF, FLINTG);
1236 result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
1237 nmod_poly_clear (FLINTF);
1238 nmod_poly_clear (FLINTG);
1239#else
1240 zz_pX NTLF= convertFacCF2NTLzzpX (F);
1241 zz_pX NTLG= convertFacCF2NTLzzpX (G);
1242 div (NTLF, NTLF, NTLG);
1243 result= convertNTLzzpX2CF(NTLF, F.mvar());
1244#endif
1245 }
1246 return result;
1247}
1248
1249// end univariate polys
1250//*************************
1251// bivariate polys
1252
1253#ifdef HAVE_FLINT
1254void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
1255{
1256 int degAy= degree (A);
1257 nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
1258 result->length= d*(degAy + 1);
1259 flint_mpn_zero (result->coeffs, d*(degAy+1));
1260
1261 nmod_poly_t buf;
1262
1263 int k;
1264 for (CFIterator i= A; i.hasTerms(); i++)
1265 {
1266 convertFacCF2nmod_poly_t (buf, i.coeff());
1267 k= i.exp()*d;
1268 flint_mpn_copyi (result->coeffs+k, buf->coeffs, nmod_poly_length(buf));
1269
1271 }
1272 _nmod_poly_normalise (result);
1273}
1274
1275#if ( __FLINT_RELEASE >= 20400)
1276void
1277kronSubFq (fq_nmod_poly_t result, const CanonicalForm& A, int d,
1278 const fq_nmod_ctx_t fq_con)
1279{
1280 int degAy= degree (A);
1281 fq_nmod_poly_init2 (result, d*(degAy + 1), fq_con);
1282 _fq_nmod_poly_set_length (result, d*(degAy + 1), fq_con);
1283 _fq_nmod_vec_zero (result->coeffs, d*(degAy + 1), fq_con);
1284
1285 fq_nmod_poly_t buf1;
1286
1287 nmod_poly_t buf2;
1288
1289 int k;
1290
1291 for (CFIterator i= A; i.hasTerms(); i++)
1292 {
1293 if (i.coeff().inCoeffDomain())
1294 {
1295 convertFacCF2nmod_poly_t (buf2, i.coeff());
1296 fq_nmod_poly_init2 (buf1, 1, fq_con);
1297 fq_nmod_poly_set_coeff (buf1, 0, buf2, fq_con);
1299 }
1300 else
1302
1303 k= i.exp()*d;
1304 _fq_nmod_vec_set (result->coeffs+k, buf1->coeffs,
1305 fq_nmod_poly_length (buf1, fq_con), fq_con);
1306
1308 }
1309
1310 _fq_nmod_poly_normalise (result, fq_con);
1311}
1312#endif
1313
1314/*void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
1315{
1316 int degAy= degree (A);
1317 fmpq_poly_init2 (result, d1*(degAy + 1));
1318
1319 fmpq_poly_t buf;
1320 fmpq_t coeff;
1321
1322 int k, l, bufRepLength;
1323 CFIterator j;
1324 for (CFIterator i= A; i.hasTerms(); i++)
1325 {
1326 if (i.coeff().inCoeffDomain())
1327 {
1328 k= d1*i.exp();
1329 convertFacCF2Fmpq_poly_t (buf, i.coeff());
1330 bufRepLength= (int) fmpq_poly_length(buf);
1331 for (l= 0; l < bufRepLength; l++)
1332 {
1333 fmpq_poly_get_coeff_fmpq (coeff, buf, l);
1334 fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
1335 }
1336 fmpq_poly_clear (buf);
1337 }
1338 else
1339 {
1340 for (j= i.coeff(); j.hasTerms(); j++)
1341 {
1342 k= d1*i.exp();
1343 k += d2*j.exp();
1344 convertFacCF2Fmpq_poly_t (buf, j.coeff());
1345 bufRepLength= (int) fmpq_poly_length(buf);
1346 for (l= 0; l < bufRepLength; l++)
1347 {
1348 fmpq_poly_get_coeff_fmpq (coeff, buf, l);
1349 fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
1350 }
1351 fmpq_poly_clear (buf);
1352 }
1353 }
1354 }
1355 fmpq_clear (coeff);
1356 _fmpq_poly_normalise (result);
1357}*/
1358
1359void kronSubQa (fmpz_poly_t result, const CanonicalForm& A, int d1, int d2)
1360{
1361 int degAy= degree (A);
1362 fmpz_poly_init2 (result, d1*(degAy + 1));
1363 _fmpz_poly_set_length (result, d1*(degAy + 1));
1364
1365 fmpz_poly_t buf;
1366
1367 int k;
1368 CFIterator j;
1369 for (CFIterator i= A; i.hasTerms(); i++)
1370 {
1371 if (i.coeff().inCoeffDomain())
1372 {
1373 k= d1*i.exp();
1374 convertFacCF2Fmpz_poly_t (buf, i.coeff());
1375 _fmpz_vec_set (result->coeffs + k, buf->coeffs, buf->length);
1376 fmpz_poly_clear (buf);
1377 }
1378 else
1379 {
1380 for (j= i.coeff(); j.hasTerms(); j++)
1381 {
1382 k= d1*i.exp();
1383 k += d2*j.exp();
1384 convertFacCF2Fmpz_poly_t (buf, j.coeff());
1385 _fmpz_vec_set (result->coeffs + k, buf->coeffs, buf->length);
1386 fmpz_poly_clear (buf);
1387 }
1388 }
1389 }
1390 _fmpz_poly_normalise (result);
1391}
1392
1393void
1394kronSubReciproFp (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
1395 int d)
1396{
1397 int degAy= degree (A);
1398 mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1399 nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
1400 nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
1401
1402 nmod_poly_t buf;
1403
1404 int k, kk, j, bufRepLength;
1405 for (CFIterator i= A; i.hasTerms(); i++)
1406 {
1407 convertFacCF2nmod_poly_t (buf, i.coeff());
1408
1409 k= i.exp()*d;
1410 kk= (degAy - i.exp())*d;
1411 bufRepLength= (int) nmod_poly_length (buf);
1412 for (j= 0; j < bufRepLength; j++)
1413 {
1414 nmod_poly_set_coeff_ui (subA1, j + k,
1415 n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
1416 nmod_poly_get_coeff_ui (buf, j),
1418 )
1419 );
1420 nmod_poly_set_coeff_ui (subA2, j + kk,
1421 n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
1422 nmod_poly_get_coeff_ui (buf, j),
1424 )
1425 );
1426 }
1428 }
1429 _nmod_poly_normalise (subA1);
1430 _nmod_poly_normalise (subA2);
1431}
1432
1433#if ( __FLINT_RELEASE >= 20400)
1434void
1435kronSubReciproFq (fq_nmod_poly_t subA1, fq_nmod_poly_t subA2,
1436 const CanonicalForm& A, int d, const fq_nmod_ctx_t fq_con)
1437{
1438 int degAy= degree (A);
1439 fq_nmod_poly_init2 (subA1, d*(degAy + 2), fq_con);
1440 fq_nmod_poly_init2 (subA2, d*(degAy + 2), fq_con);
1441
1442 _fq_nmod_poly_set_length (subA1, d*(degAy + 2), fq_con);
1443 _fq_nmod_vec_zero (subA1->coeffs, d*(degAy + 2), fq_con);
1444
1445 _fq_nmod_poly_set_length (subA2, d*(degAy + 2), fq_con);
1446 _fq_nmod_vec_zero (subA2->coeffs, d*(degAy + 2), fq_con);
1447
1448 fq_nmod_poly_t buf1;
1449
1450 nmod_poly_t buf2;
1451
1452 int k, kk;
1453 for (CFIterator i= A; i.hasTerms(); i++)
1454 {
1455 if (i.coeff().inCoeffDomain())
1456 {
1457 convertFacCF2nmod_poly_t (buf2, i.coeff());
1458 fq_nmod_poly_init2 (buf1, 1, fq_con);
1459 fq_nmod_poly_set_coeff (buf1, 0, buf2, fq_con);
1461 }
1462 else
1464
1465 k= i.exp()*d;
1466 kk= (degAy - i.exp())*d;
1467 _fq_nmod_vec_add (subA1->coeffs+k, subA1->coeffs+k, buf1->coeffs,
1468 fq_nmod_poly_length(buf1, fq_con), fq_con);
1469 _fq_nmod_vec_add (subA2->coeffs+kk, subA2->coeffs+kk, buf1->coeffs,
1470 fq_nmod_poly_length(buf1, fq_con), fq_con);
1471
1473 }
1474 _fq_nmod_poly_normalise (subA1, fq_con);
1475 _fq_nmod_poly_normalise (subA2, fq_con);
1476}
1477#endif
1478
1479void
1480kronSubReciproQ (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
1481 int d)
1482{
1483 int degAy= degree (A);
1484 fmpz_poly_init2 (subA1, d*(degAy + 2));
1485 fmpz_poly_init2 (subA2, d*(degAy + 2));
1486
1487 fmpz_poly_t buf;
1488
1489 int k, kk;
1490 for (CFIterator i= A; i.hasTerms(); i++)
1491 {
1492 convertFacCF2Fmpz_poly_t (buf, i.coeff());
1493
1494 k= i.exp()*d;
1495 kk= (degAy - i.exp())*d;
1496 _fmpz_vec_add (subA1->coeffs+k, subA1->coeffs + k, buf->coeffs, buf->length);
1497 _fmpz_vec_add (subA2->coeffs+kk, subA2->coeffs + kk, buf->coeffs, buf->length);
1498 fmpz_poly_clear (buf);
1499 }
1500
1501 _fmpz_poly_normalise (subA1);
1502 _fmpz_poly_normalise (subA2);
1503}
1504
1505CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
1506{
1507 Variable y= Variable (2);
1508 Variable x= Variable (1);
1509
1510 fmpz_poly_t buf;
1512 int i= 0;
1513 int degf= fmpz_poly_degree(F);
1514 int k= 0;
1515 int degfSubK, repLength;
1516 while (degf >= k)
1517 {
1518 degfSubK= degf - k;
1519 if (degfSubK >= d)
1520 repLength= d;
1521 else
1522 repLength= degfSubK + 1;
1523
1524 fmpz_poly_init2 (buf, repLength);
1525 _fmpz_poly_set_length (buf, repLength);
1526 _fmpz_vec_set (buf->coeffs, F->coeffs+k, repLength);
1527 _fmpz_poly_normalise (buf);
1528
1530 i++;
1531 k= d*i;
1532 fmpz_poly_clear (buf);
1533 }
1534
1535 return result;
1536}
1537
1538/*CanonicalForm
1539reverseSubstQa (const fmpq_poly_t F, int d1, int d2, const Variable& alpha,
1540 const fmpq_poly_t mipo)
1541{
1542 Variable y= Variable (2);
1543 Variable x= Variable (1);
1544
1545 fmpq_poly_t f;
1546 fmpq_poly_init (f);
1547 fmpq_poly_set (f, F);
1548
1549 fmpq_poly_t buf;
1550 CanonicalForm result= 0, result2;
1551 int i= 0;
1552 int degf= fmpq_poly_degree(f);
1553 int k= 0;
1554 int degfSubK;
1555 int repLength;
1556 fmpq_t coeff;
1557 while (degf >= k)
1558 {
1559 degfSubK= degf - k;
1560 if (degfSubK >= d1)
1561 repLength= d1;
1562 else
1563 repLength= degfSubK + 1;
1564
1565 fmpq_init (coeff);
1566 int j= 0;
1567 int l;
1568 result2= 0;
1569 while (j*d2 < repLength)
1570 {
1571 fmpq_poly_init2 (buf, d2);
1572 for (l= 0; l < d2; l++)
1573 {
1574 fmpq_poly_get_coeff_fmpq (coeff, f, k + j*d2 + l);
1575 fmpq_poly_set_coeff_fmpq (buf, l, coeff);
1576 }
1577 _fmpq_poly_normalise (buf);
1578 fmpq_poly_rem (buf, buf, mipo);
1579 result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
1580 j++;
1581 fmpq_poly_clear (buf);
1582 }
1583 if (repLength - j*d2 != 0 && j*d2 - repLength < d2)
1584 {
1585 j--;
1586 repLength -= j*d2;
1587 fmpq_poly_init2 (buf, repLength);
1588 j++;
1589 for (l= 0; l < repLength; l++)
1590 {
1591 fmpq_poly_get_coeff_fmpq (coeff, f, k + j*d2 + l);
1592 fmpq_poly_set_coeff_fmpq (buf, l, coeff);
1593 }
1594 _fmpq_poly_normalise (buf);
1595 fmpq_poly_rem (buf, buf, mipo);
1596 result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
1597 fmpq_poly_clear (buf);
1598 }
1599 fmpq_clear (coeff);
1600
1601 result += result2*power (y, i);
1602 i++;
1603 k= d1*i;
1604 }
1605
1606 fmpq_poly_clear (f);
1607 return result;
1608}*/
1609
1611reverseSubstQa (const fmpz_poly_t F, int d1, int d2, const Variable& alpha,
1612 const fmpq_poly_t mipo)
1613{
1614 Variable y= Variable (2);
1615 Variable x= Variable (1);
1616
1617 fmpq_poly_t buf;
1618 CanonicalForm result= 0, result2;
1619 int i= 0;
1620 int degf= fmpz_poly_degree(F);
1621 int k= 0;
1622 int degfSubK;
1623 int repLength;
1624 while (degf >= k)
1625 {
1626 degfSubK= degf - k;
1627 if (degfSubK >= d1)
1628 repLength= d1;
1629 else
1630 repLength= degfSubK + 1;
1631
1632 int j= 0;
1633 result2= 0;
1634 while (j*d2 < repLength)
1635 {
1636 fmpq_poly_init2 (buf, d2);
1637 _fmpq_poly_set_length (buf, d2);
1638 _fmpz_vec_set (buf->coeffs, F->coeffs + k + j*d2, d2);
1639 _fmpq_poly_normalise (buf);
1640 fmpq_poly_rem (buf, buf, mipo);
1641 result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
1642 j++;
1643 fmpq_poly_clear (buf);
1644 }
1645 if (repLength - j*d2 != 0 && j*d2 - repLength < d2)
1646 {
1647 j--;
1648 repLength -= j*d2;
1649 fmpq_poly_init2 (buf, repLength);
1650 _fmpq_poly_set_length (buf, repLength);
1651 j++;
1652 _fmpz_vec_set (buf->coeffs, F->coeffs + k + j*d2, repLength);
1653 _fmpq_poly_normalise (buf);
1654 fmpq_poly_rem (buf, buf, mipo);
1655 result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
1656 fmpq_poly_clear (buf);
1657 }
1658
1659 result += result2*power (y, i);
1660 i++;
1661 k= d1*i;
1662 }
1663
1664 return result;
1665}
1666
1668reverseSubstReciproFp (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
1669{
1670 Variable y= Variable (2);
1671 Variable x= Variable (1);
1672
1673 nmod_poly_t f, g;
1674 mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1675 nmod_poly_init_preinv (f, getCharacteristic(), ninv);
1676 nmod_poly_init_preinv (g, getCharacteristic(), ninv);
1677 nmod_poly_set (f, F);
1678 nmod_poly_set (g, G);
1679 int degf= nmod_poly_degree(f);
1680 int degg= nmod_poly_degree(g);
1681
1682
1683 nmod_poly_t buf1,buf2, buf3;
1684
1685 if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
1686 nmod_poly_fit_length (f,(long)d*(k+1));
1687
1689 int i= 0;
1690 int lf= 0;
1691 int lg= d*k;
1692 int degfSubLf= degf;
1693 int deggSubLg= degg-lg;
1694 int repLengthBuf2, repLengthBuf1, ind, tmp;
1695 while (degf >= lf || lg >= 0)
1696 {
1697 if (degfSubLf >= d)
1698 repLengthBuf1= d;
1699 else if (degfSubLf < 0)
1700 repLengthBuf1= 0;
1701 else
1702 repLengthBuf1= degfSubLf + 1;
1703 nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
1704
1705 for (ind= 0; ind < repLengthBuf1; ind++)
1706 nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
1707 _nmod_poly_normalise (buf1);
1708
1709 repLengthBuf1= nmod_poly_length (buf1);
1710
1711 if (deggSubLg >= d - 1)
1712 repLengthBuf2= d - 1;
1713 else if (deggSubLg < 0)
1714 repLengthBuf2= 0;
1715 else
1716 repLengthBuf2= deggSubLg + 1;
1717
1718 nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
1719 for (ind= 0; ind < repLengthBuf2; ind++)
1720 nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
1721
1722 _nmod_poly_normalise (buf2);
1723 repLengthBuf2= nmod_poly_length (buf2);
1724
1725 nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
1726 for (ind= 0; ind < repLengthBuf1; ind++)
1727 nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
1728 for (ind= repLengthBuf1; ind < d; ind++)
1729 nmod_poly_set_coeff_ui (buf3, ind, 0);
1730 for (ind= 0; ind < repLengthBuf2; ind++)
1731 nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
1732 _nmod_poly_normalise (buf3);
1733
1734 result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
1735 i++;
1736
1737
1738 lf= i*d;
1739 degfSubLf= degf - lf;
1740
1741 lg= d*(k-i);
1742 deggSubLg= degg - lg;
1743
1744 if (lg >= 0 && deggSubLg > 0)
1745 {
1746 if (repLengthBuf2 > degfSubLf + 1)
1747 degfSubLf= repLengthBuf2 - 1;
1748 tmp= tmin (repLengthBuf1, deggSubLg + 1);
1749 for (ind= 0; ind < tmp; ind++)
1750 nmod_poly_set_coeff_ui (g, ind + lg,
1751 n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
1752 nmod_poly_get_coeff_ui (buf1, ind),
1754 )
1755 );
1756 }
1757 if (lg < 0)
1758 {
1761 nmod_poly_clear (buf3);
1762 break;
1763 }
1764 if (degfSubLf >= 0)
1765 {
1766 for (ind= 0; ind < repLengthBuf2; ind++)
1767 nmod_poly_set_coeff_ui (f, ind + lf,
1768 n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
1769 nmod_poly_get_coeff_ui (buf2, ind),
1771 )
1772 );
1773 }
1776 nmod_poly_clear (buf3);
1777 }
1778
1781
1782 return result;
1783}
1784
1785#if ( __FLINT_RELEASE >= 20400)
1787reverseSubstReciproFq (const fq_nmod_poly_t F, const fq_nmod_poly_t G, int d,
1788 int k, const Variable& alpha, const fq_nmod_ctx_t fq_con)
1789{
1790 Variable y= Variable (2);
1791 Variable x= Variable (1);
1792
1793 fq_nmod_poly_t f, g;
1794 int degf= fq_nmod_poly_degree(F, fq_con);
1795 int degg= fq_nmod_poly_degree(G, fq_con);
1796
1797 fq_nmod_poly_t buf1,buf2, buf3;
1798
1801 fq_nmod_poly_set (f, F, fq_con);
1802 fq_nmod_poly_set (g, G, fq_con);
1803 if (fq_nmod_poly_length (f, fq_con) < (long) d*(k + 1)) //zero padding
1804 fq_nmod_poly_fit_length (f, (long) d*(k + 1), fq_con);
1805
1807 int i= 0;
1808 int lf= 0;
1809 int lg= d*k;
1810 int degfSubLf= degf;
1811 int deggSubLg= degg-lg;
1812 int repLengthBuf2, repLengthBuf1, tmp;
1813 while (degf >= lf || lg >= 0)
1814 {
1815 if (degfSubLf >= d)
1816 repLengthBuf1= d;
1817 else if (degfSubLf < 0)
1818 repLengthBuf1= 0;
1819 else
1820 repLengthBuf1= degfSubLf + 1;
1821 fq_nmod_poly_init2 (buf1, repLengthBuf1, fq_con);
1822 _fq_nmod_poly_set_length (buf1, repLengthBuf1, fq_con);
1823
1824 _fq_nmod_vec_set (buf1->coeffs, f->coeffs + lf, repLengthBuf1, fq_con);
1825 _fq_nmod_poly_normalise (buf1, fq_con);
1826
1827 repLengthBuf1= fq_nmod_poly_length (buf1, fq_con);
1828
1829 if (deggSubLg >= d - 1)
1830 repLengthBuf2= d - 1;
1831 else if (deggSubLg < 0)
1832 repLengthBuf2= 0;
1833 else
1834 repLengthBuf2= deggSubLg + 1;
1835
1836 fq_nmod_poly_init2 (buf2, repLengthBuf2, fq_con);
1837 _fq_nmod_poly_set_length (buf2, repLengthBuf2, fq_con);
1838 _fq_nmod_vec_set (buf2->coeffs, g->coeffs + lg, repLengthBuf2, fq_con);
1839
1840 _fq_nmod_poly_normalise (buf2, fq_con);
1841 repLengthBuf2= fq_nmod_poly_length (buf2, fq_con);
1842
1843 fq_nmod_poly_init2 (buf3, repLengthBuf2 + d, fq_con);
1844 _fq_nmod_poly_set_length (buf3, repLengthBuf2 + d, fq_con);
1845 _fq_nmod_vec_set (buf3->coeffs, buf1->coeffs, repLengthBuf1, fq_con);
1846 _fq_nmod_vec_set (buf3->coeffs + d, buf2->coeffs, repLengthBuf2, fq_con);
1847
1848 _fq_nmod_poly_normalise (buf3, fq_con);
1849
1851 i++;
1852
1853
1854 lf= i*d;
1855 degfSubLf= degf - lf;
1856
1857 lg= d*(k - i);
1858 deggSubLg= degg - lg;
1859
1860 if (lg >= 0 && deggSubLg > 0)
1861 {
1862 if (repLengthBuf2 > degfSubLf + 1)
1863 degfSubLf= repLengthBuf2 - 1;
1864 tmp= tmin (repLengthBuf1, deggSubLg + 1);
1865 _fq_nmod_vec_sub (g->coeffs + lg, g->coeffs + lg, buf1-> coeffs,
1866 tmp, fq_con);
1867 }
1868 if (lg < 0)
1869 {
1872 fq_nmod_poly_clear (buf3, fq_con);
1873 break;
1874 }
1875 if (degfSubLf >= 0)
1876 _fq_nmod_vec_sub (f->coeffs + lf, f->coeffs + lf, buf2->coeffs,
1877 repLengthBuf2, fq_con);
1880 fq_nmod_poly_clear (buf3, fq_con);
1881 }
1882
1885
1886 return result;
1887}
1888#endif
1889
1891reverseSubstReciproQ (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
1892{
1893 Variable y= Variable (2);
1894 Variable x= Variable (1);
1895
1896 fmpz_poly_t f, g;
1897 fmpz_poly_init (f);
1898 fmpz_poly_init (g);
1899 fmpz_poly_set (f, F);
1900 fmpz_poly_set (g, G);
1901 int degf= fmpz_poly_degree(f);
1902 int degg= fmpz_poly_degree(g);
1903
1904 fmpz_poly_t buf1,buf2, buf3;
1905
1906 if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
1907 fmpz_poly_fit_length (f,(long)d*(k+1));
1908
1910 int i= 0;
1911 int lf= 0;
1912 int lg= d*k;
1913 int degfSubLf= degf;
1914 int deggSubLg= degg-lg;
1915 int repLengthBuf2, repLengthBuf1, ind, tmp;
1916 fmpz_t tmp1, tmp2;
1917 while (degf >= lf || lg >= 0)
1918 {
1919 if (degfSubLf >= d)
1920 repLengthBuf1= d;
1921 else if (degfSubLf < 0)
1922 repLengthBuf1= 0;
1923 else
1924 repLengthBuf1= degfSubLf + 1;
1925
1926 fmpz_poly_init2 (buf1, repLengthBuf1);
1927
1928 for (ind= 0; ind < repLengthBuf1; ind++)
1929 {
1930 fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1931 fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
1932 }
1933 _fmpz_poly_normalise (buf1);
1934
1935 repLengthBuf1= fmpz_poly_length (buf1);
1936
1937 if (deggSubLg >= d - 1)
1938 repLengthBuf2= d - 1;
1939 else if (deggSubLg < 0)
1940 repLengthBuf2= 0;
1941 else
1942 repLengthBuf2= deggSubLg + 1;
1943
1944 fmpz_poly_init2 (buf2, repLengthBuf2);
1945
1946 for (ind= 0; ind < repLengthBuf2; ind++)
1947 {
1948 fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1949 fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
1950 }
1951
1952 _fmpz_poly_normalise (buf2);
1953 repLengthBuf2= fmpz_poly_length (buf2);
1954
1955 fmpz_poly_init2 (buf3, repLengthBuf2 + d);
1956 for (ind= 0; ind < repLengthBuf1; ind++)
1957 {
1958 fmpz_poly_get_coeff_fmpz (tmp1, buf1, ind);
1959 fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
1960 }
1961 for (ind= repLengthBuf1; ind < d; ind++)
1962 fmpz_poly_set_coeff_ui (buf3, ind, 0);
1963 for (ind= 0; ind < repLengthBuf2; ind++)
1964 {
1965 fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
1966 fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
1967 }
1968 _fmpz_poly_normalise (buf3);
1969
1970 result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
1971 i++;
1972
1973
1974 lf= i*d;
1975 degfSubLf= degf - lf;
1976
1977 lg= d*(k-i);
1978 deggSubLg= degg - lg;
1979
1980 if (lg >= 0 && deggSubLg > 0)
1981 {
1982 if (repLengthBuf2 > degfSubLf + 1)
1983 degfSubLf= repLengthBuf2 - 1;
1984 tmp= tmin (repLengthBuf1, deggSubLg + 1);
1985 for (ind= 0; ind < tmp; ind++)
1986 {
1987 fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1988 fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
1989 fmpz_sub (tmp1, tmp1, tmp2);
1990 fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
1991 }
1992 }
1993 if (lg < 0)
1994 {
1995 fmpz_poly_clear (buf1);
1996 fmpz_poly_clear (buf2);
1997 fmpz_poly_clear (buf3);
1998 break;
1999 }
2000 if (degfSubLf >= 0)
2001 {
2002 for (ind= 0; ind < repLengthBuf2; ind++)
2003 {
2004 fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
2005 fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
2006 fmpz_sub (tmp1, tmp1, tmp2);
2007 fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
2008 }
2009 }
2010 fmpz_poly_clear (buf1);
2011 fmpz_poly_clear (buf2);
2012 fmpz_poly_clear (buf3);
2013 }
2014
2015 fmpz_poly_clear (f);
2016 fmpz_poly_clear (g);
2017 fmpz_clear (tmp1);
2018 fmpz_clear (tmp2);
2019
2020 return result;
2021}
2022
2023#if ( __FLINT_RELEASE >= 20400)
2025reverseSubstFq (const fq_nmod_poly_t F, int d, const Variable& alpha,
2026 const fq_nmod_ctx_t fq_con)
2027{
2028 Variable y= Variable (2);
2029 Variable x= Variable (1);
2030
2031 fq_nmod_poly_t buf;
2033 int i= 0;
2034 int degf= fq_nmod_poly_degree(F, fq_con);
2035 int k= 0;
2036 int degfSubK, repLength;
2037 while (degf >= k)
2038 {
2039 degfSubK= degf - k;
2040 if (degfSubK >= d)
2041 repLength= d;
2042 else
2043 repLength= degfSubK + 1;
2044
2045 fq_nmod_poly_init2 (buf, repLength, fq_con);
2046 _fq_nmod_poly_set_length (buf, repLength, fq_con);
2047 _fq_nmod_vec_set (buf->coeffs, F->coeffs+k, repLength, fq_con);
2048 _fq_nmod_poly_normalise (buf, fq_con);
2049
2051 i++;
2052 k= d*i;
2054 }
2055
2056 return result;
2057}
2058#endif
2059
2060CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
2061{
2062 Variable y= Variable (2);
2063 Variable x= Variable (1);
2064
2065 mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
2066
2067 nmod_poly_t buf;
2069 int i= 0;
2070 int degf= nmod_poly_degree(F);
2071 int k= 0;
2072 int degfSubK, repLength, j;
2073 while (degf >= k)
2074 {
2075 degfSubK= degf - k;
2076 if (degfSubK >= d)
2077 repLength= d;
2078 else
2079 repLength= degfSubK + 1;
2080
2081 nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
2082 for (j= 0; j < repLength; j++)
2083 nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (F, j + k));
2084 _nmod_poly_normalise (buf);
2085
2087 i++;
2088 k= d*i;
2090 }
2091
2092 return result;
2093}
2094
2098{
2099 int d1= degree (F, 1) + degree (G, 1) + 1;
2100 d1 /= 2;
2101 d1 += 1;
2102
2103 nmod_poly_t F1, F2;
2104 kronSubReciproFp (F1, F2, F, d1);
2105
2106 nmod_poly_t G1, G2;
2107 kronSubReciproFp (G1, G2, G, d1);
2108
2109 int k= d1*degree (M);
2110 nmod_poly_mullow (F1, F1, G1, (long) k);
2111
2112 int degtailF= degree (tailcoeff (F), 1);;
2113 int degtailG= degree (tailcoeff (G), 1);
2114 int taildegF= taildegree (F);
2115 int taildegG= taildegree (G);
2116
2117 int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
2118 + d1*(2+taildegF + taildegG);
2119 nmod_poly_mulhigh (F2, F2, G2, b);
2120 nmod_poly_shift_right (F2, F2, b);
2121 int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
2122
2123
2125
2128 nmod_poly_clear (G1);
2129 nmod_poly_clear (G2);
2130 return result;
2131}
2132
2136{
2137 CanonicalForm A= F;
2138 CanonicalForm B= G;
2139
2140 int degAx= degree (A, 1);
2141 int degAy= degree (A, 2);
2142 int degBx= degree (B, 1);
2143 int degBy= degree (B, 2);
2144 int d1= degAx + 1 + degBx;
2145 int d2= tmax (degAy, degBy);
2146
2147 if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
2148 return mulMod2FLINTFpReci (A, B, M);
2149
2150 nmod_poly_t FLINTA, FLINTB;
2151 kronSubFp (FLINTA, A, d1);
2152 kronSubFp (FLINTB, B, d1);
2153
2154 int k= d1*degree (M);
2155 nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
2156
2157 A= reverseSubstFp (FLINTA, d1);
2158
2159 nmod_poly_clear (FLINTA);
2160 nmod_poly_clear (FLINTB);
2161 return A;
2162}
2163
2164#if ( __FLINT_RELEASE >= 20400)
2167 CanonicalForm& M, const Variable& alpha,
2168 const fq_nmod_ctx_t fq_con)
2169{
2170 int d1= degree (F, 1) + degree (G, 1) + 1;
2171 d1 /= 2;
2172 d1 += 1;
2173
2174 fq_nmod_poly_t F1, F2;
2175 kronSubReciproFq (F1, F2, F, d1, fq_con);
2176
2177 fq_nmod_poly_t G1, G2;
2178 kronSubReciproFq (G1, G2, G, d1, fq_con);
2179
2180 int k= d1*degree (M);
2181 fq_nmod_poly_mullow (F1, F1, G1, (long) k, fq_con);
2182
2183 int degtailF= degree (tailcoeff (F), 1);
2184 int degtailG= degree (tailcoeff (G), 1);
2185 int taildegF= taildegree (F);
2186 int taildegG= taildegree (G);
2187
2188 int b= k + degtailF + degtailG - d1*(2+taildegF + taildegG);
2189
2190 fq_nmod_poly_reverse (F2, F2, fq_nmod_poly_length (F2, fq_con), fq_con);
2191 fq_nmod_poly_reverse (G2, G2, fq_nmod_poly_length (G2, fq_con), fq_con);
2192 fq_nmod_poly_mullow (F2, F2, G2, b+1, fq_con);
2193 fq_nmod_poly_reverse (F2, F2, b+1, fq_con);
2194
2195 int d2= tmax (fq_nmod_poly_degree (F2, fq_con)/d1,
2196 fq_nmod_poly_degree (F1, fq_con)/d1);
2197
2199
2204 return result;
2205}
2206
2209 CanonicalForm& M, const Variable& alpha,
2210 const fq_nmod_ctx_t fq_con)
2211{
2212 CanonicalForm A= F;
2213 CanonicalForm B= G;
2214
2215 int degAx= degree (A, 1);
2216 int degAy= degree (A, 2);
2217 int degBx= degree (B, 1);
2218 int degBy= degree (B, 2);
2219 int d1= degAx + 1 + degBx;
2220 int d2= tmax (degAy, degBy);
2221
2222 if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
2223 return mulMod2FLINTFqReci (A, B, M, alpha, fq_con);
2224
2225 fq_nmod_poly_t FLINTA, FLINTB;
2226 kronSubFq (FLINTA, A, d1, fq_con);
2227 kronSubFq (FLINTB, B, d1, fq_con);
2228
2229 int k= d1*degree (M);
2230 fq_nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k, fq_con);
2231
2232 A= reverseSubstFq (FLINTA, d1, alpha, fq_con);
2233
2234 fq_nmod_poly_clear (FLINTA, fq_con);
2235 fq_nmod_poly_clear (FLINTB, fq_con);
2236 return A;
2237}
2238#endif
2239
2243{
2244 int d1= degree (F, 1) + degree (G, 1) + 1;
2245 d1 /= 2;
2246 d1 += 1;
2247
2248 fmpz_poly_t F1, F2;
2249 kronSubReciproQ (F1, F2, F, d1);
2250
2251 fmpz_poly_t G1, G2;
2252 kronSubReciproQ (G1, G2, G, d1);
2253
2254 int k= d1*degree (M);
2255 fmpz_poly_mullow (F1, F1, G1, (long) k);
2256
2257 int degtailF= degree (tailcoeff (F), 1);;
2258 int degtailG= degree (tailcoeff (G), 1);
2259 int taildegF= taildegree (F);
2260 int taildegG= taildegree (G);
2261
2262 int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
2263 + d1*(2+taildegF + taildegG);
2264 fmpz_poly_mulhigh_n (F2, F2, G2, b);
2265 fmpz_poly_shift_right (F2, F2, b);
2266 int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
2267
2269
2270 fmpz_poly_clear (F1);
2271 fmpz_poly_clear (F2);
2272 fmpz_poly_clear (G1);
2273 fmpz_poly_clear (G2);
2274 return result;
2275}
2276
2280{
2281 CanonicalForm A= F;
2282 CanonicalForm B= G;
2283
2284 int degAx= degree (A, 1);
2285 int degBx= degree (B, 1);
2286 int d1= degAx + 1 + degBx;
2287
2290 A *= f;
2291 B *= g;
2292
2293 fmpz_poly_t FLINTA, FLINTB;
2294 kronSubQa (FLINTA, A, d1);
2295 kronSubQa (FLINTB, B, d1);
2296 int k= d1*degree (M);
2297
2298 fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
2299 A= reverseSubstQ (FLINTA, d1);
2300 fmpz_poly_clear (FLINTA);
2301 fmpz_poly_clear (FLINTB);
2302 return A/(f*g);
2303}
2304
2305/*CanonicalForm
2306mulMod2FLINTQa (const CanonicalForm& F, const CanonicalForm& G,
2307 const CanonicalForm& M)
2308{
2309 Variable a;
2310 if (!hasFirstAlgVar (F,a) && !hasFirstAlgVar (G, a))
2311 return mulMod2FLINTQ (F, G, M);
2312 CanonicalForm A= F;
2313
2314 int degFx= degree (F, 1);
2315 int degFa= degree (F, a);
2316 int degGx= degree (G, 1);
2317 int degGa= degree (G, a);
2318
2319 int d2= degFa+degGa+1;
2320 int d1= degFx + 1 + degGx;
2321 d1 *= d2;
2322
2323 fmpq_poly_t FLINTF, FLINTG;
2324 kronSubQa (FLINTF, F, d1, d2);
2325 kronSubQa (FLINTG, G, d1, d2);
2326
2327 fmpq_poly_mullow (FLINTF, FLINTF, FLINTG, d1*degree (M));
2328
2329 fmpq_poly_t mipo;
2330 convertFacCF2Fmpq_poly_t (mipo, getMipo (a));
2331 CanonicalForm result= reverseSubstQa (FLINTF, d1, d2, a, mipo);
2332 fmpq_poly_clear (FLINTF);
2333 fmpq_poly_clear (FLINTG);
2334 return result;
2335}*/
2336
2339 const CanonicalForm& M)
2340{
2341 Variable a;
2342 if (!hasFirstAlgVar (F,a) && !hasFirstAlgVar (G, a))
2343 return mulMod2FLINTQ (F, G, M);
2344 CanonicalForm A= F, B= G;
2345
2346 int degFx= degree (F, 1);
2347 int degFa= degree (F, a);
2348 int degGx= degree (G, 1);
2349 int degGa= degree (G, a);
2350
2351 int d2= degFa+degGa+1;
2352 int d1= degFx + 1 + degGx;
2353 d1 *= d2;
2354
2357 A *= f;
2358 B *= g;
2359
2360 fmpz_poly_t FLINTF, FLINTG;
2361 kronSubQa (FLINTF, A, d1, d2);
2362 kronSubQa (FLINTG, B, d1, d2);
2363
2364 fmpz_poly_mullow (FLINTF, FLINTF, FLINTG, d1*degree (M));
2365
2366 fmpq_poly_t mipo;
2368 A= reverseSubstQa (FLINTF, d1, d2, a, mipo);
2369 fmpz_poly_clear (FLINTF);
2370 fmpz_poly_clear (FLINTG);
2371 return A/(f*g);
2372}
2373
2374#endif
2375
2376#ifndef HAVE_FLINT
2377zz_pX kronSubFp (const CanonicalForm& A, int d)
2378{
2379 int degAy= degree (A);
2380 zz_pX result;
2381 result.rep.SetLength (d*(degAy + 1));
2382
2383 zz_p *resultp;
2384 resultp= result.rep.elts();
2385 zz_pX buf;
2386 zz_p *bufp;
2387 int j, k, bufRepLength;
2388
2389 for (CFIterator i= A; i.hasTerms(); i++)
2390 {
2391 if (i.coeff().inCoeffDomain())
2392 buf= convertFacCF2NTLzzpX (i.coeff());
2393 else
2394 buf= convertFacCF2NTLzzpX (i.coeff());
2395
2396 k= i.exp()*d;
2397 bufp= buf.rep.elts();
2398 bufRepLength= (int) buf.rep.length();
2399 for (j= 0; j < bufRepLength; j++)
2400 resultp [j + k]= bufp [j];
2401 }
2402 result.normalize();
2403
2404 return result;
2405}
2406#endif
2407
2408#if (!(HAVE_FLINT && __FLINT_RELEASE >= 20400))
2409zz_pEX kronSubFq (const CanonicalForm& A, int d, const Variable& alpha)
2410{
2411 int degAy= degree (A);
2412 zz_pEX result;
2413 result.rep.SetLength (d*(degAy + 1));
2414
2415 Variable v;
2416 zz_pE *resultp;
2417 resultp= result.rep.elts();
2418 zz_pEX buf1;
2419 zz_pE *buf1p;
2420 zz_pX buf2;
2421 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2422 int j, k, buf1RepLength;
2423
2424 for (CFIterator i= A; i.hasTerms(); i++)
2425 {
2426 if (i.coeff().inCoeffDomain())
2427 {
2428 buf2= convertFacCF2NTLzzpX (i.coeff());
2429 buf1= to_zz_pEX (to_zz_pE (buf2));
2430 }
2431 else
2432 buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
2433
2434 k= i.exp()*d;
2435 buf1p= buf1.rep.elts();
2436 buf1RepLength= (int) buf1.rep.length();
2437 for (j= 0; j < buf1RepLength; j++)
2438 resultp [j + k]= buf1p [j];
2439 }
2440 result.normalize();
2441
2442 return result;
2443}
2444
2445void
2446kronSubReciproFq (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
2447 const Variable& alpha)
2448{
2449 int degAy= degree (A);
2450 subA1.rep.SetLength ((long) d*(degAy + 2));
2451 subA2.rep.SetLength ((long) d*(degAy + 2));
2452
2453 Variable v;
2454 zz_pE *subA1p;
2455 zz_pE *subA2p;
2456 subA1p= subA1.rep.elts();
2457 subA2p= subA2.rep.elts();
2458 zz_pEX buf;
2459 zz_pE *bufp;
2460 zz_pX buf2;
2461 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2462 int j, k, kk, bufRepLength;
2463
2464 for (CFIterator i= A; i.hasTerms(); i++)
2465 {
2466 if (i.coeff().inCoeffDomain())
2467 {
2468 buf2= convertFacCF2NTLzzpX (i.coeff());
2469 buf= to_zz_pEX (to_zz_pE (buf2));
2470 }
2471 else
2472 buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
2473
2474 k= i.exp()*d;
2475 kk= (degAy - i.exp())*d;
2476 bufp= buf.rep.elts();
2477 bufRepLength= (int) buf.rep.length();
2478 for (j= 0; j < bufRepLength; j++)
2479 {
2480 subA1p [j + k] += bufp [j];
2481 subA2p [j + kk] += bufp [j];
2482 }
2483 }
2484 subA1.normalize();
2485 subA2.normalize();
2486}
2487#endif
2488
2489#ifndef HAVE_FLINT
2490void
2491kronSubReciproFp (zz_pX& subA1, zz_pX& subA2, const CanonicalForm& A, int d)
2492{
2493 int degAy= degree (A);
2494 subA1.rep.SetLength ((long) d*(degAy + 2));
2495 subA2.rep.SetLength ((long) d*(degAy + 2));
2496
2497 zz_p *subA1p;
2498 zz_p *subA2p;
2499 subA1p= subA1.rep.elts();
2500 subA2p= subA2.rep.elts();
2501 zz_pX buf;
2502 zz_p *bufp;
2503 int j, k, kk, bufRepLength;
2504
2505 for (CFIterator i= A; i.hasTerms(); i++)
2506 {
2507 buf= convertFacCF2NTLzzpX (i.coeff());
2508
2509 k= i.exp()*d;
2510 kk= (degAy - i.exp())*d;
2511 bufp= buf.rep.elts();
2512 bufRepLength= (int) buf.rep.length();
2513 for (j= 0; j < bufRepLength; j++)
2514 {
2515 subA1p [j + k] += bufp [j];
2516 subA2p [j + kk] += bufp [j];
2517 }
2518 }
2519 subA1.normalize();
2520 subA2.normalize();
2521}
2522#endif
2523
2524#if (!(HAVE_FLINT && __FLINT_RELEASE >= 20400))
2526reverseSubstReciproFq (const zz_pEX& F, const zz_pEX& G, int d, int k,
2527 const Variable& alpha)
2528{
2529 Variable y= Variable (2);
2530 Variable x= Variable (1);
2531
2532 zz_pEX f= F;
2533 zz_pEX g= G;
2534 int degf= deg(f);
2535 int degg= deg(g);
2536
2537 zz_pEX buf1;
2538 zz_pEX buf2;
2539 zz_pEX buf3;
2540
2541 zz_pE *buf1p;
2542 zz_pE *buf2p;
2543 zz_pE *buf3p;
2544 if (f.rep.length() < (long) d*(k+1)) //zero padding
2545 f.rep.SetLength ((long)d*(k+1));
2546
2547 zz_pE *gp= g.rep.elts();
2548 zz_pE *fp= f.rep.elts();
2550 int i= 0;
2551 int lf= 0;
2552 int lg= d*k;
2553 int degfSubLf= degf;
2554 int deggSubLg= degg-lg;
2555 int repLengthBuf2, repLengthBuf1, ind, tmp;
2556 zz_pE zzpEZero= zz_pE();
2557
2558 while (degf >= lf || lg >= 0)
2559 {
2560 if (degfSubLf >= d)
2561 repLengthBuf1= d;
2562 else if (degfSubLf < 0)
2563 repLengthBuf1= 0;
2564 else
2565 repLengthBuf1= degfSubLf + 1;
2566 buf1.rep.SetLength((long) repLengthBuf1);
2567
2568 buf1p= buf1.rep.elts();
2569 for (ind= 0; ind < repLengthBuf1; ind++)
2570 buf1p [ind]= fp [ind + lf];
2571 buf1.normalize();
2572
2573 repLengthBuf1= buf1.rep.length();
2574
2575 if (deggSubLg >= d - 1)
2576 repLengthBuf2= d - 1;
2577 else if (deggSubLg < 0)
2578 repLengthBuf2= 0;
2579 else
2580 repLengthBuf2= deggSubLg + 1;
2581
2582 buf2.rep.SetLength ((long) repLengthBuf2);
2583 buf2p= buf2.rep.elts();
2584 for (ind= 0; ind < repLengthBuf2; ind++)
2585 buf2p [ind]= gp [ind + lg];
2586 buf2.normalize();
2587
2588 repLengthBuf2= buf2.rep.length();
2589
2590 buf3.rep.SetLength((long) repLengthBuf2 + d);
2591 buf3p= buf3.rep.elts();
2592 buf2p= buf2.rep.elts();
2593 buf1p= buf1.rep.elts();
2594 for (ind= 0; ind < repLengthBuf1; ind++)
2595 buf3p [ind]= buf1p [ind];
2596 for (ind= repLengthBuf1; ind < d; ind++)
2597 buf3p [ind]= zzpEZero;
2598 for (ind= 0; ind < repLengthBuf2; ind++)
2599 buf3p [ind + d]= buf2p [ind];
2600 buf3.normalize();
2601
2602 result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
2603 i++;
2604
2605
2606 lf= i*d;
2607 degfSubLf= degf - lf;
2608
2609 lg= d*(k-i);
2610 deggSubLg= degg - lg;
2611
2612 buf1p= buf1.rep.elts();
2613
2614 if (lg >= 0 && deggSubLg > 0)
2615 {
2616 if (repLengthBuf2 > degfSubLf + 1)
2617 degfSubLf= repLengthBuf2 - 1;
2618 tmp= tmin (repLengthBuf1, deggSubLg + 1);
2619 for (ind= 0; ind < tmp; ind++)
2620 gp [ind + lg] -= buf1p [ind];
2621 }
2622
2623 if (lg < 0)
2624 break;
2625
2626 buf2p= buf2.rep.elts();
2627 if (degfSubLf >= 0)
2628 {
2629 for (ind= 0; ind < repLengthBuf2; ind++)
2630 fp [ind + lf] -= buf2p [ind];
2631 }
2632 }
2633
2634 return result;
2635}
2636#endif
2637
2638#ifndef HAVE_FLINT
2640reverseSubstReciproFp (const zz_pX& F, const zz_pX& G, int d, int k)
2641{
2642 Variable y= Variable (2);
2643 Variable x= Variable (1);
2644
2645 zz_pX f= F;
2646 zz_pX g= G;
2647 int degf= deg(f);
2648 int degg= deg(g);
2649
2650 zz_pX buf1;
2651 zz_pX buf2;
2652 zz_pX buf3;
2653
2654 zz_p *buf1p;
2655 zz_p *buf2p;
2656 zz_p *buf3p;
2657
2658 if (f.rep.length() < (long) d*(k+1)) //zero padding
2659 f.rep.SetLength ((long)d*(k+1));
2660
2661 zz_p *gp= g.rep.elts();
2662 zz_p *fp= f.rep.elts();
2664 int i= 0;
2665 int lf= 0;
2666 int lg= d*k;
2667 int degfSubLf= degf;
2668 int deggSubLg= degg-lg;
2669 int repLengthBuf2, repLengthBuf1, ind, tmp;
2670 zz_p zzpZero= zz_p();
2671 while (degf >= lf || lg >= 0)
2672 {
2673 if (degfSubLf >= d)
2674 repLengthBuf1= d;
2675 else if (degfSubLf < 0)
2676 repLengthBuf1= 0;
2677 else
2678 repLengthBuf1= degfSubLf + 1;
2679 buf1.rep.SetLength((long) repLengthBuf1);
2680
2681 buf1p= buf1.rep.elts();
2682 for (ind= 0; ind < repLengthBuf1; ind++)
2683 buf1p [ind]= fp [ind + lf];
2684 buf1.normalize();
2685
2686 repLengthBuf1= buf1.rep.length();
2687
2688 if (deggSubLg >= d - 1)
2689 repLengthBuf2= d - 1;
2690 else if (deggSubLg < 0)
2691 repLengthBuf2= 0;
2692 else
2693 repLengthBuf2= deggSubLg + 1;
2694
2695 buf2.rep.SetLength ((long) repLengthBuf2);
2696 buf2p= buf2.rep.elts();
2697 for (ind= 0; ind < repLengthBuf2; ind++)
2698 buf2p [ind]= gp [ind + lg];
2699
2700 buf2.normalize();
2701
2702 repLengthBuf2= buf2.rep.length();
2703
2704
2705 buf3.rep.SetLength((long) repLengthBuf2 + d);
2706 buf3p= buf3.rep.elts();
2707 buf2p= buf2.rep.elts();
2708 buf1p= buf1.rep.elts();
2709 for (ind= 0; ind < repLengthBuf1; ind++)
2710 buf3p [ind]= buf1p [ind];
2711 for (ind= repLengthBuf1; ind < d; ind++)
2712 buf3p [ind]= zzpZero;
2713 for (ind= 0; ind < repLengthBuf2; ind++)
2714 buf3p [ind + d]= buf2p [ind];
2715 buf3.normalize();
2716
2717 result += convertNTLzzpX2CF (buf3, x)*power (y, i);
2718 i++;
2719
2720
2721 lf= i*d;
2722 degfSubLf= degf - lf;
2723
2724 lg= d*(k-i);
2725 deggSubLg= degg - lg;
2726
2727 buf1p= buf1.rep.elts();
2728
2729 if (lg >= 0 && deggSubLg > 0)
2730 {
2731 if (repLengthBuf2 > degfSubLf + 1)
2732 degfSubLf= repLengthBuf2 - 1;
2733 tmp= tmin (repLengthBuf1, deggSubLg + 1);
2734 for (ind= 0; ind < tmp; ind++)
2735 gp [ind + lg] -= buf1p [ind];
2736 }
2737 if (lg < 0)
2738 break;
2739
2740 buf2p= buf2.rep.elts();
2741 if (degfSubLf >= 0)
2742 {
2743 for (ind= 0; ind < repLengthBuf2; ind++)
2744 fp [ind + lf] -= buf2p [ind];
2745 }
2746 }
2747
2748 return result;
2749}
2750#endif
2751
2752#if (!(HAVE_FLINT && __FLINT_RELEASE >= 20400))
2753CanonicalForm reverseSubstFq (const zz_pEX& F, int d, const Variable& alpha)
2754{
2755 Variable y= Variable (2);
2756 Variable x= Variable (1);
2757
2758 zz_pEX f= F;
2759 zz_pE *fp= f.rep.elts();
2760
2761 zz_pEX buf;
2762 zz_pE *bufp;
2764 int i= 0;
2765 int degf= deg(f);
2766 int k= 0;
2767 int degfSubK, repLength, j;
2768 while (degf >= k)
2769 {
2770 degfSubK= degf - k;
2771 if (degfSubK >= d)
2772 repLength= d;
2773 else
2774 repLength= degfSubK + 1;
2775
2776 buf.rep.SetLength ((long) repLength);
2777 bufp= buf.rep.elts();
2778 for (j= 0; j < repLength; j++)
2779 bufp [j]= fp [j + k];
2780 buf.normalize();
2781
2783 i++;
2784 k= d*i;
2785 }
2786
2787 return result;
2788}
2789#endif
2790
2791#ifndef HAVE_FLINT
2792CanonicalForm reverseSubstFp (const zz_pX& F, int d)
2793{
2794 Variable y= Variable (2);
2795 Variable x= Variable (1);
2796
2797 zz_pX f= F;
2798 zz_p *fp= f.rep.elts();
2799
2800 zz_pX buf;
2801 zz_p *bufp;
2803 int i= 0;
2804 int degf= deg(f);
2805 int k= 0;
2806 int degfSubK, repLength, j;
2807 while (degf >= k)
2808 {
2809 degfSubK= degf - k;
2810 if (degfSubK >= d)
2811 repLength= d;
2812 else
2813 repLength= degfSubK + 1;
2814
2815 buf.rep.SetLength ((long) repLength);
2816 bufp= buf.rep.elts();
2817 for (j= 0; j < repLength; j++)
2818 bufp [j]= fp [j + k];
2819 buf.normalize();
2820
2822 i++;
2823 k= d*i;
2824 }
2825
2826 return result;
2827}
2828
2829// assumes input to be reduced mod M and to be an element of Fp
2831mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
2833{
2834 int d1= degree (F, 1) + degree (G, 1) + 1;
2835 d1 /= 2;
2836 d1 += 1;
2837
2838 zz_pX F1, F2;
2839 kronSubReciproFp (F1, F2, F, d1);
2840 zz_pX G1, G2;
2841 kronSubReciproFp (G1, G2, G, d1);
2842
2843 int k= d1*degree (M);
2844 MulTrunc (F1, F1, G1, (long) k);
2845
2846 int degtailF= degree (tailcoeff (F), 1);
2847 int degtailG= degree (tailcoeff (G), 1);
2848 int taildegF= taildegree (F);
2849 int taildegG= taildegree (G);
2850 int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
2851
2852 reverse (F2, F2);
2853 reverse (G2, G2);
2854 MulTrunc (F2, F2, G2, b + 1);
2855 reverse (F2, F2, b);
2856
2857 int d2= tmax (deg (F2)/d1, deg (F1)/d1);
2858 return reverseSubstReciproFp (F1, F2, d1, d2);
2859}
2860
2861//Kronecker substitution
2863mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
2865{
2866 CanonicalForm A= F;
2867 CanonicalForm B= G;
2868
2869 int degAx= degree (A, 1);
2870 int degAy= degree (A, 2);
2871 int degBx= degree (B, 1);
2872 int degBy= degree (B, 2);
2873 int d1= degAx + 1 + degBx;
2874 int d2= tmax (degAy, degBy);
2875
2876 if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
2877 return mulMod2NTLFpReci (A, B, M);
2878
2879 zz_pX NTLA= kronSubFp (A, d1);
2880 zz_pX NTLB= kronSubFp (B, d1);
2881
2882 int k= d1*degree (M);
2883 MulTrunc (NTLA, NTLA, NTLB, (long) k);
2884
2885 A= reverseSubstFp (NTLA, d1);
2886
2887 return A;
2888}
2889#endif
2890
2891#if (!(HAVE_FLINT && __FLINT_RELEASE >= 20400))
2892// assumes input to be reduced mod M and to be an element of Fq not Fp
2894mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
2895 CanonicalForm& M, const Variable& alpha)
2896{
2897 int d1= degree (F, 1) + degree (G, 1) + 1;
2898 d1 /= 2;
2899 d1 += 1;
2900
2901 zz_pEX F1, F2;
2902 kronSubReciproFq (F1, F2, F, d1, alpha);
2903 zz_pEX G1, G2;
2904 kronSubReciproFq (G1, G2, G, d1, alpha);
2905
2906 int k= d1*degree (M);
2907 MulTrunc (F1, F1, G1, (long) k);
2908
2909 int degtailF= degree (tailcoeff (F), 1);
2910 int degtailG= degree (tailcoeff (G), 1);
2911 int taildegF= taildegree (F);
2912 int taildegG= taildegree (G);
2913 int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
2914
2915 reverse (F2, F2);
2916 reverse (G2, G2);
2917 MulTrunc (F2, F2, G2, b + 1);
2918 reverse (F2, F2, b);
2919
2920 int d2= tmax (deg (F2)/d1, deg (F1)/d1);
2921 return reverseSubstReciproFq (F1, F2, d1, d2, alpha);
2922}
2923#endif
2924
2925#ifdef HAVE_FLINT
2927mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
2928 CanonicalForm& M);
2929#endif
2930
2934{
2936 CanonicalForm A= F;
2937 CanonicalForm B= G;
2938
2940 {
2941#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
2942 nmod_poly_t FLINTmipo;
2943 convertFacCF2nmod_poly_t (FLINTmipo, getMipo (alpha));
2944
2945 fq_nmod_ctx_t fq_con;
2946 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
2947
2948 A= mulMod2FLINTFq (A, B, M, alpha, fq_con);
2949 nmod_poly_clear (FLINTmipo);
2951#else
2952 int degAx= degree (A, 1);
2953 int degAy= degree (A, 2);
2954 int degBx= degree (B, 1);
2955 int degBy= degree (B, 2);
2956 int d1= degAx + degBx + 1;
2957 int d2= tmax (degAy, degBy);
2959 {
2961 zz_p::init (getCharacteristic());
2962 }
2963 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2964 zz_pE::init (NTLMipo);
2965
2966 int degMipo= degree (getMipo (alpha));
2967 if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
2968 (2*degAy > degree (M)))
2969 return mulMod2NTLFqReci (A, B, M, alpha);
2970
2971 zz_pEX NTLA= kronSubFq (A, d1, alpha);
2972 zz_pEX NTLB= kronSubFq (B, d1, alpha);
2973
2974 int k= d1*degree (M);
2975
2976 MulTrunc (NTLA, NTLA, NTLB, (long) k);
2977
2978 A= reverseSubstFq (NTLA, d1, alpha);
2979#endif
2980 }
2981 else
2982 {
2983#ifdef HAVE_FLINT
2984 A= mulMod2FLINTFp (A, B, M);
2985#else
2986 A= mulMod2NTLFp (A, B, M);
2987#endif
2988 }
2989 return A;
2990}
2991
2993 const CanonicalForm& M)
2994{
2995 if (A.isZero() || B.isZero())
2996 return 0;
2997
2998 ASSERT (M.isUnivariate(), "M must be univariate");
2999
3000 CanonicalForm F= mod (A, M);
3001 CanonicalForm G= mod (B, M);
3002 if (F.inCoeffDomain())
3003 return G*F;
3004 if (G.inCoeffDomain())
3005 return F*G;
3006
3007 Variable y= M.mvar();
3008 int degF= degree (F, y);
3009 int degG= degree (G, y);
3010
3011 if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
3012 (F.level() == G.level()))
3013 {
3015 return mod (result, M);
3016 }
3017 else if (degF <= 1 && degG <= 1)
3018 {
3020 return mod (result, M);
3021 }
3022
3023 int sizeF= size (F);
3024 int sizeG= size (G);
3025
3026 int fallBackToNaive= 50;
3027 if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
3028 {
3029 if (sizeF < sizeG)
3030 return mod (G*F, M);
3031 else
3032 return mod (F*G, M);
3033 }
3034
3035#ifdef HAVE_FLINT
3036 if (getCharacteristic() == 0)
3037 return mulMod2FLINTQa (F, G, M);
3038#endif
3039
3041 (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
3042 return mulMod2NTLFq (F, G, M);
3043
3044 int m= (int) ceil (degree (M)/2.0);
3045 if (degF >= m || degG >= m)
3046 {
3047 CanonicalForm MLo= power (y, m);
3048 CanonicalForm MHi= power (y, degree (M) - m);
3049 CanonicalForm F0= mod (F, MLo);
3050 CanonicalForm F1= div (F, MLo);
3051 CanonicalForm G0= mod (G, MLo);
3052 CanonicalForm G1= div (G, MLo);
3053 CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
3054 CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
3055 CanonicalForm F0G0= mulMod2 (F0, G0, M);
3056 return F0G0 + MLo*(F0G1 + F1G0);
3057 }
3058 else
3059 {
3060 m= (int) ceil (tmax (degF, degG)/2.0);
3061 CanonicalForm yToM= power (y, m);
3062 CanonicalForm F0= mod (F, yToM);
3063 CanonicalForm F1= div (F, yToM);
3064 CanonicalForm G0= mod (G, yToM);
3065 CanonicalForm G1= div (G, yToM);
3066 CanonicalForm H00= mulMod2 (F0, G0, M);
3067 CanonicalForm H11= mulMod2 (F1, G1, M);
3068 CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
3069 return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
3070 }
3071 DEBOUTLN (cerr, "fatal end in mulMod2");
3072}
3073
3074// end bivariate polys
3075//**********************
3076// multivariate polys
3077
3079{
3080 CanonicalForm A= F;
3081 for (CFListIterator i= M; i.hasItem(); i++)
3082 A= mod (A, i.getItem());
3083 return A;
3084}
3085
3087 const CFList& MOD)
3088{
3089 if (A.isZero() || B.isZero())
3090 return 0;
3091
3092 if (MOD.length() == 1)
3093 return mulMod2 (A, B, MOD.getLast());
3094
3095 CanonicalForm M= MOD.getLast();
3096 CanonicalForm F= mod (A, M);
3097 CanonicalForm G= mod (B, M);
3098 if (F.inCoeffDomain())
3099 return G*F;
3100 if (G.inCoeffDomain())
3101 return F*G;
3102
3103 int sizeF= size (F);
3104 int sizeG= size (G);
3105
3106 if (sizeF / MOD.length() < 100 || sizeG / MOD.length() < 100)
3107 {
3108 if (sizeF < sizeG)
3109 return mod (G*F, MOD);
3110 else
3111 return mod (F*G, MOD);
3112 }
3113
3114 Variable y= M.mvar();
3115 int degF= degree (F, y);
3116 int degG= degree (G, y);
3117
3118 if ((degF <= 1 && F.level() <= M.level()) &&
3119 (degG <= 1 && G.level() <= M.level()))
3120 {
3121 CFList buf= MOD;
3122 buf.removeLast();
3123 if (degF == 1 && degG == 1)
3124 {
3125 CanonicalForm F0= mod (F, y);
3126 CanonicalForm F1= div (F, y);
3127 CanonicalForm G0= mod (G, y);
3128 CanonicalForm G1= div (G, y);
3129 if (degree (M) > 2)
3130 {
3131 CanonicalForm H00= mulMod (F0, G0, buf);
3132 CanonicalForm H11= mulMod (F1, G1, buf);
3133 CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
3134 return H11*y*y + (H01 - H00 - H11)*y + H00;
3135 }
3136 else //here degree (M) == 2
3137 {
3138 buf.append (y);
3139 CanonicalForm F0G1= mulMod (F0, G1, buf);
3140 CanonicalForm F1G0= mulMod (F1, G0, buf);
3141 CanonicalForm F0G0= mulMod (F0, G0, MOD);
3142 CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
3143 return result;
3144 }
3145 }
3146 else if (degF == 1 && degG == 0)
3147 return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
3148 else if (degF == 0 && degG == 1)
3149 return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
3150 else
3151 return mulMod (F, G, buf);
3152 }
3153 int m= (int) ceil (degree (M)/2.0);
3154 if (degF >= m || degG >= m)
3155 {
3156 CanonicalForm MLo= power (y, m);
3157 CanonicalForm MHi= power (y, degree (M) - m);
3158 CanonicalForm F0= mod (F, MLo);
3159 CanonicalForm F1= div (F, MLo);
3160 CanonicalForm G0= mod (G, MLo);
3161 CanonicalForm G1= div (G, MLo);
3162 CFList buf= MOD;
3163 buf.removeLast();
3164 buf.append (MHi);
3165 CanonicalForm F0G1= mulMod (F0, G1, buf);
3166 CanonicalForm F1G0= mulMod (F1, G0, buf);
3167 CanonicalForm F0G0= mulMod (F0, G0, MOD);
3168 return F0G0 + MLo*(F0G1 + F1G0);
3169 }
3170 else
3171 {
3172 m= (tmax(degF, degG)+1)/2;
3173 CanonicalForm yToM= power (y, m);
3174 CanonicalForm F0= mod (F, yToM);
3175 CanonicalForm F1= div (F, yToM);
3176 CanonicalForm G0= mod (G, yToM);
3177 CanonicalForm G1= div (G, yToM);
3178 CanonicalForm H00= mulMod (F0, G0, MOD);
3179 CanonicalForm H11= mulMod (F1, G1, MOD);
3180 CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
3181 return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
3182 }
3183 DEBOUTLN (cerr, "fatal end in mulMod");
3184}
3185
3187{
3188 if (L.isEmpty())
3189 return 1;
3190 int l= L.length();
3191 if (l == 1)
3192 return mod (L.getFirst(), M);
3193 else if (l == 2) {
3195 return result;
3196 }
3197 else
3198 {
3199 l /= 2;
3200 CFList tmp1, tmp2;
3201 CFListIterator i= L;
3203 for (int j= 1; j <= l; j++, i++)
3204 tmp1.append (i.getItem());
3205 tmp2= Difference (L, tmp1);
3206 buf1= prodMod (tmp1, M);
3207 buf2= prodMod (tmp2, M);
3209 return result;
3210 }
3211}
3212
3214{
3215 if (L.isEmpty())
3216 return 1;
3217 else if (L.length() == 1)
3218 return L.getFirst();
3219 else if (L.length() == 2)
3220 return mulMod (L.getFirst(), L.getLast(), M);
3221 else
3222 {
3223 int l= L.length()/2;
3224 CFListIterator i= L;
3225 CFList tmp1, tmp2;
3227 for (int j= 1; j <= l; j++, i++)
3228 tmp1.append (i.getItem());
3229 tmp2= Difference (L, tmp1);
3230 buf1= prodMod (tmp1, M);
3231 buf2= prodMod (tmp2, M);
3232 return mulMod (buf1, buf2, M);
3233 }
3234}
3235
3236// end multivariate polys
3237//***************************
3238// division
3239
3241{
3242 if (d == 0)
3243 return F;
3244 CanonicalForm A= F;
3245 Variable y= Variable (2);
3246 Variable x= Variable (1);
3247 if (degree (A, x) > 0)
3248 {
3249 A= swapvar (A, x, y);
3251 CFIterator i= A;
3252 while (d - i.exp() < 0)
3253 i++;
3254
3255 for (; i.hasTerms() && (d - i.exp() >= 0); i++)
3256 result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
3257 return result;
3258 }
3259 else
3260 return A*power (x, d);
3261}
3262
3264newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
3265{
3266 int l= ilog2(n);
3267
3268 CanonicalForm g= mod (F, M)[0] [0];
3269
3270 ASSERT (!g.isZero(), "expected a unit");
3271
3273
3274 if (!g.isOne())
3275 g = 1/g;
3276 Variable x= Variable (1);
3278 int exp= 0;
3279 if (n & 1)
3280 {
3281 result= g;
3282 exp= 1;
3283 }
3285
3286 for (int i= 1; i <= l; i++)
3287 {
3288 h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
3289 h= mod (h, power (x, (1 << i)) - 1);
3290 h= div (h, power (x, (1 << (i - 1))));
3291 h= mod (h, M);
3292 g -= power (x, (1 << (i - 1)))*
3293 mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
3294
3295 if (n & (1 << i))
3296 {
3297 if (exp)
3298 {
3299 h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
3300 h= mod (h, power (x, exp + (1 << i)) - 1);
3301 h= div (h, power (x, exp));
3302 h= mod (h, M);
3303 result -= power(x, exp)*mod (mulMod2 (g, h, M),
3304 power (x, (1 << i)));
3305 exp += (1 << i);
3306 }
3307 else
3308 {
3309 exp= (1 << i);
3310 result= g;
3311 }
3312 }
3313 }
3314
3315 return result;
3316}
3317
3320 M)
3321{
3322 ASSERT (getCharacteristic() > 0, "positive characteristic expected");
3323
3324 CanonicalForm A= mod (F, M);
3325 CanonicalForm B= mod (G, M);
3326
3327 Variable x= Variable (1);
3328 int degA= degree (A, x);
3329 int degB= degree (B, x);
3330 int m= degA - degB;
3331 if (m < 0)
3332 return 0;
3333
3334 Variable v;
3336 if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
3337 {
3339 divrem2 (A, B, Q, R, M);
3340 }
3341 else
3342 {
3343 if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
3344 {
3345 CanonicalForm R= reverse (A, degA);
3346 CanonicalForm revB= reverse (B, degB);
3347 revB= newtonInverse (revB, m + 1, M);
3348 Q= mulMod2 (R, revB, M);
3349 Q= mod (Q, power (x, m + 1));
3350 Q= reverse (Q, m);
3351 }
3352 else
3353 {
3354 Variable y= Variable (2);
3355#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
3356 nmod_poly_t FLINTmipo;
3357 fq_nmod_ctx_t fq_con;
3358
3359 nmod_poly_init (FLINTmipo, getCharacteristic());
3360 convertFacCF2nmod_poly_t (FLINTmipo, M);
3361
3362 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
3363
3364
3365 fq_nmod_poly_t FLINTA, FLINTB;
3368
3369 fq_nmod_poly_divrem (FLINTA, FLINTB, FLINTA, FLINTB, fq_con);
3370
3371 Q= convertFq_nmod_poly_t2FacCF (FLINTA, x, y, fq_con);
3372
3373 fq_nmod_poly_clear (FLINTA, fq_con);
3374 fq_nmod_poly_clear (FLINTB, fq_con);
3375 nmod_poly_clear (FLINTmipo);
3377#else
3378 bool zz_pEbak= zz_pE::initialized();
3379 zz_pEBak bak;
3380 if (zz_pEbak)
3381 bak.save();
3382 zz_pX mipo= convertFacCF2NTLzzpX (M);
3383 zz_pEX NTLA, NTLB;
3384 NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
3385 NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
3386 div (NTLA, NTLA, NTLB);
3387 Q= convertNTLzz_pEX2CF (NTLA, x, y);
3388 if (zz_pEbak)
3389 bak.restore();
3390#endif
3391 }
3392 }
3393
3394 return Q;
3395}
3396
3397void
3399 CanonicalForm& R, const CanonicalForm& M)
3400{
3401 CanonicalForm A= mod (F, M);
3402 CanonicalForm B= mod (G, M);
3403 Variable x= Variable (1);
3404 int degA= degree (A, x);
3405 int degB= degree (B, x);
3406 int m= degA - degB;
3407
3408 if (m < 0)
3409 {
3410 R= A;
3411 Q= 0;
3412 return;
3413 }
3414
3415 Variable v;
3416 if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
3417 {
3418 divrem2 (A, B, Q, R, M);
3419 }
3420 else
3421 {
3422 if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
3423 {
3424 R= reverse (A, degA);
3425
3426 CanonicalForm revB= reverse (B, degB);
3427 revB= newtonInverse (revB, m + 1, M);
3428 Q= mulMod2 (R, revB, M);
3429
3430 Q= mod (Q, power (x, m + 1));
3431 Q= reverse (Q, m);
3432
3433 R= A - mulMod2 (Q, B, M);
3434 }
3435 else
3436 {
3437 Variable y= Variable (2);
3438#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
3439 nmod_poly_t FLINTmipo;
3440 fq_nmod_ctx_t fq_con;
3441
3442 nmod_poly_init (FLINTmipo, getCharacteristic());
3443 convertFacCF2nmod_poly_t (FLINTmipo, M);
3444
3445 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
3446
3447 fq_nmod_poly_t FLINTA, FLINTB;
3450
3451 fq_nmod_poly_divrem (FLINTA, FLINTB, FLINTA, FLINTB, fq_con);
3452
3453 Q= convertFq_nmod_poly_t2FacCF (FLINTA, x, y, fq_con);
3454 R= convertFq_nmod_poly_t2FacCF (FLINTB, x, y, fq_con);
3455
3456 fq_nmod_poly_clear (FLINTA, fq_con);
3457 fq_nmod_poly_clear (FLINTB, fq_con);
3458 nmod_poly_clear (FLINTmipo);
3460#else
3461 zz_pX mipo= convertFacCF2NTLzzpX (M);
3462 zz_pEX NTLA, NTLB;
3463 NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
3464 NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
3465 zz_pEX NTLQ, NTLR;
3466 DivRem (NTLQ, NTLR, NTLA, NTLB);
3467 Q= convertNTLzz_pEX2CF (NTLQ, x, y);
3468 R= convertNTLzz_pEX2CF (NTLR, x, y);
3469#endif
3470 }
3471 }
3472}
3473
3474static inline
3475CFList split (const CanonicalForm& F, const int m, const Variable& x)
3476{
3477 CanonicalForm A= F;
3478 CanonicalForm buf= 0;
3479 bool swap= false;
3480 if (degree (A, x) <= 0)
3481 return CFList(A);
3482 else if (x.level() != A.level())
3483 {
3484 swap= true;
3485 A= swapvar (A, x, A.mvar());
3486 }
3487
3488 int j= (int) floor ((double) degree (A)/ m);
3489 CFList result;
3490 CFIterator i= A;
3491 for (; j >= 0; j--)
3492 {
3493 while (i.hasTerms() && i.exp() - j*m >= 0)
3494 {
3495 if (swap)
3496 buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
3497 else
3498 buf += i.coeff()*power (x, i.exp() - j*m);
3499 i++;
3500 }
3501 if (swap)
3502 result.append (swapvar (buf, x, F.mvar()));
3503 else
3504 result.append (buf);
3505 buf= 0;
3506 }
3507 return result;
3508}
3509
3510static inline
3511void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
3512 CanonicalForm& R, const CFList& M);
3513
3514static inline
3516 CanonicalForm& R, const CFList& M)
3517{
3518 CanonicalForm A= mod (F, M);
3519 CanonicalForm B= mod (G, M);
3520 Variable x= Variable (1);
3521 int degB= degree (B, x);
3522 int degA= degree (A, x);
3523 if (degA < degB)
3524 {
3525 Q= 0;
3526 R= A;
3527 return;
3528 }
3529 if (degB < 1)
3530 {
3531 divrem (A, B, Q, R);
3532 Q= mod (Q, M);
3533 R= mod (R, M);
3534 return;
3535 }
3536 int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
3537 ASSERT (4*m >= degA, "expected degree (F, 1) < 2*degree (G, 1)");
3538 CFList splitA= split (A, m, x);
3539 if (splitA.length() == 3)
3540 splitA.insert (0);
3541 if (splitA.length() == 2)
3542 {
3543 splitA.insert (0);
3544 splitA.insert (0);
3545 }
3546 if (splitA.length() == 1)
3547 {
3548 splitA.insert (0);
3549 splitA.insert (0);
3550 splitA.insert (0);
3551 }
3552
3553 CanonicalForm xToM= power (x, m);
3554
3555 CFListIterator i= splitA;
3556 CanonicalForm H= i.getItem();
3557 i++;
3558 H *= xToM;
3559 H += i.getItem();
3560 i++;
3561 H *= xToM;
3562 H += i.getItem();
3563 i++;
3564
3565 divrem32 (H, B, Q, R, M);
3566
3567 CFList splitR= split (R, m, x);
3568 if (splitR.length() == 1)
3569 splitR.insert (0);
3570
3571 H= splitR.getFirst();
3572 H *= xToM;
3573 H += splitR.getLast();
3574 H *= xToM;
3575 H += i.getItem();
3576
3577 CanonicalForm bufQ;
3578 divrem32 (H, B, bufQ, R, M);
3579
3580 Q *= xToM;
3581 Q += bufQ;
3582 return;
3583}
3584
3585static inline
3587 CanonicalForm& R, const CFList& M)
3588{
3589 CanonicalForm A= mod (F, M);
3590 CanonicalForm B= mod (G, M);
3591 Variable x= Variable (1);
3592 int degB= degree (B, x);
3593 int degA= degree (A, x);
3594 if (degA < degB)
3595 {
3596 Q= 0;
3597 R= A;
3598 return;
3599 }
3600 if (degB < 1)
3601 {
3602 divrem (A, B, Q, R);
3603 Q= mod (Q, M);
3604 R= mod (R, M);
3605 return;
3606 }
3607 int m= (int) ceil ((double) (degB + 1)/ 2.0);
3608 ASSERT (3*m > degA, "expected degree (F, 1) < 3*degree (G, 1)");
3609 CFList splitA= split (A, m, x);
3610 CFList splitB= split (B, m, x);
3611
3612 if (splitA.length() == 2)
3613 {
3614 splitA.insert (0);
3615 }
3616 if (splitA.length() == 1)
3617 {
3618 splitA.insert (0);
3619 splitA.insert (0);
3620 }
3621 CanonicalForm xToM= power (x, m);
3622
3624 CFListIterator i= splitA;
3625 i++;
3626
3627 if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
3628 {
3629 H= splitA.getFirst()*xToM + i.getItem();
3630 divrem21 (H, splitB.getFirst(), Q, R, M);
3631 }
3632 else
3633 {
3634 R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
3635 splitB.getFirst()*xToM;
3636 Q= xToM - 1;
3637 }
3638
3639 H= mulMod (Q, splitB.getLast(), M);
3640
3641 R= R*xToM + splitA.getLast() - H;
3642
3643 while (degree (R, x) >= degB)
3644 {
3645 xToM= power (x, degree (R, x) - degB);
3646 Q += LC (R, x)*xToM;
3647 R -= mulMod (LC (R, x), B, M)*xToM;
3648 Q= mod (Q, M);
3649 R= mod (R, M);
3650 }
3651
3652 return;
3653}
3654
3656 CanonicalForm& R, const CanonicalForm& M)
3657{
3658 CanonicalForm A= mod (F, M);
3659 CanonicalForm B= mod (G, M);
3660
3661 if (B.inCoeffDomain())
3662 {
3663 divrem (A, B, Q, R);
3664 return;
3665 }
3666 if (A.inCoeffDomain() && !B.inCoeffDomain())
3667 {
3668 Q= 0;
3669 R= A;
3670 return;
3671 }
3672
3673 if (B.level() < A.level())
3674 {
3675 divrem (A, B, Q, R);
3676 return;
3677 }
3678 if (A.level() > B.level())
3679 {
3680 R= A;
3681 Q= 0;
3682 return;
3683 }
3684 if (B.level() == 1 && B.isUnivariate())
3685 {
3686 divrem (A, B, Q, R);
3687 return;
3688 }
3689
3690 Variable x= Variable (1);
3691 int degB= degree (B, x);
3692 if (degB > degree (A, x))
3693 {
3694 Q= 0;
3695 R= A;
3696 return;
3697 }
3698
3699 CFList splitA= split (A, degB, x);
3700
3701 CanonicalForm xToDegB= power (x, degB);
3702 CanonicalForm H, bufQ;
3703 Q= 0;
3704 CFListIterator i= splitA;
3705 H= i.getItem()*xToDegB;
3706 i++;
3707 H += i.getItem();
3708 CFList buf;
3709 while (i.hasItem())
3710 {
3711 buf= CFList (M);
3712 divrem21 (H, B, bufQ, R, buf);
3713 i++;
3714 if (i.hasItem())
3715 H= R*xToDegB + i.getItem();
3716 Q *= xToDegB;
3717 Q += bufQ;
3718 }
3719 return;
3720}
3721
3723 CanonicalForm& R, const CFList& MOD)
3724{
3725 CanonicalForm A= mod (F, MOD);
3726 CanonicalForm B= mod (G, MOD);
3727 Variable x= Variable (1);
3728 int degB= degree (B, x);
3729 if (degB > degree (A, x))
3730 {
3731 Q= 0;
3732 R= A;
3733 return;
3734 }
3735
3736 if (degB <= 0)
3737 {
3738 divrem (A, B, Q, R);
3739 Q= mod (Q, MOD);
3740 R= mod (R, MOD);
3741 return;
3742 }
3743 CFList splitA= split (A, degB, x);
3744
3745 CanonicalForm xToDegB= power (x, degB);
3746 CanonicalForm H, bufQ;
3747 Q= 0;
3748 CFListIterator i= splitA;
3749 H= i.getItem()*xToDegB;
3750 i++;
3751 H += i.getItem();
3752 while (i.hasItem())
3753 {
3754 divrem21 (H, B, bufQ, R, MOD);
3755 i++;
3756 if (i.hasItem())
3757 H= R*xToDegB + i.getItem();
3758 Q *= xToDegB;
3759 Q += bufQ;
3760 }
3761 return;
3762}
3763
3764bool
3766{
3767 if (B.isZero())
3768 return true;
3769 if (A.isZero())
3770 return false;
3772 return fdivides (A, B);
3773 int p= getCharacteristic();
3774 if (A.inCoeffDomain() || B.inCoeffDomain())
3775 {
3776 if (A.inCoeffDomain())
3777 return true;
3778 else
3779 return false;
3780 }
3781 if (p > 0)
3782 {
3783#if (!defined(HAVE_FLINT) || __FLINT_RELEASE < 20400)
3784 if (fac_NTL_char != p)
3785 {
3786 fac_NTL_char= p;
3787 zz_p::init (p);
3788 }
3789#endif
3792 {
3793#if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
3794 nmod_poly_t FLINTmipo;
3795 fq_nmod_ctx_t fq_con;
3796
3797 nmod_poly_init (FLINTmipo, getCharacteristic());
3798 convertFacCF2nmod_poly_t (FLINTmipo, getMipo (alpha));
3799
3800 fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
3801
3802 fq_nmod_poly_t FLINTA, FLINTB;
3805 int result= fq_nmod_poly_divides (FLINTA, FLINTB, FLINTA, fq_con);
3806 fq_nmod_poly_clear (FLINTA, fq_con);
3807 fq_nmod_poly_clear (FLINTB, fq_con);
3808 nmod_poly_clear (FLINTmipo);
3810 return result;
3811#else
3812 zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
3813 zz_pE::init (NTLMipo);
3814 zz_pEX NTLA= convertFacCF2NTLzz_pEX (A, NTLMipo);
3815 zz_pEX NTLB= convertFacCF2NTLzz_pEX (B, NTLMipo);
3816 return divide (NTLB, NTLA);
3817#endif
3818 }
3819#ifdef HAVE_FLINT
3820 nmod_poly_t FLINTA, FLINTB;
3821 convertFacCF2nmod_poly_t (FLINTA, A);
3822 convertFacCF2nmod_poly_t (FLINTB, B);
3823 nmod_poly_divrem (FLINTB, FLINTA, FLINTB, FLINTA);
3824 bool result= nmod_poly_is_zero (FLINTA);
3825 nmod_poly_clear (FLINTA);
3826 nmod_poly_clear (FLINTB);
3827 return result;
3828#else
3829 zz_pX NTLA= convertFacCF2NTLzzpX (A);
3830 zz_pX NTLB= convertFacCF2NTLzzpX (B);
3831 return divide (NTLB, NTLA);
3832#endif
3833 }
3834#ifdef HAVE_FLINT
3836 bool isRat= isOn (SW_RATIONAL);
3837 if (!isRat)
3838 On (SW_RATIONAL);
3839 if (!hasFirstAlgVar (A, alpha) && !hasFirstAlgVar (B, alpha))
3840 {
3841 fmpq_poly_t FLINTA,FLINTB;
3842 convertFacCF2Fmpq_poly_t (FLINTA, A);
3843 convertFacCF2Fmpq_poly_t (FLINTB, B);
3844 fmpq_poly_rem (FLINTA, FLINTB, FLINTA);
3845 bool result= fmpq_poly_is_zero (FLINTA);
3846 fmpq_poly_clear (FLINTA);
3847 fmpq_poly_clear (FLINTB);
3848 if (!isRat)
3849 Off (SW_RATIONAL);
3850 return result;
3851 }
3852 CanonicalForm Q, R;
3853 newtonDivrem (B, A, Q, R);
3854 if (!isRat)
3855 Off (SW_RATIONAL);
3856 return R.isZero();
3857#else
3858 bool isRat= isOn (SW_RATIONAL);
3859 if (!isRat)
3860 On (SW_RATIONAL);
3861 bool result= fdivides (A, B);
3862 if (!isRat)
3863 Off (SW_RATIONAL);
3864 return result; //maybe NTL?
3865#endif
3866}
3867
3868// end division
3869
3870#else
3872mulNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
3873{ return F*G; }
3874#endif
CanonicalForm convertFq_poly_t2FacCF(const fq_poly_t p, const Variable &x, const Variable &alpha, const fq_ctx_t ctx)
conversion of a FLINT poly over Fq (for non-word size p) to a CanonicalForm with alg....
void convertFacCF2Fq_t(fq_t result, const CanonicalForm &f, const fq_ctx_t ctx)
conversion of a factory element of F_q (for non-word size p) to a FLINT fq_t
CanonicalForm convertFq_nmod_poly_t2FacCF(const fq_nmod_poly_t p, const Variable &x, const Variable &alpha, const fq_nmod_ctx_t ctx)
conversion of a FLINT poly over Fq to a CanonicalForm with alg. variable alpha and polynomial variabl...
CanonicalForm convertFq_t2FacCF(const fq_t poly, const Variable &alpha)
conversion of a FLINT element of F_q with non-word size p to a CanonicalForm with alg....
CanonicalForm convertFmpq_poly_t2FacCF(const fmpq_poly_t p, const Variable &x)
conversion of a FLINT poly over Q to CanonicalForm
CanonicalForm convertFmpz_mod_poly_t2FacCF(const fmpz_mod_poly_t poly, const Variable &x, const modpk &b)
conversion of a FLINT poly over Z/p (for non word size p) to a CanonicalForm over Z
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
void convertFacCF2Fmpz_mod_poly_t(fmpz_mod_poly_t result, const CanonicalForm &f, const fmpz_t p)
conversion of a factory univariate poly over Z to a FLINT poly over Z/p (for non word size p)
void convertFacCF2Fq_nmod_poly_t(fq_nmod_poly_t result, const CanonicalForm &f, const fq_nmod_ctx_t ctx)
conversion of a factory univariate poly over F_q to a FLINT fq_nmod_poly_t
CanonicalForm convertFmpz_poly_t2FacCF(const fmpz_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z to CanonicalForm
void convertFacCF2Fmpq_poly_t(fmpq_poly_t result, const CanonicalForm &f)
conversion of a factory univariate polynomials over Q to fmpq_poly_t
void convertFacCF2Fmpz_poly_t(fmpz_poly_t result, const CanonicalForm &f)
conversion of a factory univariate polynomial over Z to a fmpz_poly_t
void convertCF2initFmpz(fmpz_t result, const CanonicalForm &f)
conversion of a factory integer to fmpz_t(init.)
void convertFacCF2Fq_poly_t(fq_poly_t result, const CanonicalForm &f, const fq_ctx_t ctx)
conversion of a factory univariate poly over F_q (for non-word size p) to a FLINT fq_poly_t
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
zz_pEX convertFacCF2NTLzz_pEX(const CanonicalForm &f, const zz_pX &mipo)
CanonicalForm convertNTLzz_pEX2CF(const zz_pEX &f, const Variable &x, const Variable &alpha)
ZZ_pEX convertFacCF2NTLZZ_pEX(const CanonicalForm &f, const ZZ_pX &mipo)
CanonicalForm in Z_p(a)[X] to NTL ZZ_pEX.
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
CanonicalForm convertNTLZZpX2CF(const ZZ_pX &poly, const Variable &x)
NAME: convertNTLZZpX2CF.
CanonicalForm convertNTLZZX2CF(const ZZX &polynom, const Variable &x)
CanonicalForm convertNTLZZ_pEX2CF(const ZZ_pEX &f, const Variable &x, const Variable &alpha)
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
ZZ_pX convertFacCF2NTLZZpX(const CanonicalForm &f)
NAME: convertFacCF2NTLZZpX.
Definition NTLconvert.cc:64
VAR long fac_NTL_char
Definition NTLconvert.cc:46
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Conversion to and from NTL.
#define swap(_i, _j)
bool isOn(int sw)
switches
void On(int sw)
switches
void Off(int sw)
switches
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Header for factory's main class CanonicalForm.
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition cf_ops.cc:600
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm div(const CanonicalForm &, const CanonicalForm &)
int ilog2(const CanonicalForm &a)
CanonicalForm tailcoeff(const CanonicalForm &f)
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
int taildegree(const CanonicalForm &f)
int degree(const CanonicalForm &f)
bool hasFirstAlgVar(const CanonicalForm &f, Variable &a)
check if poly f contains an algebraic variable a
Definition cf_ops.cc:679
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition cf_ops.cc:168
ListIterator< CanonicalForm > CFListIterator
CanonicalForm den(const CanonicalForm &f)
List< CanonicalForm > CFList
CanonicalForm LC(const CanonicalForm &f)
int FACTORY_PUBLIC getCharacteristic()
Definition cf_char.cc:70
int l
Definition cfEzgcd.cc:100
int m
Definition cfEzgcd.cc:128
int i
Definition cfEzgcd.cc:132
int k
Definition cfEzgcd.cc:99
Variable x
Definition cfModGcd.cc:4090
int p
Definition cfModGcd.cc:4086
g
Definition cfModGcd.cc:4098
CanonicalForm fp
Definition cfModGcd.cc:4110
CanonicalForm cd(bCommonDen(FF))
Definition cfModGcd.cc:4097
CanonicalForm gp
Definition cfModGcd.cc:4110
CanonicalForm b
Definition cfModGcd.cc:4111
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
declarations of higher level algorithms.
#define ASSERT(expression, message)
Definition cf_assert.h:99
static const int SW_RATIONAL
set to 1 for computations over Q
Definition cf_defs.h:31
#define GaloisFieldDomain
Definition cf_defs.h:18
Iterators for CanonicalForm's.
FILE * f
Definition checklibs.c:9
static int gettype()
Definition cf_factory.h:28
class to iterate through CanonicalForm's
Definition cf_iter.h:44
factory's main class
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
bool inCoeffDomain() const
int level() const
level() returns the level of CO.
bool inBaseDomain() const
bool isUnivariate() const
T getFirst() const
int length() const
T getLast() const
void insert(const T &)
int isEmpty() const
factory's class for variables
Definition factory.h:127
class to do operations mod p^k for int's p and k
Definition fac_util.h:23
functions to print debug output
#define DEBOUTLN(stream, objects)
Definition debug.h:49
Variable alpha
return result
const CanonicalForm int const CFList const Variable & y
Definition facAbsFact.cc:53
CanonicalForm H
Definition facAbsFact.cc:60
int degg
Definition facAlgExt.cc:64
CanonicalForm mipo
Definition facAlgExt.cc:57
CanonicalForm divide(const CanonicalForm &ff, const CanonicalForm &f, const CFList &as)
b *CanonicalForm B
Definition facBivar.cc:52
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
CFList tmp1
Definition facFqBivar.cc:75
CanonicalForm buf2
Definition facFqBivar.cc:76
CFList tmp2
Definition facFqBivar.cc:75
CanonicalForm buf1
Definition facFqBivar.cc:76
fq_nmod_ctx_t fq_con
Definition facHensel.cc:99
int j
Definition facHensel.cc:110
fq_nmod_ctx_clear(fq_con)
nmod_poly_init(FLINTmipo, getCharacteristic())
fq_nmod_ctx_init_modulus(fq_con, FLINTmipo, "Z")
fq_nmod_poly_init(prod, fq_con)
convertFacCF2nmod_poly_t(FLINTmipo, M)
nmod_poly_clear(FLINTmipo)
fq_nmod_poly_clear(prod, fq_con)
CanonicalForm mod(const CanonicalForm &F, const CFList &M)
reduce F modulo elements in M.
Definition facMul.cc:3078
CanonicalForm uniReverse(const CanonicalForm &F, int d, const Variable &x)
Definition facMul.cc:280
void newtonDivrem(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Q, CanonicalForm &R)
division with remainder of univariate polynomials over Q and Q(a) using Newton inversion,...
Definition facMul.cc:352
void kronSubFq(fq_nmod_poly_t result, const CanonicalForm &A, int d, const fq_nmod_ctx_t fq_con)
Definition facMul.cc:1277
CanonicalForm mulNTL(const CanonicalForm &F, const CanonicalForm &G, const modpk &b)
multiplication of univariate polys using FLINT/NTL over F_p, F_q, Z/p^k, Z/p^k[t]/(f),...
Definition facMul.cc:417
void divrem(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Q, CanonicalForm &R, const CFList &MOD)
division with remainder of F by G wrt Variable (1) modulo MOD. Uses an algorithm based on Burnikel,...
Definition facMul.cc:3722
bool uniFdivides(const CanonicalForm &A, const CanonicalForm &B)
divisibility test for univariate polys
Definition facMul.cc:3765
CanonicalForm divFLINTQ(const CanonicalForm &F, const CanonicalForm &G)
Definition facMul.cc:180
void kronSubQa(fmpz_poly_t result, const CanonicalForm &A, int d)
Definition facMul.cc:52
CanonicalForm reverseSubstFp(const nmod_poly_t F, int d)
Definition facMul.cc:2060
static CFList split(const CanonicalForm &F, const int m, const Variable &x)
Definition facMul.cc:3475
static void divrem32(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Q, CanonicalForm &R, const CFList &M)
Definition facMul.cc:3586
CanonicalForm mulMod2FLINTQ(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M)
Definition facMul.cc:2278
CanonicalForm reverse(const CanonicalForm &F, int d)
Definition facMul.cc:3240
CanonicalForm mulMod2(const CanonicalForm &A, const CanonicalForm &B, const CanonicalForm &M)
Karatsuba style modular multiplication for bivariate polynomials.
Definition facMul.cc:2992
CanonicalForm mulMod2FLINTFqReci(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M, const Variable &alpha, const fq_nmod_ctx_t fq_con)
Definition facMul.cc:2166
CanonicalForm mulFLINTQ(const CanonicalForm &F, const CanonicalForm &G)
Definition facMul.cc:139
CanonicalForm mulMod2FLINTFpReci(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M)
Definition facMul.cc:2096
CanonicalForm mulMod2FLINTFq(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M, const Variable &alpha, const fq_nmod_ctx_t fq_con)
Definition facMul.cc:2208
CanonicalForm reverseSubstReciproFq(const fq_nmod_poly_t F, const fq_nmod_poly_t G, int d, int k, const Variable &alpha, const fq_nmod_ctx_t fq_con)
Definition facMul.cc:1787
CanonicalForm modFLINTQ(const CanonicalForm &F, const CanonicalForm &G)
Definition facMul.cc:198
void kronSubReciproQ(fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm &A, int d)
Definition facMul.cc:1480
void kronSubReciproFq(fq_nmod_poly_t subA1, fq_nmod_poly_t subA2, const CanonicalForm &A, int d, const fq_nmod_ctx_t fq_con)
Definition facMul.cc:1435
CanonicalForm reverseSubstQ(const fmpz_poly_t F, int d)
Definition facMul.cc:1505
CanonicalForm mulMod2FLINTQReci(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M)
Definition facMul.cc:2241
CanonicalForm reverseSubstFq(const fq_nmod_poly_t F, int d, const Variable &alpha, const fq_nmod_ctx_t fq_con)
Definition facMul.cc:2025
CanonicalForm mulMod(const CanonicalForm &A, const CanonicalForm &B, const CFList &MOD)
Karatsuba style modular multiplication for multivariate polynomials.
Definition facMul.cc:3086
CanonicalForm mulFLINTQTrunc(const CanonicalForm &F, const CanonicalForm &G, int m)
Definition facMul.cc:247
CanonicalForm mulFLINTQa(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha)
Definition facMul.cc:109
CanonicalForm reverseSubstReciproQ(const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
Definition facMul.cc:1891
static void divrem21(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Q, CanonicalForm &R, const CFList &M)
Definition facMul.cc:3515
CanonicalForm newtonInverse(const CanonicalForm &F, const int n, const Variable &x)
Definition facMul.cc:297
void newtonDiv(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Q)
Definition facMul.cc:386
void divrem2(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Q, CanonicalForm &R, const CanonicalForm &M)
division with remainder of F by G wrt Variable (1) modulo M. Uses an algorithm based on Burnikel,...
Definition facMul.cc:3655
CanonicalForm reverseSubstQa(const fmpz_poly_t F, int d, const Variable &x, const Variable &alpha, const CanonicalForm &den)
Definition facMul.cc:72
void kronSubReciproFp(nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm &A, int d)
Definition facMul.cc:1394
CanonicalForm divNTL(const CanonicalForm &F, const CanonicalForm &G, const modpk &b)
division of univariate polys using FLINT/NTL over F_p, F_q, Z/p^k, Z/p^k[t]/(f), Z,...
Definition facMul.cc:942
CanonicalForm mulFLINTQaTrunc(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, int m)
Definition facMul.cc:216
CanonicalForm modNTL(const CanonicalForm &F, const CanonicalForm &G, const modpk &b)
mod of univariate polys using FLINT/NTL over F_p, F_q, Z/p^k, Z/p^k[t]/(f), Z, Q, Q(a),...
Definition facMul.cc:737
CanonicalForm prodMod(const CFList &L, const CanonicalForm &M)
product of all elements in L modulo M via divide-and-conquer.
Definition facMul.cc:3186
CanonicalForm reverseSubstReciproFp(const nmod_poly_t F, const nmod_poly_t G, int d, int k)
Definition facMul.cc:1668
void kronSubFp(nmod_poly_t result, const CanonicalForm &A, int d)
Definition facMul.cc:1254
CanonicalForm mulMod2FLINTFp(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M)
Definition facMul.cc:2134
CanonicalForm mulMod2NTLFq(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M)
Definition facMul.cc:2932
CanonicalForm mulMod2FLINTQa(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &M)
Definition facMul.cc:2338
This file defines functions for fast multiplication and division with remainder.
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition variable.cc:207
some useful template functions.
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
template List< Variable > Difference(const List< Variable > &, const List< Variable > &)
STATIC_VAR TreeM * G
Definition janet.cc:31
STATIC_VAR Poly * h
Definition janet.cc:971
void rem(unsigned long *a, unsigned long *q, unsigned long p, int &dega, int degq)
Definition minpoly.cc:572
gmp_float exp(const gmp_float &a)
The main handler for Singular numbers which are suitable for Singular polynomials.
int status int void * buf
Definition si_signals.h:69
#define R
Definition sirandom.c:27
#define A
Definition sirandom.c:24
#define M
Definition sirandom.c:25
#define Q
Definition sirandom.c:26
int F1(int a1, int &r1)
F1.
void F2(int a2, int &r2)
F2.
bool getReduce(const Variable &alpha)
Definition variable.cc:232