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