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