source: git/factory/facMul.cc @ e016ba

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