source: git/factory/facMul.cc @ bffe62

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