source: git/factory/facMul.cc @ 3e0c11d

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