source: git/factory/facMul.cc @ 42af505

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