source: git/factory/facMul.cc @ 81d96c

spielwiese
Last change on this file since 81d96c was 81d96c, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: more reorganization
  • Property mode set to 100644
File size: 58.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)
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);  //oder fmpz_set (fmpz_poly_get_coeff_ptr (buf3, ind),fmpz_poly_get_coeff_ptr (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) && ! hasFirstAlgVar (G, alpha))
1829    return mulMod2FLINTQ (F, G, M);
1830#endif
1831
1832  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
1833      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
1834    return mulMod2NTLFq (F, G, M);
1835
1836  int m= (int) ceil (degree (M)/2.0);
1837  if (degF >= m || degG >= m)
1838  {
1839    CanonicalForm MLo= power (y, m);
1840    CanonicalForm MHi= power (y, degree (M) - m);
1841    CanonicalForm F0= mod (F, MLo);
1842    CanonicalForm F1= div (F, MLo);
1843    CanonicalForm G0= mod (G, MLo);
1844    CanonicalForm G1= div (G, MLo);
1845    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
1846    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
1847    CanonicalForm F0G0= mulMod2 (F0, G0, M);
1848    return F0G0 + MLo*(F0G1 + F1G0);
1849  }
1850  else
1851  {
1852    m= (int) ceil (tmax (degF, degG)/2.0);
1853    CanonicalForm yToM= power (y, m);
1854    CanonicalForm F0= mod (F, yToM);
1855    CanonicalForm F1= div (F, yToM);
1856    CanonicalForm G0= mod (G, yToM);
1857    CanonicalForm G1= div (G, yToM);
1858    CanonicalForm H00= mulMod2 (F0, G0, M);
1859    CanonicalForm H11= mulMod2 (F1, G1, M);
1860    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
1861    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
1862  }
1863  DEBOUTLN (cerr, "fatal end in mulMod2");
1864}
1865
1866// end bivariate polys
1867//**********************
1868// multivariate polys
1869
1870CanonicalForm mod (const CanonicalForm& F, const CFList& M)
1871{
1872  CanonicalForm A= F;
1873  for (CFListIterator i= M; i.hasItem(); i++)
1874    A= mod (A, i.getItem());
1875  return A;
1876}
1877
1878CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
1879                      const CFList& MOD)
1880{
1881  if (A.isZero() || B.isZero())
1882    return 0;
1883
1884  if (MOD.length() == 1)
1885    return mulMod2 (A, B, MOD.getLast());
1886
1887  CanonicalForm M= MOD.getLast();
1888  CanonicalForm F= mod (A, M);
1889  CanonicalForm G= mod (B, M);
1890  if (F.inCoeffDomain() || G.inCoeffDomain())
1891    return F*G;
1892  Variable y= M.mvar();
1893  int degF= degree (F, y);
1894  int degG= degree (G, y);
1895
1896  if ((degF <= 1 && F.level() <= M.level()) &&
1897      (degG <= 1 && G.level() <= M.level()))
1898  {
1899    CFList buf= MOD;
1900    buf.removeLast();
1901    if (degF == 1 && degG == 1)
1902    {
1903      CanonicalForm F0= mod (F, y);
1904      CanonicalForm F1= div (F, y);
1905      CanonicalForm G0= mod (G, y);
1906      CanonicalForm G1= div (G, y);
1907      if (degree (M) > 2)
1908      {
1909        CanonicalForm H00= mulMod (F0, G0, buf);
1910        CanonicalForm H11= mulMod (F1, G1, buf);
1911        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
1912        return H11*y*y + (H01 - H00 - H11)*y + H00;
1913      }
1914      else //here degree (M) == 2
1915      {
1916        buf.append (y);
1917        CanonicalForm F0G1= mulMod (F0, G1, buf);
1918        CanonicalForm F1G0= mulMod (F1, G0, buf);
1919        CanonicalForm F0G0= mulMod (F0, G0, MOD);
1920        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
1921        return result;
1922      }
1923    }
1924    else if (degF == 1 && degG == 0)
1925      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
1926    else if (degF == 0 && degG == 1)
1927      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
1928    else
1929      return mulMod (F, G, buf);
1930  }
1931  int m= (int) ceil (degree (M)/2.0);
1932  if (degF >= m || degG >= m)
1933  {
1934    CanonicalForm MLo= power (y, m);
1935    CanonicalForm MHi= power (y, degree (M) - m);
1936    CanonicalForm F0= mod (F, MLo);
1937    CanonicalForm F1= div (F, MLo);
1938    CanonicalForm G0= mod (G, MLo);
1939    CanonicalForm G1= div (G, MLo);
1940    CFList buf= MOD;
1941    buf.removeLast();
1942    buf.append (MHi);
1943    CanonicalForm F0G1= mulMod (F0, G1, buf);
1944    CanonicalForm F1G0= mulMod (F1, G0, buf);
1945    CanonicalForm F0G0= mulMod (F0, G0, MOD);
1946    return F0G0 + MLo*(F0G1 + F1G0);
1947  }
1948  else
1949  {
1950    m= (int) ceil (tmax (degF, degG)/2.0);
1951    CanonicalForm yToM= power (y, m);
1952    CanonicalForm F0= mod (F, yToM);
1953    CanonicalForm F1= div (F, yToM);
1954    CanonicalForm G0= mod (G, yToM);
1955    CanonicalForm G1= div (G, yToM);
1956    CanonicalForm H00= mulMod (F0, G0, MOD);
1957    CanonicalForm H11= mulMod (F1, G1, MOD);
1958    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
1959    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
1960  }
1961  DEBOUTLN (cerr, "fatal end in mulMod");
1962}
1963
1964CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
1965{
1966  if (L.isEmpty())
1967    return 1;
1968  int l= L.length();
1969  if (l == 1)
1970    return mod (L.getFirst(), M);
1971  else if (l == 2) {
1972    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
1973    return result;
1974  }
1975  else
1976  {
1977    l /= 2;
1978    CFList tmp1, tmp2;
1979    CFListIterator i= L;
1980    CanonicalForm buf1, buf2;
1981    for (int j= 1; j <= l; j++, i++)
1982      tmp1.append (i.getItem());
1983    tmp2= Difference (L, tmp1);
1984    buf1= prodMod (tmp1, M);
1985    buf2= prodMod (tmp2, M);
1986    CanonicalForm result= mulMod2 (buf1, buf2, M);
1987    return result;
1988  }
1989}
1990
1991CanonicalForm prodMod (const CFList& L, const CFList& M)
1992{
1993  if (L.isEmpty())
1994    return 1;
1995  else if (L.length() == 1)
1996    return L.getFirst();
1997  else if (L.length() == 2)
1998    return mulMod (L.getFirst(), L.getLast(), M);
1999  else
2000  {
2001    int l= L.length()/2;
2002    CFListIterator i= L;
2003    CFList tmp1, tmp2;
2004    CanonicalForm buf1, buf2;
2005    for (int j= 1; j <= l; j++, i++)
2006      tmp1.append (i.getItem());
2007    tmp2= Difference (L, tmp1);
2008    buf1= prodMod (tmp1, M);
2009    buf2= prodMod (tmp2, M);
2010    return mulMod (buf1, buf2, M);
2011  }
2012}
2013
2014// end multivariate polys
2015//***************************
2016// division
2017
2018CanonicalForm reverse (const CanonicalForm& F, int d)
2019{
2020  if (d == 0)
2021    return F;
2022  CanonicalForm A= F;
2023  Variable y= Variable (2);
2024  Variable x= Variable (1);
2025  if (degree (A, x) > 0)
2026  {
2027    A= swapvar (A, x, y);
2028    CanonicalForm result= 0;
2029    CFIterator i= A;
2030    while (d - i.exp() < 0)
2031      i++;
2032
2033    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
2034      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
2035    return result;
2036  }
2037  else
2038    return A*power (x, d);
2039}
2040
2041CanonicalForm
2042newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
2043{
2044  int l= ilog2(n);
2045
2046  CanonicalForm g= mod (F, M)[0] [0];
2047
2048  ASSERT (!g.isZero(), "expected a unit");
2049
2050  Variable alpha;
2051
2052  if (!g.isOne())
2053    g = 1/g;
2054  Variable x= Variable (1);
2055  CanonicalForm result;
2056  int exp= 0;
2057  if (n & 1)
2058  {
2059    result= g;
2060    exp= 1;
2061  }
2062  CanonicalForm h;
2063
2064  for (int i= 1; i <= l; i++)
2065  {
2066    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
2067    h= mod (h, power (x, (1 << i)) - 1);
2068    h= div (h, power (x, (1 << (i - 1))));
2069    h= mod (h, M);
2070    g -= power (x, (1 << (i - 1)))*
2071         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
2072
2073    if (n & (1 << i))
2074    {
2075      if (exp)
2076      {
2077        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
2078        h= mod (h, power (x, exp + (1 << i)) - 1);
2079        h= div (h, power (x, exp));
2080        h= mod (h, M);
2081        result -= power(x, exp)*mod (mulMod2 (g, h, M),
2082                                       power (x, (1 << i)));
2083        exp += (1 << i);
2084      }
2085      else
2086      {
2087        exp= (1 << i);
2088        result= g;
2089      }
2090    }
2091  }
2092
2093  return result;
2094}
2095
2096CanonicalForm
2097newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
2098           M)
2099{
2100  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
2101  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
2102
2103  CanonicalForm A= mod (F, M);
2104  CanonicalForm B= mod (G, M);
2105
2106  Variable x= Variable (1);
2107  int degA= degree (A, x);
2108  int degB= degree (B, x);
2109  int m= degA - degB;
2110  if (m < 0)
2111    return 0;
2112
2113  Variable v;
2114  CanonicalForm Q;
2115  if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
2116  {
2117    CanonicalForm R;
2118    divrem2 (A, B, Q, R, M);
2119  }
2120  else
2121  {
2122    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
2123    {
2124      CanonicalForm R= reverse (A, degA);
2125      CanonicalForm revB= reverse (B, degB);
2126      revB= newtonInverse (revB, m + 1, M);
2127      Q= mulMod2 (R, revB, M);
2128      Q= mod (Q, power (x, m + 1));
2129      Q= reverse (Q, m);
2130    }
2131    else
2132    {
2133      zz_pX mipo= convertFacCF2NTLzzpX (M);
2134      Variable y= Variable (2);
2135      zz_pEX NTLA, NTLB;
2136      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
2137      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
2138      div (NTLA, NTLA, NTLB);
2139      Q= convertNTLzz_pEX2CF (NTLA, x, y);
2140    }
2141  }
2142
2143  return Q;
2144}
2145
2146void
2147newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2148              CanonicalForm& R, const CanonicalForm& M)
2149{
2150  CanonicalForm A= mod (F, M);
2151  CanonicalForm B= mod (G, M);
2152  Variable x= Variable (1);
2153  int degA= degree (A, x);
2154  int degB= degree (B, x);
2155  int m= degA - degB;
2156
2157  if (m < 0)
2158  {
2159    R= A;
2160    Q= 0;
2161    return;
2162  }
2163
2164  Variable v;
2165  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
2166  {
2167     divrem2 (A, B, Q, R, M);
2168  }
2169  else
2170  {
2171    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
2172    {
2173      R= reverse (A, degA);
2174
2175      CanonicalForm revB= reverse (B, degB);
2176      revB= newtonInverse (revB, m + 1, M);
2177      Q= mulMod2 (R, revB, M);
2178
2179      Q= mod (Q, power (x, m + 1));
2180      Q= reverse (Q, m);
2181
2182      R= A - mulMod2 (Q, B, M);
2183    }
2184    else
2185    {
2186      zz_pX mipo= convertFacCF2NTLzzpX (M);
2187      Variable y= Variable (2);
2188      zz_pEX NTLA, NTLB;
2189      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
2190      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
2191      zz_pEX NTLQ, NTLR;
2192      DivRem (NTLQ, NTLR, NTLA, NTLB);
2193      Q= convertNTLzz_pEX2CF (NTLQ, x, y);
2194      R= convertNTLzz_pEX2CF (NTLR, x, y);
2195    }
2196  }
2197}
2198
2199static inline
2200CFList split (const CanonicalForm& F, const int m, const Variable& x)
2201{
2202  CanonicalForm A= F;
2203  CanonicalForm buf= 0;
2204  bool swap= false;
2205  if (degree (A, x) <= 0)
2206    return CFList(A);
2207  else if (x.level() != A.level())
2208  {
2209    swap= true;
2210    A= swapvar (A, x, A.mvar());
2211  }
2212
2213  int j= (int) floor ((double) degree (A)/ m);
2214  CFList result;
2215  CFIterator i= A;
2216  for (; j >= 0; j--)
2217  {
2218    while (i.hasTerms() && i.exp() - j*m >= 0)
2219    {
2220      if (swap)
2221        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
2222      else
2223        buf += i.coeff()*power (x, i.exp() - j*m);
2224      i++;
2225    }
2226    if (swap)
2227      result.append (swapvar (buf, x, F.mvar()));
2228    else
2229      result.append (buf);
2230    buf= 0;
2231  }
2232  return result;
2233}
2234
2235static inline
2236void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2237               CanonicalForm& R, const CFList& M);
2238
2239static inline
2240void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2241               CanonicalForm& R, const CFList& M)
2242{
2243  CanonicalForm A= mod (F, M);
2244  CanonicalForm B= mod (G, M);
2245  Variable x= Variable (1);
2246  int degB= degree (B, x);
2247  int degA= degree (A, x);
2248  if (degA < degB)
2249  {
2250    Q= 0;
2251    R= A;
2252    return;
2253  }
2254  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
2255  if (degB < 1)
2256  {
2257    divrem (A, B, Q, R);
2258    Q= mod (Q, M);
2259    R= mod (R, M);
2260    return;
2261  }
2262
2263  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
2264  CFList splitA= split (A, m, x);
2265  if (splitA.length() == 3)
2266    splitA.insert (0);
2267  if (splitA.length() == 2)
2268  {
2269    splitA.insert (0);
2270    splitA.insert (0);
2271  }
2272  if (splitA.length() == 1)
2273  {
2274    splitA.insert (0);
2275    splitA.insert (0);
2276    splitA.insert (0);
2277  }
2278
2279  CanonicalForm xToM= power (x, m);
2280
2281  CFListIterator i= splitA;
2282  CanonicalForm H= i.getItem();
2283  i++;
2284  H *= xToM;
2285  H += i.getItem();
2286  i++;
2287  H *= xToM;
2288  H += i.getItem();
2289  i++;
2290
2291  divrem32 (H, B, Q, R, M);
2292
2293  CFList splitR= split (R, m, x);
2294  if (splitR.length() == 1)
2295    splitR.insert (0);
2296
2297  H= splitR.getFirst();
2298  H *= xToM;
2299  H += splitR.getLast();
2300  H *= xToM;
2301  H += i.getItem();
2302
2303  CanonicalForm bufQ;
2304  divrem32 (H, B, bufQ, R, M);
2305
2306  Q *= xToM;
2307  Q += bufQ;
2308  return;
2309}
2310
2311static inline
2312void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2313               CanonicalForm& R, const CFList& M)
2314{
2315  CanonicalForm A= mod (F, M);
2316  CanonicalForm B= mod (G, M);
2317  Variable x= Variable (1);
2318  int degB= degree (B, x);
2319  int degA= degree (A, x);
2320  if (degA < degB)
2321  {
2322    Q= 0;
2323    R= A;
2324    return;
2325  }
2326  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
2327  if (degB < 1)
2328  {
2329    divrem (A, B, Q, R);
2330    Q= mod (Q, M);
2331    R= mod (R, M);
2332    return;
2333  }
2334  int m= (int) ceil ((double) (degB + 1)/ 2.0);
2335
2336  CFList splitA= split (A, m, x);
2337  CFList splitB= split (B, m, x);
2338
2339  if (splitA.length() == 2)
2340  {
2341    splitA.insert (0);
2342  }
2343  if (splitA.length() == 1)
2344  {
2345    splitA.insert (0);
2346    splitA.insert (0);
2347  }
2348  CanonicalForm xToM= power (x, m);
2349
2350  CanonicalForm H;
2351  CFListIterator i= splitA;
2352  i++;
2353
2354  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
2355  {
2356    H= splitA.getFirst()*xToM + i.getItem();
2357    divrem21 (H, splitB.getFirst(), Q, R, M);
2358  }
2359  else
2360  {
2361    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
2362       splitB.getFirst()*xToM;
2363    Q= xToM - 1;
2364  }
2365
2366  H= mulMod (Q, splitB.getLast(), M);
2367
2368  R= R*xToM + splitA.getLast() - H;
2369
2370  while (degree (R, x) >= degB)
2371  {
2372    xToM= power (x, degree (R, x) - degB);
2373    Q += LC (R, x)*xToM;
2374    R -= mulMod (LC (R, x), B, M)*xToM;
2375    Q= mod (Q, M);
2376    R= mod (R, M);
2377  }
2378
2379  return;
2380}
2381
2382void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2383              CanonicalForm& R, const CanonicalForm& M)
2384{
2385  CanonicalForm A= mod (F, M);
2386  CanonicalForm B= mod (G, M);
2387
2388  if (B.inCoeffDomain())
2389  {
2390    divrem (A, B, Q, R);
2391    return;
2392  }
2393  if (A.inCoeffDomain() && !B.inCoeffDomain())
2394  {
2395    Q= 0;
2396    R= A;
2397    return;
2398  }
2399
2400  if (B.level() < A.level())
2401  {
2402    divrem (A, B, Q, R);
2403    return;
2404  }
2405  if (A.level() > B.level())
2406  {
2407    R= A;
2408    Q= 0;
2409    return;
2410  }
2411  if (B.level() == 1 && B.isUnivariate())
2412  {
2413    divrem (A, B, Q, R);
2414    return;
2415  }
2416  if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
2417  {
2418    Q= 0;
2419    R= A;
2420    return;
2421  }
2422
2423  Variable x= Variable (1);
2424  int degB= degree (B, x);
2425  if (degB > degree (A, x))
2426  {
2427    Q= 0;
2428    R= A;
2429    return;
2430  }
2431
2432  CFList splitA= split (A, degB, x);
2433
2434  CanonicalForm xToDegB= power (x, degB);
2435  CanonicalForm H, bufQ;
2436  Q= 0;
2437  CFListIterator i= splitA;
2438  H= i.getItem()*xToDegB;
2439  i++;
2440  H += i.getItem();
2441  CFList buf;
2442  while (i.hasItem())
2443  {
2444    buf= CFList (M);
2445    divrem21 (H, B, bufQ, R, buf);
2446    i++;
2447    if (i.hasItem())
2448      H= R*xToDegB + i.getItem();
2449    Q *= xToDegB;
2450    Q += bufQ;
2451  }
2452  return;
2453}
2454
2455void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2456             CanonicalForm& R, const CFList& MOD)
2457{
2458  CanonicalForm A= mod (F, MOD);
2459  CanonicalForm B= mod (G, MOD);
2460  Variable x= Variable (1);
2461  int degB= degree (B, x);
2462  if (degB > degree (A, x))
2463  {
2464    Q= 0;
2465    R= A;
2466    return;
2467  }
2468
2469  if (degB <= 0)
2470  {
2471    divrem (A, B, Q, R);
2472    Q= mod (Q, MOD);
2473    R= mod (R, MOD);
2474    return;
2475  }
2476  CFList splitA= split (A, degB, x);
2477
2478  CanonicalForm xToDegB= power (x, degB);
2479  CanonicalForm H, bufQ;
2480  Q= 0;
2481  CFListIterator i= splitA;
2482  H= i.getItem()*xToDegB;
2483  i++;
2484  H += i.getItem();
2485  while (i.hasItem())
2486  {
2487    divrem21 (H, B, bufQ, R, MOD);
2488    i++;
2489    if (i.hasItem())
2490      H= R*xToDegB + i.getItem();
2491    Q *= xToDegB;
2492    Q += bufQ;
2493  }
2494  return;
2495}
2496
2497// end division
2498
2499#endif
Note: See TracBrowser for help on using the repository browser.