source: git/factory/facMul.cc @ aba90e8

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