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