source: git/factory/facMul.cc @ 68a3479

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