source: git/factory/facMul.cc @ 237c42

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