source: git/factory/facMul.cc @ e785e9

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