source: git/factory/facMul.cc @ b78a13

spielwiese
Last change on this file since b78a13 was b78a13, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: faster reversion of Kronecker substitution
  • Property mode set to 100644
File size: 64.5 KB
Line 
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 * @author Martin Lee
10 *
11 **/
12/*****************************************************************************/
13
14#include "debug.h"
15
16#include "canonicalform.h"
17#include "facMul.h"
18#include "algext.h"
19#include "cf_util.h"
20#include "templates/ftmpl_functions.h"
21
22#ifdef HAVE_NTL
23#include <NTL/lzz_pEX.h>
24#include "NTLconvert.h"
25
26#ifdef HAVE_FLINT
27#include "FLINTconvert.h"
28#endif
29
30// univariate polys
31
32#ifdef HAVE_FLINT
33void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
34{
35  int degAy= degree (A);
36  fmpz_poly_init2 (result, d*(degAy + 1));
37  _fmpz_poly_set_length (result, d*(degAy + 1));
38  CFIterator j;
39  for (CFIterator i= A; i.hasTerms(); i++)
40  {
41    if (i.coeff().inBaseDomain())
42      convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
43    else
44      for (j= i.coeff(); j.hasTerms(); j++)
45        convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
46                        j.coeff());
47  }
48  _fmpz_poly_normalise(result);
49}
50
51
52CanonicalForm
53reverseSubstQa (const fmpz_poly_t F, int d, const Variable& alpha,
54                const CanonicalForm& den)
55{
56  Variable x= Variable (1);
57
58  CanonicalForm result= 0;
59  int i= 0;
60  int degf= fmpz_poly_degree (F);
61  int k= 0;
62  int degfSubK;
63  int repLength, j;
64  CanonicalForm coeff, ff;
65  fmpz* tmp;
66  while (degf >= k)
67  {
68    coeff= 0;
69    degfSubK= degf - k;
70    if (degfSubK >= d)
71      repLength= d;
72    else
73      repLength= degfSubK + 1;
74
75    for (j= 0; j < repLength; j++)
76    {
77      tmp= fmpz_poly_get_coeff_ptr (F, j+k);
78      if (!fmpz_is_zero (tmp))
79      {
80        ff= convertFmpz2CF (tmp);
81        coeff += ff*power (alpha, j); //TODO faster reduction mod alpha
82      }
83    }
84    result += coeff*power (x, i);
85    i++;
86    k= d*i;
87  }
88  result /= den;
89  return result;
90}
91
92CanonicalForm
93mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
94            const Variable& alpha)
95{
96  CanonicalForm A= F;
97  CanonicalForm B= G;
98
99  CanonicalForm denA= bCommonDen (A);
100  CanonicalForm denB= bCommonDen (B);
101
102  A *= denA;
103  B *= denB;
104  int degAa= degree (A, alpha);
105  int degBa= degree (B, alpha);
106  int d= degAa + 1 + degBa;
107
108  fmpz_poly_t FLINTA,FLINTB;
109  fmpz_poly_init (FLINTA);
110  fmpz_poly_init (FLINTB);
111  kronSub (FLINTA, A, d);
112  kronSub (FLINTB, B, d);
113
114  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
115
116  denA *= denB;
117  A= reverseSubstQa (FLINTA, d, alpha, denA);
118
119  fmpz_poly_clear (FLINTA);
120  fmpz_poly_clear (FLINTB);
121  return A;
122}
123
124CanonicalForm
125mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
126{
127  CanonicalForm A= F;
128  CanonicalForm B= G;
129
130  CanonicalForm denA= bCommonDen (A);
131  CanonicalForm denB= bCommonDen (B);
132
133  A *= denA;
134  B *= denB;
135  fmpz_poly_t FLINTA,FLINTB;
136  convertFacCF2Fmpz_poly_t (FLINTA, A);
137  convertFacCF2Fmpz_poly_t (FLINTB, B);
138  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
139  denA *= denB;
140  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
141  A /= denA;
142  fmpz_poly_clear (FLINTA);
143  fmpz_poly_clear (FLINTB);
144
145  return A;
146}
147
148/*CanonicalForm
149mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
150{
151  CanonicalForm A= F;
152  CanonicalForm B= G;
153
154  fmpq_poly_t FLINTA,FLINTB;
155  convertFacCF2Fmpq_poly_t (FLINTA, A);
156  convertFacCF2Fmpq_poly_t (FLINTB, B);
157
158  fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
159  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
160  fmpq_poly_clear (FLINTA);
161  fmpq_poly_clear (FLINTB);
162  return A;
163}*/
164
165CanonicalForm
166divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
167{
168  CanonicalForm A= F;
169  CanonicalForm B= G;
170
171  fmpq_poly_t FLINTA,FLINTB;
172  convertFacCF2Fmpq_poly_t (FLINTA, A);
173  convertFacCF2Fmpq_poly_t (FLINTB, B);
174
175  fmpq_poly_div (FLINTA, FLINTA, FLINTB);
176  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
177
178  fmpq_poly_clear (FLINTA);
179  fmpq_poly_clear (FLINTB);
180  return A;
181}
182
183CanonicalForm
184modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
185{
186  CanonicalForm A= F;
187  CanonicalForm B= G;
188
189  fmpq_poly_t FLINTA,FLINTB;
190  convertFacCF2Fmpq_poly_t (FLINTA, A);
191  convertFacCF2Fmpq_poly_t (FLINTB, B);
192
193  fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
194  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
195
196  fmpq_poly_clear (FLINTA);
197  fmpq_poly_clear (FLINTB);
198  return A;
199}
200
201CanonicalForm
202mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
203                 const Variable& alpha, int m)
204{
205  CanonicalForm A= F;
206  CanonicalForm B= G;
207
208  CanonicalForm denA= bCommonDen (A);
209  CanonicalForm denB= bCommonDen (B);
210
211  A *= denA;
212  B *= denB;
213
214  int degAa= degree (A, alpha);
215  int degBa= degree (B, alpha);
216  int d= degAa + 1 + degBa;
217
218  fmpz_poly_t FLINTA,FLINTB;
219  fmpz_poly_init (FLINTA);
220  fmpz_poly_init (FLINTB);
221  kronSub (FLINTA, A, d);
222  kronSub (FLINTB, B, d);
223
224  int k= d*m;
225  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
226
227  denA *= denB;
228  A= reverseSubstQa (FLINTA, d, alpha, denA);
229  fmpz_poly_clear (FLINTA);
230  fmpz_poly_clear (FLINTB);
231  return A;
232}
233
234CanonicalForm
235mulFLINTQTrunc (const CanonicalForm& F, const CanonicalForm& G, int m)
236{
237  if (F.inCoeffDomain() || G.inCoeffDomain())
238    return mod (F*G, power (Variable (1), m));
239  Variable alpha;
240  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
241    return mulFLINTQaTrunc (F, G, alpha, m);
242
243  CanonicalForm A= F;
244  CanonicalForm B= G;
245
246  CanonicalForm denA= bCommonDen (A);
247  CanonicalForm denB= bCommonDen (B);
248
249  A *= denA;
250  B *= denB;
251  fmpz_poly_t FLINTA,FLINTB;
252  convertFacCF2Fmpz_poly_t (FLINTA, A);
253  convertFacCF2Fmpz_poly_t (FLINTB, B);
254  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, m);
255  denA *= denB;
256  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
257  A /= denA;
258  fmpz_poly_clear (FLINTA);
259  fmpz_poly_clear (FLINTB);
260
261  return A;
262}
263
264CanonicalForm uniReverse (const CanonicalForm& F, int d)
265{
266  if (d == 0)
267    return F;
268  if (F.inCoeffDomain())
269    return F*power (Variable (1),d);
270  Variable x= Variable (1);
271  CanonicalForm result= 0;
272  CFIterator i= F;
273  while (d - i.exp() < 0)
274    i++;
275
276  for (; i.hasTerms() && (d - i.exp() >= 0); i++)
277    result += i.coeff()*power (x, d - i.exp());
278  return result;
279}
280
281CanonicalForm
282newtonInverse (const CanonicalForm& F, const int n)
283{
284  int l= ilog2(n);
285
286  CanonicalForm g= F [0];
287
288  ASSERT (!g.isZero(), "expected a unit");
289
290  if (!g.isOne())
291    g = 1/g;
292  Variable x= Variable (1);
293  CanonicalForm result;
294  int exp= 0;
295  if (n & 1)
296  {
297    result= g;
298    exp= 1;
299  }
300  CanonicalForm h;
301
302  for (int i= 1; i <= l; i++)
303  {
304    h= mulNTL (g, mod (F, power (x, (1 << i))));
305    h= mod (h, power (x, (1 << i)) - 1);
306    h= div (h, power (x, (1 << (i - 1))));
307    g -= power (x, (1 << (i - 1)))*
308         mulFLINTQTrunc (g, h, 1 << (i-1));
309
310    if (n & (1 << i))
311    {
312      if (exp)
313      {
314        h= mulNTL (result, mod (F, power (x, exp + (1 << i))));
315        h= mod (h, power (x, exp + (1 << i)) - 1);
316        h= div (h, power (x, exp));
317        result -= power(x, exp)*mulFLINTQTrunc (g, h, 1 << i);
318        exp += (1 << i);
319      }
320      else
321      {
322        exp= (1 << i);
323        result= g;
324      }
325    }
326  }
327
328  return result;
329}
330
331void
332newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
333              CanonicalForm& R)
334{
335  CanonicalForm A= F;
336  CanonicalForm B= G;
337  Variable x= Variable (1);
338  int degA= degree (A, x);
339  int degB= degree (B, x);
340  int m= degA - degB;
341
342  if (m < 0)
343  {
344    R= A;
345    Q= 0;
346    return;
347  }
348
349  if (degB <= 1)
350    divrem (A, B, Q, R);
351  else
352  {
353    R= uniReverse (A, degA);
354
355    CanonicalForm revB= uniReverse (B, degB);
356    CanonicalForm buf= revB;
357    revB= newtonInverse (revB, m + 1);
358    Q= mulFLINTQTrunc (R, revB, m + 1);
359    Q= uniReverse (Q, m);
360
361    R= A - mulNTL (Q, B);
362  }
363}
364
365void
366newtonDiv (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q)
367{
368  CanonicalForm A= F;
369  CanonicalForm B= G;
370  Variable x= Variable (1);
371  int degA= degree (A, x);
372  int degB= degree (B, x);
373  int m= degA - degB;
374
375  if (m < 0)
376  {
377    Q= 0;
378    return;
379  }
380
381  if (degB <= 1)
382    Q= div (A, B);
383  else
384  {
385    CanonicalForm R= uniReverse (A, degA);
386
387    CanonicalForm revB= uniReverse (B, degB);
388    revB= newtonInverse (revB, m + 1);
389    Q= mulFLINTQTrunc (R, revB, m + 1);
390    Q= uniReverse (Q, m);
391  }
392}
393
394#endif
395
396CanonicalForm
397mulNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
398{
399  if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
400  {
401    Variable alpha;
402#ifdef HAVE_FLINT
403    if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
404        (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
405    {
406      CanonicalForm result= mulFLINTQa (F, G, alpha);
407      return result;
408    }
409    else if (!F.inCoeffDomain() && !G.inCoeffDomain())
410    {
411      if (b.getp() != 0)
412      {
413        fmpz_t FLINTpk;
414        fmpz_init_set_ui (FLINTpk, b.getp());
415        fmpz_pow_ui (FLINTpk, FLINTpk, b.getk());
416        fmpz_mod_poly_t FLINTF, FLINTG;
417        convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
418        convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
419        fmpz_mod_poly_mul (FLINTF, FLINTF, FLINTG);
420        CanonicalForm result= convertFmpz_mod_poly_t2FacCF (FLINTF, F.mvar(), b);
421        fmpz_mod_poly_clear (FLINTG);
422        fmpz_mod_poly_clear (FLINTF);
423        return result;
424      }
425      return mulFLINTQ (F, G);
426    }
427#endif
428    if (b.getp() != 0)
429      return b (F*G);
430    return F*G;
431  }
432  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
433  ASSERT (F.level() == G.level(), "expected polys of same level");
434  if (CFFactory::gettype() == GaloisFieldDomain)
435    return F*G;
436  zz_p::init (getCharacteristic());
437  Variable alpha;
438  CanonicalForm result;
439  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
440  {
441    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
442    zz_pE::init (NTLMipo);
443    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
444    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
445    mul (NTLF, NTLF, NTLG);
446    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
447  }
448  else
449  {
450#ifdef HAVE_FLINT
451    nmod_poly_t FLINTF, FLINTG;
452    convertFacCF2nmod_poly_t (FLINTF, F);
453    convertFacCF2nmod_poly_t (FLINTG, G);
454    nmod_poly_mul (FLINTF, FLINTF, FLINTG);
455    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
456    nmod_poly_clear (FLINTF);
457    nmod_poly_clear (FLINTG);
458#else
459    zz_pX NTLF= convertFacCF2NTLzzpX (F);
460    zz_pX NTLG= convertFacCF2NTLzzpX (G);
461    mul (NTLF, NTLF, NTLG);
462    result= convertNTLzzpX2CF(NTLF, F.mvar());
463#endif
464  }
465  return result;
466}
467
468CanonicalForm
469modNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
470{
471  if (F.inCoeffDomain() && G.isUnivariate())
472  {
473    if (b.getp() != 0)
474      return b(F);
475    return F;
476  }
477  else if (F.inCoeffDomain() && G.inCoeffDomain())
478  {
479    if (b.getp() != 0)
480      return b(F%G);
481    return mod (F, G);
482  }
483  else if (F.isUnivariate() && G.inCoeffDomain())
484  {
485    if (b.getp() != 0)
486      return b(F%G);
487    return mod (F,G);
488  }
489
490  if (getCharacteristic() == 0)
491  {
492#ifdef HAVE_FLINT
493    Variable alpha;
494    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
495    {
496      if (b.getp() != 0)
497      {
498        fmpz_t FLINTpk;
499        fmpz_init_set_ui (FLINTpk, b.getp());
500        fmpz_pow_ui (FLINTpk, FLINTpk, b.getk());
501        fmpz_mod_poly_t FLINTF, FLINTG;
502        convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
503        convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
504        fmpz_mod_poly_rem (FLINTF, FLINTF, FLINTG);
505        CanonicalForm result= convertFmpz_mod_poly_t2FacCF (FLINTF,F.mvar(),b);
506        fmpz_mod_poly_clear (FLINTG);
507        fmpz_mod_poly_clear (FLINTF);
508        return result;
509      }
510      return modFLINTQ (F, G);
511    }
512    else
513    {
514      CanonicalForm Q, R;
515      newtonDivrem (F, G, Q, R);
516      return R;
517    }
518#else
519    if (b.getp() != 0)
520    {
521      ZZ NTLpk= power_ZZ (b.getp(), b.getk());
522      ZZ_p::init (NTLpk);
523      ZZX ZZf= convertFacCF2NTLZZX (F);
524      ZZX ZZg= convertFacCF2NTLZZX (G);
525      ZZ_pX NTLf= to_ZZ_pX (ZZf);
526      ZZ_pX NTLg= to_ZZ_pX (ZZg);
527      rem (NTLf, NTLf, NTLg);
528      return b (convertNTLZZX2CF (to_ZZX (NTLf), F.mvar()));
529    }
530    return mod (F, G);
531#endif
532  }
533
534  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
535  ASSERT (F.level() == G.level(), "expected polys of same level");
536  if (CFFactory::gettype() == GaloisFieldDomain)
537    return mod (F, G);
538  zz_p::init (getCharacteristic());
539  Variable alpha;
540  CanonicalForm result;
541  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
542  {
543    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
544    zz_pE::init (NTLMipo);
545    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
546    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
547    rem (NTLF, NTLF, NTLG);
548    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
549  }
550  else
551  {
552#ifdef HAVE_FLINT
553    nmod_poly_t FLINTF, FLINTG;
554    convertFacCF2nmod_poly_t (FLINTF, F);
555    convertFacCF2nmod_poly_t (FLINTG, G);
556    nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
557    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
558    nmod_poly_clear (FLINTF);
559    nmod_poly_clear (FLINTG);
560#else
561    zz_pX NTLF= convertFacCF2NTLzzpX (F);
562    zz_pX NTLG= convertFacCF2NTLzzpX (G);
563    rem (NTLF, NTLF, NTLG);
564    result= convertNTLzzpX2CF(NTLF, F.mvar());
565#endif
566  }
567  return result;
568}
569
570CanonicalForm
571divNTL (const CanonicalForm& F, const CanonicalForm& G, const modpk& b)
572{
573  if (F.inCoeffDomain() && G.isUnivariate())
574  {
575    if (b.getp() != 0)
576      return b(F);
577    return F;
578  }
579  else if (F.inCoeffDomain() && G.inCoeffDomain())
580  {
581    if (b.getp() != 0)
582      return b(div (F,G));
583    return div (F, G);
584  }
585  else if (F.isUnivariate() && G.inCoeffDomain())
586  {
587    if (b.getp() != 0)
588      return b(div (F,G));
589    return div (F, G);
590  }
591
592  if (getCharacteristic() == 0)
593  {
594#ifdef HAVE_FLINT
595    Variable alpha;
596    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
597    {
598      if (b.getp() != 0)
599      {
600        fmpz_t FLINTpk;
601        fmpz_init_set_ui (FLINTpk, b.getp());
602        fmpz_pow_ui (FLINTpk, FLINTpk, b.getk());
603        fmpz_mod_poly_t FLINTF, FLINTG;
604        convertFacCF2Fmpz_mod_poly_t (FLINTF, F, FLINTpk);
605        convertFacCF2Fmpz_mod_poly_t (FLINTG, G, FLINTpk);
606        fmpz_mod_poly_divrem (FLINTF, FLINTG, FLINTF, FLINTG);
607        CanonicalForm result= convertFmpz_mod_poly_t2FacCF (FLINTF,F.mvar(),b);
608        fmpz_mod_poly_clear (FLINTG);
609        fmpz_mod_poly_clear (FLINTF);
610        return result;
611      }
612      return divFLINTQ (F,G);
613    }
614    else
615    {
616      CanonicalForm Q;
617      newtonDiv (F, G, Q);
618      return Q;
619    }
620#else
621    if (b.getp() != 0)
622    {
623      ZZ NTLpk= power_ZZ (b.getp(), b.getk());
624      ZZ_p::init (NTLpk);
625      ZZX ZZf= convertFacCF2NTLZZX (F);
626      ZZX ZZg= convertFacCF2NTLZZX (G);
627      ZZ_pX NTLf= to_ZZ_pX (ZZf);
628      ZZ_pX NTLg= to_ZZ_pX (ZZg);
629      div (NTLf, NTLf, NTLg);
630      return b (convertNTLZZX2CF (to_ZZX (NTLf), F.mvar()));
631    }
632    return div (F, G);
633#endif
634  }
635
636  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
637  ASSERT (F.level() == G.level(), "expected polys of same level");
638  if (CFFactory::gettype() == GaloisFieldDomain)
639    return div (F, G);
640  zz_p::init (getCharacteristic());
641  Variable alpha;
642  CanonicalForm result;
643  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
644  {
645    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
646    zz_pE::init (NTLMipo);
647    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
648    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
649    div (NTLF, NTLF, NTLG);
650    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
651  }
652  else
653  {
654#ifdef HAVE_FLINT
655    nmod_poly_t FLINTF, FLINTG;
656    convertFacCF2nmod_poly_t (FLINTF, F);
657    convertFacCF2nmod_poly_t (FLINTG, G);
658    nmod_poly_div (FLINTF, FLINTF, FLINTG);
659    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
660    nmod_poly_clear (FLINTF);
661    nmod_poly_clear (FLINTG);
662#else
663    zz_pX NTLF= convertFacCF2NTLzzpX (F);
664    zz_pX NTLG= convertFacCF2NTLzzpX (G);
665    div (NTLF, NTLF, NTLG);
666    result= convertNTLzzpX2CF(NTLF, F.mvar());
667#endif
668  }
669  return result;
670}
671
672// end univariate polys
673//*************************
674// bivariate polys
675
676#ifdef HAVE_FLINT
677void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
678{
679  int degAy= degree (A);
680  nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
681
682  nmod_poly_t buf;
683
684  int j, k, bufRepLength;
685  for (CFIterator i= A; i.hasTerms(); i++)
686  {
687    convertFacCF2nmod_poly_t (buf, i.coeff());
688
689    k= i.exp()*d;
690    bufRepLength= (int) nmod_poly_length (buf);
691    for (j= 0; j < bufRepLength; j++)
692      nmod_poly_set_coeff_ui (result, j + k, nmod_poly_get_coeff_ui (buf, j));
693    nmod_poly_clear (buf);
694  }
695  _nmod_poly_normalise (result);
696}
697
698void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
699{
700  int degAy= degree (A);
701  fmpq_poly_init2 (result, d1*(degAy + 1));
702
703  fmpq_poly_t buf;
704  fmpq_t coeff;
705
706  int k, l, bufRepLength;
707  CFIterator j;
708  for (CFIterator i= A; i.hasTerms(); i++)
709  {
710    if (i.coeff().inCoeffDomain())
711    {
712      k= d1*i.exp();
713      convertFacCF2Fmpq_poly_t (buf, i.coeff());
714      bufRepLength= (int) fmpq_poly_length(buf);
715      for (l= 0; l < bufRepLength; l++)
716      {
717        fmpq_poly_get_coeff_fmpq (coeff, buf, l);
718        fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
719      }
720      fmpq_poly_clear (buf);
721    }
722    else
723    {
724      for (j= i.coeff(); j.hasTerms(); j++)
725      {
726        k= d1*i.exp();
727        k += d2*j.exp();
728        convertFacCF2Fmpq_poly_t (buf, j.coeff());
729        bufRepLength= (int) fmpq_poly_length(buf);
730        for (l= 0; l < bufRepLength; l++)
731        {
732          fmpq_poly_get_coeff_fmpq (coeff, buf, l);
733          fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
734        }
735        fmpq_poly_clear (buf);
736      }
737    }
738  }
739  fmpq_clear (coeff);
740  _fmpq_poly_normalise (result);
741}
742
743void
744kronSubReciproFp (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
745                  int d)
746{
747  int degAy= degree (A);
748  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
749  nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
750  nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
751
752  nmod_poly_t buf;
753
754  int k, kk, j, bufRepLength;
755  for (CFIterator i= A; i.hasTerms(); i++)
756  {
757    convertFacCF2nmod_poly_t (buf, i.coeff());
758
759    k= i.exp()*d;
760    kk= (degAy - i.exp())*d;
761    bufRepLength= (int) nmod_poly_length (buf);
762    for (j= 0; j < bufRepLength; j++)
763    {
764      nmod_poly_set_coeff_ui (subA1, j + k,
765                              n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
766                                        nmod_poly_get_coeff_ui (buf, j),
767                                        getCharacteristic()
768                                       )
769                             );
770      nmod_poly_set_coeff_ui (subA2, j + kk,
771                              n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
772                                        nmod_poly_get_coeff_ui (buf, j),
773                                        getCharacteristic()
774                                       )
775                             );
776    }
777    nmod_poly_clear (buf);
778  }
779  _nmod_poly_normalise (subA1);
780  _nmod_poly_normalise (subA2);
781}
782
783void
784kronSubReciproQ (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
785                 int d)
786{
787  int degAy= degree (A);
788  fmpz_poly_init2 (subA1, d*(degAy + 2));
789  fmpz_poly_init2 (subA2, d*(degAy + 2));
790
791  fmpz_poly_t buf;
792  fmpz_t coeff1, coeff2;
793
794  int k, kk, j, bufRepLength;
795  for (CFIterator i= A; i.hasTerms(); i++)
796  {
797    convertFacCF2Fmpz_poly_t (buf, i.coeff());
798
799    k= i.exp()*d;
800    kk= (degAy - i.exp())*d;
801    bufRepLength= (int) fmpz_poly_length (buf);
802    for (j= 0; j < bufRepLength; j++)
803    {
804      fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k);
805      fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
806      fmpz_add (coeff1, coeff1, coeff2);
807      fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
808      fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
809      fmpz_add (coeff1, coeff1, coeff2);
810      fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
811    }
812    fmpz_poly_clear (buf);
813  }
814  fmpz_clear (coeff1);
815  fmpz_clear (coeff2);
816  _fmpz_poly_normalise (subA1);
817  _fmpz_poly_normalise (subA2);
818}
819
820CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
821{
822  Variable y= Variable (2);
823  Variable x= Variable (1);
824
825  fmpz_poly_t f;
826  fmpz_poly_init (f);
827  fmpz_poly_set (f, F);
828
829  fmpz_poly_t buf;
830  CanonicalForm result= 0;
831  int i= 0;
832  int degf= fmpz_poly_degree(f);
833  int k= 0;
834  int degfSubK, repLength, j;
835  fmpz_t coeff;
836  while (degf >= k)
837  {
838    degfSubK= degf - k;
839    if (degfSubK >= d)
840      repLength= d;
841    else
842      repLength= degfSubK + 1;
843
844    fmpz_poly_init2 (buf, repLength);
845    fmpz_init (coeff);
846    for (j= 0; j < repLength; j++)
847    {
848      fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
849      fmpz_poly_set_coeff_fmpz (buf, j, coeff);
850    }
851    _fmpz_poly_normalise (buf);
852
853    result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
854    i++;
855    k= d*i;
856    fmpz_poly_clear (buf);
857    fmpz_clear (coeff);
858  }
859  fmpz_poly_clear (f);
860
861  return result;
862}
863
864CanonicalForm
865reverseSubstQa (const fmpq_poly_t F, int d1, int d2, const Variable& alpha,
866                const fmpq_poly_t mipo)
867{
868  Variable y= Variable (2);
869  Variable x= Variable (1);
870
871  fmpq_poly_t f;
872  fmpq_poly_init (f);
873  fmpq_poly_set (f, F);
874
875  fmpq_poly_t buf;
876  CanonicalForm result= 0, result2;
877  int i= 0;
878  int degf= fmpq_poly_degree(f);
879  int k= 0;
880  int degfSubK;
881  int repLength;
882  fmpq_t coeff;
883  while (degf >= k)
884  {
885    degfSubK= degf - k;
886    if (degfSubK >= d1)
887      repLength= d1;
888    else
889      repLength= degfSubK + 1;
890
891    fmpq_init (coeff);
892    int j= 0;
893    int l;
894    result2= 0;
895    while (j*d2 < repLength)
896    {
897      fmpq_poly_init2 (buf, d2);
898      for (l= 0; l < d2; l++)
899      {
900        fmpq_poly_get_coeff_fmpq (coeff, f, k + j*d2 + l);
901        fmpq_poly_set_coeff_fmpq (buf, l, coeff);
902      }
903      _fmpq_poly_normalise (buf);
904      fmpq_poly_rem (buf, buf, mipo);
905      result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
906      j++;
907      fmpq_poly_clear (buf);
908    }
909    if (repLength - j*d2 != 0 && j*d2 - repLength < d2)
910    {
911      j--;
912      repLength -= j*d2;
913      fmpq_poly_init2 (buf, repLength);
914      j++;
915      for (l= 0; l < repLength; l++)
916      {
917        fmpq_poly_get_coeff_fmpq (coeff, f, k + j*d2 + l);
918        fmpq_poly_set_coeff_fmpq (buf, l, coeff);
919      }
920      _fmpq_poly_normalise (buf);
921      fmpq_poly_rem (buf, buf, mipo);
922      result2 += convertFmpq_poly_t2FacCF (buf, alpha)*power (x, j);
923      fmpq_poly_clear (buf);
924    }
925    fmpq_clear (coeff);
926
927    result += result2*power (y, i);
928    i++;
929    k= d1*i;
930  }
931
932  fmpq_poly_clear (f);
933  return result;
934}
935
936CanonicalForm
937reverseSubstReciproFp (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
938{
939  Variable y= Variable (2);
940  Variable x= Variable (1);
941
942  nmod_poly_t f, g;
943  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
944  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
945  nmod_poly_init_preinv (g, getCharacteristic(), ninv);
946  nmod_poly_set (f, F);
947  nmod_poly_set (g, G);
948  int degf= nmod_poly_degree(f);
949  int degg= nmod_poly_degree(g);
950
951
952  nmod_poly_t buf1,buf2, buf3;
953
954  if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
955    nmod_poly_fit_length (f,(long)d*(k+1));
956
957  CanonicalForm result= 0;
958  int i= 0;
959  int lf= 0;
960  int lg= d*k;
961  int degfSubLf= degf;
962  int deggSubLg= degg-lg;
963  int repLengthBuf2, repLengthBuf1, ind, tmp;
964  while (degf >= lf || lg >= 0)
965  {
966    if (degfSubLf >= d)
967      repLengthBuf1= d;
968    else if (degfSubLf < 0)
969      repLengthBuf1= 0;
970    else
971      repLengthBuf1= degfSubLf + 1;
972    nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
973
974    for (ind= 0; ind < repLengthBuf1; ind++)
975      nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
976    _nmod_poly_normalise (buf1);
977
978    repLengthBuf1= nmod_poly_length (buf1);
979
980    if (deggSubLg >= d - 1)
981      repLengthBuf2= d - 1;
982    else if (deggSubLg < 0)
983      repLengthBuf2= 0;
984    else
985      repLengthBuf2= deggSubLg + 1;
986
987    nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
988    for (ind= 0; ind < repLengthBuf2; ind++)
989      nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
990
991    _nmod_poly_normalise (buf2);
992    repLengthBuf2= nmod_poly_length (buf2);
993
994    nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
995    for (ind= 0; ind < repLengthBuf1; ind++)
996      nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
997    for (ind= repLengthBuf1; ind < d; ind++)
998      nmod_poly_set_coeff_ui (buf3, ind, 0);
999    for (ind= 0; ind < repLengthBuf2; ind++)
1000      nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
1001    _nmod_poly_normalise (buf3);
1002
1003    result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
1004    i++;
1005
1006
1007    lf= i*d;
1008    degfSubLf= degf - lf;
1009
1010    lg= d*(k-i);
1011    deggSubLg= degg - lg;
1012
1013    if (lg >= 0 && deggSubLg > 0)
1014    {
1015      if (repLengthBuf2 > degfSubLf + 1)
1016        degfSubLf= repLengthBuf2 - 1;
1017      tmp= tmin (repLengthBuf1, deggSubLg + 1);
1018      for (ind= 0; ind < tmp; ind++)
1019        nmod_poly_set_coeff_ui (g, ind + lg,
1020                                n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
1021                                          nmod_poly_get_coeff_ui (buf1, ind),
1022                                          getCharacteristic()
1023                                         )
1024                               );
1025    }
1026    if (lg < 0)
1027    {
1028      nmod_poly_clear (buf1);
1029      nmod_poly_clear (buf2);
1030      nmod_poly_clear (buf3);
1031      break;
1032    }
1033    if (degfSubLf >= 0)
1034    {
1035      for (ind= 0; ind < repLengthBuf2; ind++)
1036        nmod_poly_set_coeff_ui (f, ind + lf,
1037                                n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
1038                                          nmod_poly_get_coeff_ui (buf2, ind),
1039                                          getCharacteristic()
1040                                         )
1041                               );
1042    }
1043    nmod_poly_clear (buf1);
1044    nmod_poly_clear (buf2);
1045    nmod_poly_clear (buf3);
1046  }
1047
1048  nmod_poly_clear (f);
1049  nmod_poly_clear (g);
1050
1051  return result;
1052}
1053
1054CanonicalForm
1055reverseSubstReciproQ (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
1056{
1057  Variable y= Variable (2);
1058  Variable x= Variable (1);
1059
1060  fmpz_poly_t f, g;
1061  fmpz_poly_init (f);
1062  fmpz_poly_init (g);
1063  fmpz_poly_set (f, F);
1064  fmpz_poly_set (g, G);
1065  int degf= fmpz_poly_degree(f);
1066  int degg= fmpz_poly_degree(g);
1067
1068
1069  fmpz_poly_t buf1,buf2, buf3;
1070
1071  if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
1072    fmpz_poly_fit_length (f,(long)d*(k+1));
1073
1074  CanonicalForm result= 0;
1075  int i= 0;
1076  int lf= 0;
1077  int lg= d*k;
1078  int degfSubLf= degf;
1079  int deggSubLg= degg-lg;
1080  int repLengthBuf2, repLengthBuf1, ind, tmp;
1081  fmpz_t tmp1, tmp2;
1082  while (degf >= lf || lg >= 0)
1083  {
1084    if (degfSubLf >= d)
1085      repLengthBuf1= d;
1086    else if (degfSubLf < 0)
1087      repLengthBuf1= 0;
1088    else
1089      repLengthBuf1= degfSubLf + 1;
1090
1091    fmpz_poly_init2 (buf1, repLengthBuf1);
1092
1093    for (ind= 0; ind < repLengthBuf1; ind++)
1094    {
1095      fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1096      fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
1097    }
1098    _fmpz_poly_normalise (buf1);
1099
1100    repLengthBuf1= fmpz_poly_length (buf1);
1101
1102    if (deggSubLg >= d - 1)
1103      repLengthBuf2= d - 1;
1104    else if (deggSubLg < 0)
1105      repLengthBuf2= 0;
1106    else
1107      repLengthBuf2= deggSubLg + 1;
1108
1109    fmpz_poly_init2 (buf2, repLengthBuf2);
1110
1111    for (ind= 0; ind < repLengthBuf2; ind++)
1112    {
1113      fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1114      fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
1115    }
1116
1117    _fmpz_poly_normalise (buf2);
1118    repLengthBuf2= fmpz_poly_length (buf2);
1119
1120    fmpz_poly_init2 (buf3, repLengthBuf2 + d);
1121    for (ind= 0; ind < repLengthBuf1; ind++)
1122    {
1123      fmpz_poly_get_coeff_fmpz (tmp1, buf1, ind);
1124      fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
1125    }
1126    for (ind= repLengthBuf1; ind < d; ind++)
1127      fmpz_poly_set_coeff_ui (buf3, ind, 0);
1128    for (ind= 0; ind < repLengthBuf2; ind++)
1129    {
1130      fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
1131      fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
1132    }
1133    _fmpz_poly_normalise (buf3);
1134
1135    result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
1136    i++;
1137
1138
1139    lf= i*d;
1140    degfSubLf= degf - lf;
1141
1142    lg= d*(k-i);
1143    deggSubLg= degg - lg;
1144
1145    if (lg >= 0 && deggSubLg > 0)
1146    {
1147      if (repLengthBuf2 > degfSubLf + 1)
1148        degfSubLf= repLengthBuf2 - 1;
1149      tmp= tmin (repLengthBuf1, deggSubLg + 1);
1150      for (ind= 0; ind < tmp; ind++)
1151      {
1152        fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1153        fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
1154        fmpz_sub (tmp1, tmp1, tmp2);
1155        fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
1156      }
1157    }
1158    if (lg < 0)
1159    {
1160      fmpz_poly_clear (buf1);
1161      fmpz_poly_clear (buf2);
1162      fmpz_poly_clear (buf3);
1163      break;
1164    }
1165    if (degfSubLf >= 0)
1166    {
1167      for (ind= 0; ind < repLengthBuf2; ind++)
1168      {
1169        fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1170        fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
1171        fmpz_sub (tmp1, tmp1, tmp2);
1172        fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
1173      }
1174    }
1175    fmpz_poly_clear (buf1);
1176    fmpz_poly_clear (buf2);
1177    fmpz_poly_clear (buf3);
1178  }
1179
1180  fmpz_poly_clear (f);
1181  fmpz_poly_clear (g);
1182  fmpz_clear (tmp1);
1183  fmpz_clear (tmp2);
1184
1185  return result;
1186}
1187
1188CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
1189{
1190  Variable y= Variable (2);
1191  Variable x= Variable (1);
1192
1193  nmod_poly_t f;
1194  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1195  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
1196  nmod_poly_set (f, F);
1197
1198  nmod_poly_t buf;
1199  CanonicalForm result= 0;
1200  int i= 0;
1201  int degf= nmod_poly_degree(f);
1202  int k= 0;
1203  int degfSubK, repLength, j;
1204  while (degf >= k)
1205  {
1206    degfSubK= degf - k;
1207    if (degfSubK >= d)
1208      repLength= d;
1209    else
1210      repLength= degfSubK + 1;
1211
1212    nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
1213    for (j= 0; j < repLength; j++)
1214      nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
1215    _nmod_poly_normalise (buf);
1216
1217    result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
1218    i++;
1219    k= d*i;
1220    nmod_poly_clear (buf);
1221  }
1222  nmod_poly_clear (f);
1223
1224  return result;
1225}
1226
1227CanonicalForm
1228mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
1229                    CanonicalForm& M)
1230{
1231  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
1232  d1 /= 2;
1233  d1 += 1;
1234
1235  nmod_poly_t F1, F2;
1236  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1237  nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
1238  nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
1239  kronSubReciproFp (F1, F2, F, d1);
1240
1241  nmod_poly_t G1, G2;
1242  nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
1243  nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
1244  kronSubReciproFp (G1, G2, G, d1);
1245
1246  int k= d1*degree (M);
1247  nmod_poly_mullow (F1, F1, G1, (long) k);
1248
1249  int degtailF= degree (tailcoeff (F), 1);;
1250  int degtailG= degree (tailcoeff (G), 1);
1251  int taildegF= taildegree (F);
1252  int taildegG= taildegree (G);
1253
1254  int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
1255         + d1*(2+taildegF + taildegG);
1256  nmod_poly_mulhigh (F2, F2, G2, b);
1257  nmod_poly_shift_right (F2, F2, b);
1258  int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
1259
1260
1261  CanonicalForm result= reverseSubstReciproFp (F1, F2, d1, d2);
1262
1263  nmod_poly_clear (F1);
1264  nmod_poly_clear (F2);
1265  nmod_poly_clear (G1);
1266  nmod_poly_clear (G2);
1267  return result;
1268}
1269
1270CanonicalForm
1271mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
1272                CanonicalForm& M)
1273{
1274  CanonicalForm A= F;
1275  CanonicalForm B= G;
1276
1277  int degAx= degree (A, 1);
1278  int degAy= degree (A, 2);
1279  int degBx= degree (B, 1);
1280  int degBy= degree (B, 2);
1281  int d1= degAx + 1 + degBx;
1282  int d2= tmax (degAy, degBy);
1283
1284  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
1285    return mulMod2FLINTFpReci (A, B, M);
1286
1287  nmod_poly_t FLINTA, FLINTB;
1288  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1289  nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
1290  nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
1291  kronSubFp (FLINTA, A, d1);
1292  kronSubFp (FLINTB, B, d1);
1293
1294  int k= d1*degree (M);
1295  nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
1296
1297  A= reverseSubstFp (FLINTA, d1);
1298
1299  nmod_poly_clear (FLINTA);
1300  nmod_poly_clear (FLINTB);
1301  return A;
1302}
1303
1304CanonicalForm
1305mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
1306                    CanonicalForm& M)
1307{
1308  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
1309  d1 /= 2;
1310  d1 += 1;
1311
1312  fmpz_poly_t F1, F2;
1313  fmpz_poly_init (F1);
1314  fmpz_poly_init (F2);
1315  kronSubReciproQ (F1, F2, F, d1);
1316
1317  fmpz_poly_t G1, G2;
1318  fmpz_poly_init (G1);
1319  fmpz_poly_init (G2);
1320  kronSubReciproQ (G1, G2, G, d1);
1321
1322  int k= d1*degree (M);
1323  fmpz_poly_mullow (F1, F1, G1, (long) k);
1324
1325  int degtailF= degree (tailcoeff (F), 1);;
1326  int degtailG= degree (tailcoeff (G), 1);
1327  int taildegF= taildegree (F);
1328  int taildegG= taildegree (G);
1329
1330  int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
1331         + d1*(2+taildegF + taildegG);
1332  fmpz_poly_mulhigh_n (F2, F2, G2, b);
1333  fmpz_poly_shift_right (F2, F2, b);
1334  int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
1335
1336  CanonicalForm result= reverseSubstReciproQ (F1, F2, d1, d2);
1337
1338  fmpz_poly_clear (F1);
1339  fmpz_poly_clear (F2);
1340  fmpz_poly_clear (G1);
1341  fmpz_poly_clear (G2);
1342  return result;
1343}
1344
1345CanonicalForm
1346mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
1347                CanonicalForm& M)
1348{
1349  CanonicalForm A= F;
1350  CanonicalForm B= G;
1351
1352  int degAx= degree (A, 1);
1353  int degBx= degree (B, 1);
1354  int d1= degAx + 1 + degBx;
1355
1356  CanonicalForm f= bCommonDen (F);
1357  CanonicalForm g= bCommonDen (G);
1358  A *= f;
1359  B *= g;
1360
1361  fmpz_poly_t FLINTA, FLINTB;
1362  fmpz_poly_init (FLINTA);
1363  fmpz_poly_init (FLINTB);
1364  kronSub (FLINTA, A, d1);
1365  kronSub (FLINTB, B, d1);
1366  int k= d1*degree (M);
1367
1368  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
1369  A= reverseSubstQ (FLINTA, d1);
1370  fmpz_poly_clear (FLINTA);
1371  fmpz_poly_clear (FLINTB);
1372  return A/(f*g);
1373}
1374
1375CanonicalForm
1376mulMod2FLINTQa (const CanonicalForm& F, const CanonicalForm& G,
1377                const CanonicalForm& M)
1378{
1379  Variable a;
1380  if (!hasFirstAlgVar (F,a) && !hasFirstAlgVar (G, a))
1381    return mulMod2FLINTQ (F, G, M);
1382  CanonicalForm A= F;
1383
1384  int degFx= degree (F, 1);
1385  int degFa= degree (F, a);
1386  int degGx= degree (G, 1);
1387  int degGa= degree (G, a);
1388
1389  int d2= degFa+degGa+1;
1390  int d1= degFx + 1 + degGx;
1391  d1 *= d2;
1392
1393  fmpq_poly_t FLINTF, FLINTG;
1394  kronSubQa (FLINTF, F, d1, d2);
1395  kronSubQa (FLINTG, G, d1, d2);
1396
1397  fmpq_poly_mullow (FLINTF, FLINTF, FLINTG, d1*degree (M));
1398
1399  fmpq_poly_t mipo;
1400  convertFacCF2Fmpq_poly_t (mipo, getMipo (a));
1401  CanonicalForm result= reverseSubstQa (FLINTF, d1, d2, a, mipo);
1402  fmpq_poly_clear (FLINTF);
1403  fmpq_poly_clear (FLINTG);
1404  return result;
1405}
1406
1407#endif
1408
1409zz_pX kronSubFp (const CanonicalForm& A, int d)
1410{
1411  int degAy= degree (A);
1412  zz_pX result;
1413  result.rep.SetLength (d*(degAy + 1));
1414
1415  zz_p *resultp;
1416  resultp= result.rep.elts();
1417  zz_pX buf;
1418  zz_p *bufp;
1419  int j, k, bufRepLength;
1420
1421  for (CFIterator i= A; i.hasTerms(); i++)
1422  {
1423    if (i.coeff().inCoeffDomain())
1424      buf= convertFacCF2NTLzzpX (i.coeff());
1425    else
1426      buf= convertFacCF2NTLzzpX (i.coeff());
1427
1428    k= i.exp()*d;
1429    bufp= buf.rep.elts();
1430    bufRepLength= (int) buf.rep.length();
1431    for (j= 0; j < bufRepLength; j++)
1432      resultp [j + k]= bufp [j];
1433  }
1434  result.normalize();
1435
1436  return result;
1437}
1438
1439zz_pEX kronSubFq (const CanonicalForm& A, int d, const Variable& alpha)
1440{
1441  int degAy= degree (A);
1442  zz_pEX result;
1443  result.rep.SetLength (d*(degAy + 1));
1444
1445  Variable v;
1446  zz_pE *resultp;
1447  resultp= result.rep.elts();
1448  zz_pEX buf1;
1449  zz_pE *buf1p;
1450  zz_pX buf2;
1451  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1452  int j, k, buf1RepLength;
1453
1454  for (CFIterator i= A; i.hasTerms(); i++)
1455  {
1456    if (i.coeff().inCoeffDomain())
1457    {
1458      buf2= convertFacCF2NTLzzpX (i.coeff());
1459      buf1= to_zz_pEX (to_zz_pE (buf2));
1460    }
1461    else
1462      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
1463
1464    k= i.exp()*d;
1465    buf1p= buf1.rep.elts();
1466    buf1RepLength= (int) buf1.rep.length();
1467    for (j= 0; j < buf1RepLength; j++)
1468      resultp [j + k]= buf1p [j];
1469  }
1470  result.normalize();
1471
1472  return result;
1473}
1474
1475void
1476kronSubReciproFq (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
1477                  const Variable& alpha)
1478{
1479  int degAy= degree (A);
1480  subA1.rep.SetLength ((long) d*(degAy + 2));
1481  subA2.rep.SetLength ((long) d*(degAy + 2));
1482
1483  Variable v;
1484  zz_pE *subA1p;
1485  zz_pE *subA2p;
1486  subA1p= subA1.rep.elts();
1487  subA2p= subA2.rep.elts();
1488  zz_pEX buf;
1489  zz_pE *bufp;
1490  zz_pX buf2;
1491  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1492  int j, k, kk, bufRepLength;
1493
1494  for (CFIterator i= A; i.hasTerms(); i++)
1495  {
1496    if (i.coeff().inCoeffDomain())
1497    {
1498      buf2= convertFacCF2NTLzzpX (i.coeff());
1499      buf= to_zz_pEX (to_zz_pE (buf2));
1500    }
1501    else
1502      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
1503
1504    k= i.exp()*d;
1505    kk= (degAy - i.exp())*d;
1506    bufp= buf.rep.elts();
1507    bufRepLength= (int) buf.rep.length();
1508    for (j= 0; j < bufRepLength; j++)
1509    {
1510      subA1p [j + k] += bufp [j];
1511      subA2p [j + kk] += bufp [j];
1512    }
1513  }
1514  subA1.normalize();
1515  subA2.normalize();
1516}
1517
1518void
1519kronSubReciproFp (zz_pX& subA1, zz_pX& subA2, const CanonicalForm& A, int d)
1520{
1521  int degAy= degree (A);
1522  subA1.rep.SetLength ((long) d*(degAy + 2));
1523  subA2.rep.SetLength ((long) d*(degAy + 2));
1524
1525  zz_p *subA1p;
1526  zz_p *subA2p;
1527  subA1p= subA1.rep.elts();
1528  subA2p= subA2.rep.elts();
1529  zz_pX buf;
1530  zz_p *bufp;
1531  int j, k, kk, bufRepLength;
1532
1533  for (CFIterator i= A; i.hasTerms(); i++)
1534  {
1535    buf= convertFacCF2NTLzzpX (i.coeff());
1536
1537    k= i.exp()*d;
1538    kk= (degAy - i.exp())*d;
1539    bufp= buf.rep.elts();
1540    bufRepLength= (int) buf.rep.length();
1541    for (j= 0; j < bufRepLength; j++)
1542    {
1543      subA1p [j + k] += bufp [j];
1544      subA2p [j + kk] += bufp [j];
1545    }
1546  }
1547  subA1.normalize();
1548  subA2.normalize();
1549}
1550
1551CanonicalForm
1552reverseSubstReciproFq (const zz_pEX& F, const zz_pEX& G, int d, int k,
1553                       const Variable& alpha)
1554{
1555  Variable y= Variable (2);
1556  Variable x= Variable (1);
1557
1558  zz_pEX f= F;
1559  zz_pEX g= G;
1560  int degf= deg(f);
1561  int degg= deg(g);
1562
1563  zz_pEX buf1;
1564  zz_pEX buf2;
1565  zz_pEX buf3;
1566
1567  zz_pE *buf1p;
1568  zz_pE *buf2p;
1569  zz_pE *buf3p;
1570  if (f.rep.length() < (long) d*(k+1)) //zero padding
1571    f.rep.SetLength ((long)d*(k+1));
1572
1573  zz_pE *gp= g.rep.elts();
1574  zz_pE *fp= f.rep.elts();
1575  CanonicalForm result= 0;
1576  int i= 0;
1577  int lf= 0;
1578  int lg= d*k;
1579  int degfSubLf= degf;
1580  int deggSubLg= degg-lg;
1581  int repLengthBuf2, repLengthBuf1, ind, tmp;
1582  zz_pE zzpEZero= zz_pE();
1583
1584  while (degf >= lf || lg >= 0)
1585  {
1586    if (degfSubLf >= d)
1587      repLengthBuf1= d;
1588    else if (degfSubLf < 0)
1589      repLengthBuf1= 0;
1590    else
1591      repLengthBuf1= degfSubLf + 1;
1592    buf1.rep.SetLength((long) repLengthBuf1);
1593
1594    buf1p= buf1.rep.elts();
1595    for (ind= 0; ind < repLengthBuf1; ind++)
1596      buf1p [ind]= fp [ind + lf];
1597    buf1.normalize();
1598
1599    repLengthBuf1= buf1.rep.length();
1600
1601    if (deggSubLg >= d - 1)
1602      repLengthBuf2= d - 1;
1603    else if (deggSubLg < 0)
1604      repLengthBuf2= 0;
1605    else
1606      repLengthBuf2= deggSubLg + 1;
1607
1608    buf2.rep.SetLength ((long) repLengthBuf2);
1609    buf2p= buf2.rep.elts();
1610    for (ind= 0; ind < repLengthBuf2; ind++)
1611      buf2p [ind]= gp [ind + lg];
1612    buf2.normalize();
1613
1614    repLengthBuf2= buf2.rep.length();
1615
1616    buf3.rep.SetLength((long) repLengthBuf2 + d);
1617    buf3p= buf3.rep.elts();
1618    buf2p= buf2.rep.elts();
1619    buf1p= buf1.rep.elts();
1620    for (ind= 0; ind < repLengthBuf1; ind++)
1621      buf3p [ind]= buf1p [ind];
1622    for (ind= repLengthBuf1; ind < d; ind++)
1623      buf3p [ind]= zzpEZero;
1624    for (ind= 0; ind < repLengthBuf2; ind++)
1625      buf3p [ind + d]= buf2p [ind];
1626    buf3.normalize();
1627
1628    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
1629    i++;
1630
1631
1632    lf= i*d;
1633    degfSubLf= degf - lf;
1634
1635    lg= d*(k-i);
1636    deggSubLg= degg - lg;
1637
1638    buf1p= buf1.rep.elts();
1639
1640    if (lg >= 0 && deggSubLg > 0)
1641    {
1642      if (repLengthBuf2 > degfSubLf + 1)
1643        degfSubLf= repLengthBuf2 - 1;
1644      tmp= tmin (repLengthBuf1, deggSubLg + 1);
1645      for (ind= 0; ind < tmp; ind++)
1646        gp [ind + lg] -= buf1p [ind];
1647    }
1648
1649    if (lg < 0)
1650      break;
1651
1652    buf2p= buf2.rep.elts();
1653    if (degfSubLf >= 0)
1654    {
1655      for (ind= 0; ind < repLengthBuf2; ind++)
1656        fp [ind + lf] -= buf2p [ind];
1657    }
1658  }
1659
1660  return result;
1661}
1662
1663CanonicalForm
1664reverseSubstReciproFp (const zz_pX& F, const zz_pX& G, int d, int k)
1665{
1666  Variable y= Variable (2);
1667  Variable x= Variable (1);
1668
1669  zz_pX f= F;
1670  zz_pX g= G;
1671  int degf= deg(f);
1672  int degg= deg(g);
1673
1674  zz_pX buf1;
1675  zz_pX buf2;
1676  zz_pX buf3;
1677
1678  zz_p *buf1p;
1679  zz_p *buf2p;
1680  zz_p *buf3p;
1681
1682  if (f.rep.length() < (long) d*(k+1)) //zero padding
1683    f.rep.SetLength ((long)d*(k+1));
1684
1685  zz_p *gp= g.rep.elts();
1686  zz_p *fp= f.rep.elts();
1687  CanonicalForm result= 0;
1688  int i= 0;
1689  int lf= 0;
1690  int lg= d*k;
1691  int degfSubLf= degf;
1692  int deggSubLg= degg-lg;
1693  int repLengthBuf2, repLengthBuf1, ind, tmp;
1694  zz_p zzpZero= zz_p();
1695  while (degf >= lf || lg >= 0)
1696  {
1697    if (degfSubLf >= d)
1698      repLengthBuf1= d;
1699    else if (degfSubLf < 0)
1700      repLengthBuf1= 0;
1701    else
1702      repLengthBuf1= degfSubLf + 1;
1703    buf1.rep.SetLength((long) repLengthBuf1);
1704
1705    buf1p= buf1.rep.elts();
1706    for (ind= 0; ind < repLengthBuf1; ind++)
1707      buf1p [ind]= fp [ind + lf];
1708    buf1.normalize();
1709
1710    repLengthBuf1= buf1.rep.length();
1711
1712    if (deggSubLg >= d - 1)
1713      repLengthBuf2= d - 1;
1714    else if (deggSubLg < 0)
1715      repLengthBuf2= 0;
1716    else
1717      repLengthBuf2= deggSubLg + 1;
1718
1719    buf2.rep.SetLength ((long) repLengthBuf2);
1720    buf2p= buf2.rep.elts();
1721    for (ind= 0; ind < repLengthBuf2; ind++)
1722      buf2p [ind]= gp [ind + lg];
1723
1724    buf2.normalize();
1725
1726    repLengthBuf2= buf2.rep.length();
1727
1728
1729    buf3.rep.SetLength((long) repLengthBuf2 + d);
1730    buf3p= buf3.rep.elts();
1731    buf2p= buf2.rep.elts();
1732    buf1p= buf1.rep.elts();
1733    for (ind= 0; ind < repLengthBuf1; ind++)
1734      buf3p [ind]= buf1p [ind];
1735    for (ind= repLengthBuf1; ind < d; ind++)
1736      buf3p [ind]= zzpZero;
1737    for (ind= 0; ind < repLengthBuf2; ind++)
1738      buf3p [ind + d]= buf2p [ind];
1739    buf3.normalize();
1740
1741    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
1742    i++;
1743
1744
1745    lf= i*d;
1746    degfSubLf= degf - lf;
1747
1748    lg= d*(k-i);
1749    deggSubLg= degg - lg;
1750
1751    buf1p= buf1.rep.elts();
1752
1753    if (lg >= 0 && deggSubLg > 0)
1754    {
1755      if (repLengthBuf2 > degfSubLf + 1)
1756        degfSubLf= repLengthBuf2 - 1;
1757      tmp= tmin (repLengthBuf1, deggSubLg + 1);
1758      for (ind= 0; ind < tmp; ind++)
1759        gp [ind + lg] -= buf1p [ind];
1760    }
1761    if (lg < 0)
1762      break;
1763
1764    buf2p= buf2.rep.elts();
1765    if (degfSubLf >= 0)
1766    {
1767      for (ind= 0; ind < repLengthBuf2; ind++)
1768        fp [ind + lf] -= buf2p [ind];
1769    }
1770  }
1771
1772  return result;
1773}
1774
1775CanonicalForm reverseSubstFq (const zz_pEX& F, int d, const Variable& alpha)
1776{
1777  Variable y= Variable (2);
1778  Variable x= Variable (1);
1779
1780  zz_pEX f= F;
1781  zz_pE *fp= f.rep.elts();
1782
1783  zz_pEX buf;
1784  zz_pE *bufp;
1785  CanonicalForm result= 0;
1786  int i= 0;
1787  int degf= deg(f);
1788  int k= 0;
1789  int degfSubK, repLength, j;
1790  while (degf >= k)
1791  {
1792    degfSubK= degf - k;
1793    if (degfSubK >= d)
1794      repLength= d;
1795    else
1796      repLength= degfSubK + 1;
1797
1798    buf.rep.SetLength ((long) repLength);
1799    bufp= buf.rep.elts();
1800    for (j= 0; j < repLength; j++)
1801      bufp [j]= fp [j + k];
1802    buf.normalize();
1803
1804    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
1805    i++;
1806    k= d*i;
1807  }
1808
1809  return result;
1810}
1811
1812CanonicalForm reverseSubstFp (const zz_pX& F, int d)
1813{
1814  Variable y= Variable (2);
1815  Variable x= Variable (1);
1816
1817  zz_pX f= F;
1818  zz_p *fp= f.rep.elts();
1819
1820  zz_pX buf;
1821  zz_p *bufp;
1822  CanonicalForm result= 0;
1823  int i= 0;
1824  int degf= deg(f);
1825  int k= 0;
1826  int degfSubK, repLength, j;
1827  while (degf >= k)
1828  {
1829    degfSubK= degf - k;
1830    if (degfSubK >= d)
1831      repLength= d;
1832    else
1833      repLength= degfSubK + 1;
1834
1835    buf.rep.SetLength ((long) repLength);
1836    bufp= buf.rep.elts();
1837    for (j= 0; j < repLength; j++)
1838      bufp [j]= fp [j + k];
1839    buf.normalize();
1840
1841    result += convertNTLzzpX2CF (buf, x)*power (y, i);
1842    i++;
1843    k= d*i;
1844  }
1845
1846  return result;
1847}
1848
1849// assumes input to be reduced mod M and to be an element of Fq not Fp
1850CanonicalForm
1851mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
1852                  CanonicalForm& M)
1853{
1854  int d1= degree (F, 1) + degree (G, 1) + 1;
1855  d1 /= 2;
1856  d1 += 1;
1857
1858  zz_pX F1, F2;
1859  kronSubReciproFp (F1, F2, F, d1);
1860  zz_pX G1, G2;
1861  kronSubReciproFp (G1, G2, G, d1);
1862
1863  int k= d1*degree (M);
1864  MulTrunc (F1, F1, G1, (long) k);
1865
1866  int degtailF= degree (tailcoeff (F), 1);
1867  int degtailG= degree (tailcoeff (G), 1);
1868  int taildegF= taildegree (F);
1869  int taildegG= taildegree (G);
1870  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
1871
1872  reverse (F2, F2);
1873  reverse (G2, G2);
1874  MulTrunc (F2, F2, G2, b + 1);
1875  reverse (F2, F2, b);
1876
1877  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
1878  return reverseSubstReciproFp (F1, F2, d1, d2);
1879}
1880
1881//Kronecker substitution
1882CanonicalForm
1883mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
1884              CanonicalForm& M)
1885{
1886  CanonicalForm A= F;
1887  CanonicalForm B= G;
1888
1889  int degAx= degree (A, 1);
1890  int degAy= degree (A, 2);
1891  int degBx= degree (B, 1);
1892  int degBy= degree (B, 2);
1893  int d1= degAx + 1 + degBx;
1894  int d2= tmax (degAy, degBy);
1895
1896  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
1897    return mulMod2NTLFpReci (A, B, M);
1898
1899  zz_pX NTLA= kronSubFp (A, d1);
1900  zz_pX NTLB= kronSubFp (B, d1);
1901
1902  int k= d1*degree (M);
1903  MulTrunc (NTLA, NTLA, NTLB, (long) k);
1904
1905  A= reverseSubstFp (NTLA, d1);
1906
1907  return A;
1908}
1909
1910// assumes input to be reduced mod M and to be an element of Fq not Fp
1911CanonicalForm
1912mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
1913                  CanonicalForm& M, const Variable& alpha)
1914{
1915  int d1= degree (F, 1) + degree (G, 1) + 1;
1916  d1 /= 2;
1917  d1 += 1;
1918
1919  zz_pEX F1, F2;
1920  kronSubReciproFq (F1, F2, F, d1, alpha);
1921  zz_pEX G1, G2;
1922  kronSubReciproFq (G1, G2, G, d1, alpha);
1923
1924  int k= d1*degree (M);
1925  MulTrunc (F1, F1, G1, (long) k);
1926
1927  int degtailF= degree (tailcoeff (F), 1);
1928  int degtailG= degree (tailcoeff (G), 1);
1929  int taildegF= taildegree (F);
1930  int taildegG= taildegree (G);
1931  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
1932
1933  reverse (F2, F2);
1934  reverse (G2, G2);
1935  MulTrunc (F2, F2, G2, b + 1);
1936  reverse (F2, F2, b);
1937
1938  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
1939  return reverseSubstReciproFq (F1, F2, d1, d2, alpha);
1940}
1941
1942#ifdef HAVE_FLINT
1943CanonicalForm
1944mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
1945                CanonicalForm& M);
1946#endif
1947
1948CanonicalForm
1949mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
1950              CanonicalForm& M)
1951{
1952  Variable alpha;
1953  CanonicalForm A= F;
1954  CanonicalForm B= G;
1955
1956  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
1957  {
1958    int degAx= degree (A, 1);
1959    int degAy= degree (A, 2);
1960    int degBx= degree (B, 1);
1961    int degBy= degree (B, 2);
1962    int d1= degAx + degBx + 1;
1963    int d2= tmax (degAy, degBy);
1964    zz_p::init (getCharacteristic());
1965    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1966    zz_pE::init (NTLMipo);
1967
1968    int degMipo= degree (getMipo (alpha));
1969    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
1970        (2*degAy > degree (M)))
1971      return mulMod2NTLFqReci (A, B, M, alpha);
1972
1973    zz_pEX NTLA= kronSubFq (A, d1, alpha);
1974    zz_pEX NTLB= kronSubFq (B, d1, alpha);
1975
1976    int k= d1*degree (M);
1977
1978    MulTrunc (NTLA, NTLA, NTLB, (long) k);
1979
1980    A= reverseSubstFq (NTLA, d1, alpha);
1981
1982    return A;
1983  }
1984  else
1985#ifdef HAVE_FLINT
1986    return mulMod2FLINTFp (A, B, M);
1987#else
1988    return mulMod2NTLFp (A, B, M);
1989#endif
1990}
1991
1992CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
1993                       const CanonicalForm& M)
1994{
1995  if (A.isZero() || B.isZero())
1996    return 0;
1997
1998  ASSERT (M.isUnivariate(), "M must be univariate");
1999
2000  CanonicalForm F= mod (A, M);
2001  CanonicalForm G= mod (B, M);
2002  if (F.inCoeffDomain() || G.inCoeffDomain())
2003    return F*G;
2004  Variable y= M.mvar();
2005  int degF= degree (F, y);
2006  int degG= degree (G, y);
2007
2008  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
2009      (F.level() == G.level()))
2010  {
2011    CanonicalForm result= mulNTL (F, G);
2012    return mod (result, M);
2013  }
2014  else if (degF <= 1 && degG <= 1)
2015  {
2016    CanonicalForm result= F*G;
2017    return mod (result, M);
2018  }
2019
2020  int sizeF= size (F);
2021  int sizeG= size (G);
2022
2023  int fallBackToNaive= 50;
2024  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
2025    return mod (F*G, M);
2026
2027#ifdef HAVE_FLINT
2028  if (getCharacteristic() == 0)
2029    return mulMod2FLINTQa (F, G, M);
2030#endif
2031
2032  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
2033      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
2034    return mulMod2NTLFq (F, G, M);
2035
2036  int m= (int) ceil (degree (M)/2.0);
2037  if (degF >= m || degG >= m)
2038  {
2039    CanonicalForm MLo= power (y, m);
2040    CanonicalForm MHi= power (y, degree (M) - m);
2041    CanonicalForm F0= mod (F, MLo);
2042    CanonicalForm F1= div (F, MLo);
2043    CanonicalForm G0= mod (G, MLo);
2044    CanonicalForm G1= div (G, MLo);
2045    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
2046    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
2047    CanonicalForm F0G0= mulMod2 (F0, G0, M);
2048    return F0G0 + MLo*(F0G1 + F1G0);
2049  }
2050  else
2051  {
2052    m= (int) ceil (tmax (degF, degG)/2.0);
2053    CanonicalForm yToM= power (y, m);
2054    CanonicalForm F0= mod (F, yToM);
2055    CanonicalForm F1= div (F, yToM);
2056    CanonicalForm G0= mod (G, yToM);
2057    CanonicalForm G1= div (G, yToM);
2058    CanonicalForm H00= mulMod2 (F0, G0, M);
2059    CanonicalForm H11= mulMod2 (F1, G1, M);
2060    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
2061    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
2062  }
2063  DEBOUTLN (cerr, "fatal end in mulMod2");
2064}
2065
2066// end bivariate polys
2067//**********************
2068// multivariate polys
2069
2070CanonicalForm mod (const CanonicalForm& F, const CFList& M)
2071{
2072  CanonicalForm A= F;
2073  for (CFListIterator i= M; i.hasItem(); i++)
2074    A= mod (A, i.getItem());
2075  return A;
2076}
2077
2078CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
2079                      const CFList& MOD)
2080{
2081  if (A.isZero() || B.isZero())
2082    return 0;
2083
2084  if (MOD.length() == 1)
2085    return mulMod2 (A, B, MOD.getLast());
2086
2087  CanonicalForm M= MOD.getLast();
2088  CanonicalForm F= mod (A, M);
2089  CanonicalForm G= mod (B, M);
2090  if (F.inCoeffDomain() || G.inCoeffDomain())
2091    return F*G;
2092  Variable y= M.mvar();
2093  int degF= degree (F, y);
2094  int degG= degree (G, y);
2095
2096  if ((degF <= 1 && F.level() <= M.level()) &&
2097      (degG <= 1 && G.level() <= M.level()))
2098  {
2099    CFList buf= MOD;
2100    buf.removeLast();
2101    if (degF == 1 && degG == 1)
2102    {
2103      CanonicalForm F0= mod (F, y);
2104      CanonicalForm F1= div (F, y);
2105      CanonicalForm G0= mod (G, y);
2106      CanonicalForm G1= div (G, y);
2107      if (degree (M) > 2)
2108      {
2109        CanonicalForm H00= mulMod (F0, G0, buf);
2110        CanonicalForm H11= mulMod (F1, G1, buf);
2111        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
2112        return H11*y*y + (H01 - H00 - H11)*y + H00;
2113      }
2114      else //here degree (M) == 2
2115      {
2116        buf.append (y);
2117        CanonicalForm F0G1= mulMod (F0, G1, buf);
2118        CanonicalForm F1G0= mulMod (F1, G0, buf);
2119        CanonicalForm F0G0= mulMod (F0, G0, MOD);
2120        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
2121        return result;
2122      }
2123    }
2124    else if (degF == 1 && degG == 0)
2125      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
2126    else if (degF == 0 && degG == 1)
2127      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
2128    else
2129      return mulMod (F, G, buf);
2130  }
2131  int m= (int) ceil (degree (M)/2.0);
2132  if (degF >= m || degG >= m)
2133  {
2134    CanonicalForm MLo= power (y, m);
2135    CanonicalForm MHi= power (y, degree (M) - m);
2136    CanonicalForm F0= mod (F, MLo);
2137    CanonicalForm F1= div (F, MLo);
2138    CanonicalForm G0= mod (G, MLo);
2139    CanonicalForm G1= div (G, MLo);
2140    CFList buf= MOD;
2141    buf.removeLast();
2142    buf.append (MHi);
2143    CanonicalForm F0G1= mulMod (F0, G1, buf);
2144    CanonicalForm F1G0= mulMod (F1, G0, buf);
2145    CanonicalForm F0G0= mulMod (F0, G0, MOD);
2146    return F0G0 + MLo*(F0G1 + F1G0);
2147  }
2148  else
2149  {
2150    m= (int) ceil (tmax (degF, degG)/2.0);
2151    CanonicalForm yToM= power (y, m);
2152    CanonicalForm F0= mod (F, yToM);
2153    CanonicalForm F1= div (F, yToM);
2154    CanonicalForm G0= mod (G, yToM);
2155    CanonicalForm G1= div (G, yToM);
2156    CanonicalForm H00= mulMod (F0, G0, MOD);
2157    CanonicalForm H11= mulMod (F1, G1, MOD);
2158    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
2159    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
2160  }
2161  DEBOUTLN (cerr, "fatal end in mulMod");
2162}
2163
2164CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
2165{
2166  if (L.isEmpty())
2167    return 1;
2168  int l= L.length();
2169  if (l == 1)
2170    return mod (L.getFirst(), M);
2171  else if (l == 2) {
2172    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
2173    return result;
2174  }
2175  else
2176  {
2177    l /= 2;
2178    CFList tmp1, tmp2;
2179    CFListIterator i= L;
2180    CanonicalForm buf1, buf2;
2181    for (int j= 1; j <= l; j++, i++)
2182      tmp1.append (i.getItem());
2183    tmp2= Difference (L, tmp1);
2184    buf1= prodMod (tmp1, M);
2185    buf2= prodMod (tmp2, M);
2186    CanonicalForm result= mulMod2 (buf1, buf2, M);
2187    return result;
2188  }
2189}
2190
2191CanonicalForm prodMod (const CFList& L, const CFList& M)
2192{
2193  if (L.isEmpty())
2194    return 1;
2195  else if (L.length() == 1)
2196    return L.getFirst();
2197  else if (L.length() == 2)
2198    return mulMod (L.getFirst(), L.getLast(), M);
2199  else
2200  {
2201    int l= L.length()/2;
2202    CFListIterator i= L;
2203    CFList tmp1, tmp2;
2204    CanonicalForm buf1, buf2;
2205    for (int j= 1; j <= l; j++, i++)
2206      tmp1.append (i.getItem());
2207    tmp2= Difference (L, tmp1);
2208    buf1= prodMod (tmp1, M);
2209    buf2= prodMod (tmp2, M);
2210    return mulMod (buf1, buf2, M);
2211  }
2212}
2213
2214// end multivariate polys
2215//***************************
2216// division
2217
2218CanonicalForm reverse (const CanonicalForm& F, int d)
2219{
2220  if (d == 0)
2221    return F;
2222  CanonicalForm A= F;
2223  Variable y= Variable (2);
2224  Variable x= Variable (1);
2225  if (degree (A, x) > 0)
2226  {
2227    A= swapvar (A, x, y);
2228    CanonicalForm result= 0;
2229    CFIterator i= A;
2230    while (d - i.exp() < 0)
2231      i++;
2232
2233    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
2234      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
2235    return result;
2236  }
2237  else
2238    return A*power (x, d);
2239}
2240
2241CanonicalForm
2242newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
2243{
2244  int l= ilog2(n);
2245
2246  CanonicalForm g= mod (F, M)[0] [0];
2247
2248  ASSERT (!g.isZero(), "expected a unit");
2249
2250  Variable alpha;
2251
2252  if (!g.isOne())
2253    g = 1/g;
2254  Variable x= Variable (1);
2255  CanonicalForm result;
2256  int exp= 0;
2257  if (n & 1)
2258  {
2259    result= g;
2260    exp= 1;
2261  }
2262  CanonicalForm h;
2263
2264  for (int i= 1; i <= l; i++)
2265  {
2266    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
2267    h= mod (h, power (x, (1 << i)) - 1);
2268    h= div (h, power (x, (1 << (i - 1))));
2269    h= mod (h, M);
2270    g -= power (x, (1 << (i - 1)))*
2271         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
2272
2273    if (n & (1 << i))
2274    {
2275      if (exp)
2276      {
2277        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
2278        h= mod (h, power (x, exp + (1 << i)) - 1);
2279        h= div (h, power (x, exp));
2280        h= mod (h, M);
2281        result -= power(x, exp)*mod (mulMod2 (g, h, M),
2282                                       power (x, (1 << i)));
2283        exp += (1 << i);
2284      }
2285      else
2286      {
2287        exp= (1 << i);
2288        result= g;
2289      }
2290    }
2291  }
2292
2293  return result;
2294}
2295
2296CanonicalForm
2297newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
2298           M)
2299{
2300  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
2301  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
2302
2303  CanonicalForm A= mod (F, M);
2304  CanonicalForm B= mod (G, M);
2305
2306  Variable x= Variable (1);
2307  int degA= degree (A, x);
2308  int degB= degree (B, x);
2309  int m= degA - degB;
2310  if (m < 0)
2311    return 0;
2312
2313  Variable v;
2314  CanonicalForm Q;
2315  if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
2316  {
2317    CanonicalForm R;
2318    divrem2 (A, B, Q, R, M);
2319  }
2320  else
2321  {
2322    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
2323    {
2324      CanonicalForm R= reverse (A, degA);
2325      CanonicalForm revB= reverse (B, degB);
2326      revB= newtonInverse (revB, m + 1, M);
2327      Q= mulMod2 (R, revB, M);
2328      Q= mod (Q, power (x, m + 1));
2329      Q= reverse (Q, m);
2330    }
2331    else
2332    {
2333      zz_pX mipo= convertFacCF2NTLzzpX (M);
2334      Variable y= Variable (2);
2335      zz_pEX NTLA, NTLB;
2336      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
2337      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
2338      div (NTLA, NTLA, NTLB);
2339      Q= convertNTLzz_pEX2CF (NTLA, x, y);
2340    }
2341  }
2342
2343  return Q;
2344}
2345
2346void
2347newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2348              CanonicalForm& R, const CanonicalForm& M)
2349{
2350  CanonicalForm A= mod (F, M);
2351  CanonicalForm B= mod (G, M);
2352  Variable x= Variable (1);
2353  int degA= degree (A, x);
2354  int degB= degree (B, x);
2355  int m= degA - degB;
2356
2357  if (m < 0)
2358  {
2359    R= A;
2360    Q= 0;
2361    return;
2362  }
2363
2364  Variable v;
2365  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
2366  {
2367     divrem2 (A, B, Q, R, M);
2368  }
2369  else
2370  {
2371    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
2372    {
2373      R= reverse (A, degA);
2374
2375      CanonicalForm revB= reverse (B, degB);
2376      revB= newtonInverse (revB, m + 1, M);
2377      Q= mulMod2 (R, revB, M);
2378
2379      Q= mod (Q, power (x, m + 1));
2380      Q= reverse (Q, m);
2381
2382      R= A - mulMod2 (Q, B, M);
2383    }
2384    else
2385    {
2386      zz_pX mipo= convertFacCF2NTLzzpX (M);
2387      Variable y= Variable (2);
2388      zz_pEX NTLA, NTLB;
2389      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
2390      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
2391      zz_pEX NTLQ, NTLR;
2392      DivRem (NTLQ, NTLR, NTLA, NTLB);
2393      Q= convertNTLzz_pEX2CF (NTLQ, x, y);
2394      R= convertNTLzz_pEX2CF (NTLR, x, y);
2395    }
2396  }
2397}
2398
2399static inline
2400CFList split (const CanonicalForm& F, const int m, const Variable& x)
2401{
2402  CanonicalForm A= F;
2403  CanonicalForm buf= 0;
2404  bool swap= false;
2405  if (degree (A, x) <= 0)
2406    return CFList(A);
2407  else if (x.level() != A.level())
2408  {
2409    swap= true;
2410    A= swapvar (A, x, A.mvar());
2411  }
2412
2413  int j= (int) floor ((double) degree (A)/ m);
2414  CFList result;
2415  CFIterator i= A;
2416  for (; j >= 0; j--)
2417  {
2418    while (i.hasTerms() && i.exp() - j*m >= 0)
2419    {
2420      if (swap)
2421        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
2422      else
2423        buf += i.coeff()*power (x, i.exp() - j*m);
2424      i++;
2425    }
2426    if (swap)
2427      result.append (swapvar (buf, x, F.mvar()));
2428    else
2429      result.append (buf);
2430    buf= 0;
2431  }
2432  return result;
2433}
2434
2435static inline
2436void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2437               CanonicalForm& R, const CFList& M);
2438
2439static inline
2440void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2441               CanonicalForm& R, const CFList& M)
2442{
2443  CanonicalForm A= mod (F, M);
2444  CanonicalForm B= mod (G, M);
2445  Variable x= Variable (1);
2446  int degB= degree (B, x);
2447  int degA= degree (A, x);
2448  if (degA < degB)
2449  {
2450    Q= 0;
2451    R= A;
2452    return;
2453  }
2454  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
2455  if (degB < 1)
2456  {
2457    divrem (A, B, Q, R);
2458    Q= mod (Q, M);
2459    R= mod (R, M);
2460    return;
2461  }
2462
2463  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
2464  CFList splitA= split (A, m, x);
2465  if (splitA.length() == 3)
2466    splitA.insert (0);
2467  if (splitA.length() == 2)
2468  {
2469    splitA.insert (0);
2470    splitA.insert (0);
2471  }
2472  if (splitA.length() == 1)
2473  {
2474    splitA.insert (0);
2475    splitA.insert (0);
2476    splitA.insert (0);
2477  }
2478
2479  CanonicalForm xToM= power (x, m);
2480
2481  CFListIterator i= splitA;
2482  CanonicalForm H= i.getItem();
2483  i++;
2484  H *= xToM;
2485  H += i.getItem();
2486  i++;
2487  H *= xToM;
2488  H += i.getItem();
2489  i++;
2490
2491  divrem32 (H, B, Q, R, M);
2492
2493  CFList splitR= split (R, m, x);
2494  if (splitR.length() == 1)
2495    splitR.insert (0);
2496
2497  H= splitR.getFirst();
2498  H *= xToM;
2499  H += splitR.getLast();
2500  H *= xToM;
2501  H += i.getItem();
2502
2503  CanonicalForm bufQ;
2504  divrem32 (H, B, bufQ, R, M);
2505
2506  Q *= xToM;
2507  Q += bufQ;
2508  return;
2509}
2510
2511static inline
2512void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2513               CanonicalForm& R, const CFList& M)
2514{
2515  CanonicalForm A= mod (F, M);
2516  CanonicalForm B= mod (G, M);
2517  Variable x= Variable (1);
2518  int degB= degree (B, x);
2519  int degA= degree (A, x);
2520  if (degA < degB)
2521  {
2522    Q= 0;
2523    R= A;
2524    return;
2525  }
2526  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
2527  if (degB < 1)
2528  {
2529    divrem (A, B, Q, R);
2530    Q= mod (Q, M);
2531    R= mod (R, M);
2532    return;
2533  }
2534  int m= (int) ceil ((double) (degB + 1)/ 2.0);
2535
2536  CFList splitA= split (A, m, x);
2537  CFList splitB= split (B, m, x);
2538
2539  if (splitA.length() == 2)
2540  {
2541    splitA.insert (0);
2542  }
2543  if (splitA.length() == 1)
2544  {
2545    splitA.insert (0);
2546    splitA.insert (0);
2547  }
2548  CanonicalForm xToM= power (x, m);
2549
2550  CanonicalForm H;
2551  CFListIterator i= splitA;
2552  i++;
2553
2554  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
2555  {
2556    H= splitA.getFirst()*xToM + i.getItem();
2557    divrem21 (H, splitB.getFirst(), Q, R, M);
2558  }
2559  else
2560  {
2561    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
2562       splitB.getFirst()*xToM;
2563    Q= xToM - 1;
2564  }
2565
2566  H= mulMod (Q, splitB.getLast(), M);
2567
2568  R= R*xToM + splitA.getLast() - H;
2569
2570  while (degree (R, x) >= degB)
2571  {
2572    xToM= power (x, degree (R, x) - degB);
2573    Q += LC (R, x)*xToM;
2574    R -= mulMod (LC (R, x), B, M)*xToM;
2575    Q= mod (Q, M);
2576    R= mod (R, M);
2577  }
2578
2579  return;
2580}
2581
2582void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2583              CanonicalForm& R, const CanonicalForm& M)
2584{
2585  CanonicalForm A= mod (F, M);
2586  CanonicalForm B= mod (G, M);
2587
2588  if (B.inCoeffDomain())
2589  {
2590    divrem (A, B, Q, R);
2591    return;
2592  }
2593  if (A.inCoeffDomain() && !B.inCoeffDomain())
2594  {
2595    Q= 0;
2596    R= A;
2597    return;
2598  }
2599
2600  if (B.level() < A.level())
2601  {
2602    divrem (A, B, Q, R);
2603    return;
2604  }
2605  if (A.level() > B.level())
2606  {
2607    R= A;
2608    Q= 0;
2609    return;
2610  }
2611  if (B.level() == 1 && B.isUnivariate())
2612  {
2613    divrem (A, B, Q, R);
2614    return;
2615  }
2616  if (!(B.level() == 1 && B.isUnivariate()) &&
2617      (A.level() == 1 && A.isUnivariate()))
2618  {
2619    Q= 0;
2620    R= A;
2621    return;
2622  }
2623
2624  Variable x= Variable (1);
2625  int degB= degree (B, x);
2626  if (degB > degree (A, x))
2627  {
2628    Q= 0;
2629    R= A;
2630    return;
2631  }
2632
2633  CFList splitA= split (A, degB, x);
2634
2635  CanonicalForm xToDegB= power (x, degB);
2636  CanonicalForm H, bufQ;
2637  Q= 0;
2638  CFListIterator i= splitA;
2639  H= i.getItem()*xToDegB;
2640  i++;
2641  H += i.getItem();
2642  CFList buf;
2643  while (i.hasItem())
2644  {
2645    buf= CFList (M);
2646    divrem21 (H, B, bufQ, R, buf);
2647    i++;
2648    if (i.hasItem())
2649      H= R*xToDegB + i.getItem();
2650    Q *= xToDegB;
2651    Q += bufQ;
2652  }
2653  return;
2654}
2655
2656void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2657             CanonicalForm& R, const CFList& MOD)
2658{
2659  CanonicalForm A= mod (F, MOD);
2660  CanonicalForm B= mod (G, MOD);
2661  Variable x= Variable (1);
2662  int degB= degree (B, x);
2663  if (degB > degree (A, x))
2664  {
2665    Q= 0;
2666    R= A;
2667    return;
2668  }
2669
2670  if (degB <= 0)
2671  {
2672    divrem (A, B, Q, R);
2673    Q= mod (Q, MOD);
2674    R= mod (R, MOD);
2675    return;
2676  }
2677  CFList splitA= split (A, degB, x);
2678
2679  CanonicalForm xToDegB= power (x, degB);
2680  CanonicalForm H, bufQ;
2681  Q= 0;
2682  CFListIterator i= splitA;
2683  H= i.getItem()*xToDegB;
2684  i++;
2685  H += i.getItem();
2686  while (i.hasItem())
2687  {
2688    divrem21 (H, B, bufQ, R, MOD);
2689    i++;
2690    if (i.hasItem())
2691      H= R*xToDegB + i.getItem();
2692    Q *= xToDegB;
2693    Q += bufQ;
2694  }
2695  return;
2696}
2697
2698bool
2699uniFdivides (const CanonicalForm& A, const CanonicalForm& B)
2700{
2701  int p= getCharacteristic();
2702  if (p > 0)
2703  {
2704    zz_p::init (p);
2705    Variable alpha;
2706    if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
2707    {
2708      zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2709      zz_pE::init (NTLMipo);
2710      zz_pEX NTLA= convertFacCF2NTLzz_pEX (A, NTLMipo);
2711      zz_pEX NTLB= convertFacCF2NTLzz_pEX (B, NTLMipo);
2712      return divide (NTLB, NTLA);
2713    }
2714#ifdef HAVE_FLINT
2715    nmod_poly_t FLINTA, FLINTB;
2716    convertFacCF2nmod_poly_t (FLINTA, A);
2717    convertFacCF2nmod_poly_t (FLINTB, B);
2718    nmod_poly_rem (FLINTA, FLINTB, FLINTA);
2719    bool result= nmod_poly_is_zero (FLINTA);
2720    nmod_poly_clear (FLINTA);
2721    nmod_poly_clear (FLINTB);
2722    return result;
2723#else
2724    zz_pX NTLA= convertFacCF2NTLzzpX (A);
2725    zz_pX NTLB= convertFacCF2NTLzzpX (B);
2726    return divide (NTLB, NTLA);
2727#endif
2728  }
2729#ifdef HAVE_FLINT
2730  Variable alpha;
2731  if (!hasFirstAlgVar (A, alpha) && !hasFirstAlgVar (B, alpha))
2732  {
2733    fmpq_poly_t FLINTA,FLINTB;
2734    fmpq_poly_init (FLINTA);
2735    fmpq_poly_init (FLINTB);
2736    convertFacCF2Fmpq_poly_t (FLINTA, A);
2737    convertFacCF2Fmpq_poly_t (FLINTB, B);
2738    fmpq_poly_rem (FLINTA, FLINTB, FLINTA);
2739    bool result= fmpq_poly_is_zero (FLINTA);
2740    fmpq_poly_clear (FLINTA);
2741    fmpq_poly_clear (FLINTB);
2742    return result;
2743  }
2744  else
2745    return true;
2746    //return fdivides (A, B);
2747#else
2748  return fdivides (A, B); //maybe NTL?
2749#endif
2750}
2751
2752// end division
2753
2754#endif
Note: See TracBrowser for help on using the repository browser.