source: git/factory/facMul.cc @ ca058c

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