source: git/factory/facMul.cc @ 2773fe

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