source: git/factory/facMul.cc @ 7a1151

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