source: git/factory/facMul.cc @ f9bd3d

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