source: git/factory/facMul.cc @ ad0177

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