source: git/factory/facHensel.cc @ 0e6668

spielwiese
Last change on this file since 0e6668 was 0e6668, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: more changes in order to replace old factorization over integers
  • Property mode set to 100644
File size: 123.0 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facHensel.cc
5 *
6 * This file implements functions to lift factors via Hensel lifting and
7 * functions for modular multiplication and division with remainder.
8 *
9 * ABSTRACT: Hensel lifting is described in "Efficient Multivariate
10 * Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with
11 * remainder is described in "Fast Recursive Division" by C. Burnikel and
12 * J. Ziegler. Karatsuba multiplication is described in "Modern Computer
13 * Algebra" by J. von zur Gathen and J. Gerhard.
14 *
15 * @author Martin Lee
16 *
17 * @internal @version \$Id$
18 *
19 **/
20/*****************************************************************************/
21
22#include "config.h"
23
24#include "cf_assert.h"
25#include "debug.h"
26#include "timing.h"
27
28#include "facHensel.h"
29#include "cf_util.h"
30#include "fac_util.h"
31#include "cf_algorithm.h"
32#include "cf_primes.h"
33
34#ifdef HAVE_NTL
35#include <NTL/lzz_pEX.h>
36#include "NTLconvert.h"
37
38#ifdef HAVE_FLINT
39#include "FLINTconvert.h"
40#endif
41
42#ifdef HAVE_FLINT
43void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
44{
45  int degAy= degree (A);
46  fmpz_poly_init2 (result, d*(degAy + 1));
47  _fmpz_poly_set_length (result, d*(degAy + 1));
48  CFIterator j;
49  for (CFIterator i= A; i.hasTerms(); i++)
50  {
51    if (i.coeff().inBaseDomain())
52      convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
53    else
54      for (j=i.coeff(); j.hasTerms(); j++)
55        convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result,i.exp()*d+j.exp()), j.coeff());
56  }
57  _fmpz_poly_normalise(result);
58
59
60  /*fmpz_t coeff;
61  for (CFIterator i= A; i.hasTerms(); i++)
62  {
63    if (i.coeff().inBaseDomain())
64    {
65      convertCF2Fmpz (coeff, i.coeff());
66      fmpz_poly_set_coeff_fmpz (result, i.exp()*d, coeff);
67    }
68    else
69    {
70      for (CFIterator j= i.coeff();j.hasTerms(); j++)
71      {
72        convertCF2Fmpz (coeff, j.coeff());
73        fmpz_poly_set_coeff_fmpz (result, i.exp()*d+j.exp(), coeff);
74      }
75    }
76  }
77  fmpz_clear (coeff);
78  _fmpz_poly_normalise (result);*/
79}
80
81
82CanonicalForm reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha, const CanonicalForm& den)
83{
84  Variable x= Variable (1);
85
86  CanonicalForm result= 0;
87  int i= 0;
88  int degf= fmpz_poly_degree (F);
89  int k= 0;
90  int degfSubK;
91  int repLength;
92  CanonicalForm coeff;
93  fmpz* tmp;
94  while (degf >= k)
95  {
96    coeff= 0;
97    degfSubK= degf - k;
98    if (degfSubK >= d)
99      repLength= d;
100    else
101      repLength= degfSubK + 1;
102
103    for (int j= 0; j < repLength; j++)
104    {
105      tmp= fmpz_poly_get_coeff_ptr (F, j+k);
106      if (!fmpz_is_zero (tmp))
107      {
108        CanonicalForm ff= convertFmpz2CF (tmp)/den;
109        coeff += ff*power (alpha, j);
110      }
111    }
112    result += coeff*power (x, i);
113    i++;
114    k= d*i;
115  }
116  return result;
117}
118
119CanonicalForm mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G, const Variable& alpha)
120{
121  CanonicalForm A= F;
122  CanonicalForm B= G;
123
124  CanonicalForm denA= bCommonDen (A);
125  CanonicalForm denB= bCommonDen (B);
126
127  A *= denA;
128  B *= denB;
129  int degAa= degree (A, alpha);
130  int degBa= degree (B, alpha);
131  int d= degAa + 1 + degBa;
132
133  fmpz_poly_t FLINTA,FLINTB;
134  fmpz_poly_init (FLINTA);
135  fmpz_poly_init (FLINTB);
136  kronSub (FLINTA, A, d);
137  kronSub (FLINTB, B, d);
138
139  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
140
141  denA *= denB;
142  A= reverseSubst (FLINTA, d, alpha, denA);
143
144  fmpz_poly_clear (FLINTA);
145  fmpz_poly_clear (FLINTB);
146  return A;
147}
148
149CanonicalForm
150mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
151{
152  CanonicalForm A= F;
153  CanonicalForm B= G;
154
155  CanonicalForm denA= bCommonDen (A);
156  CanonicalForm denB= bCommonDen (B);
157
158  A *= denA;
159  B *= denB;
160  fmpz_poly_t FLINTA,FLINTB;
161  convertFacCF2Fmpz_poly_t (FLINTA, A);
162  convertFacCF2Fmpz_poly_t (FLINTB, B);
163  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
164  denA *= denB;
165  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
166  A /= denA;
167  fmpz_poly_clear (FLINTA);
168  fmpz_poly_clear (FLINTB);
169
170  return A;
171}
172
173CanonicalForm
174mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
175{
176  CanonicalForm A= F;
177  CanonicalForm B= G;
178
179  fmpq_poly_t FLINTA,FLINTB;
180  convertFacCF2Fmpq_poly_t (FLINTA, A);
181  convertFacCF2Fmpq_poly_t (FLINTB, B);
182
183  fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
184  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
185  fmpq_poly_clear (FLINTA);
186  fmpq_poly_clear (FLINTB);
187  return A;
188}
189
190CanonicalForm
191divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
192{
193  CanonicalForm A= F;
194  CanonicalForm B= G;
195
196  fmpq_poly_t FLINTA,FLINTB;
197  convertFacCF2Fmpq_poly_t (FLINTA, A);
198  convertFacCF2Fmpq_poly_t (FLINTB, B);
199
200  fmpq_poly_div (FLINTA, FLINTA, FLINTB);
201  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
202
203  fmpq_poly_clear (FLINTA);
204  fmpq_poly_clear (FLINTB);
205  return A;
206}
207
208CanonicalForm
209modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
210{
211  CanonicalForm A= F;
212  CanonicalForm B= G;
213
214  fmpq_poly_t FLINTA,FLINTB;
215  convertFacCF2Fmpq_poly_t (FLINTA, A);
216  convertFacCF2Fmpq_poly_t (FLINTB, B);
217
218  fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
219  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
220
221  fmpq_poly_clear (FLINTA);
222  fmpq_poly_clear (FLINTB);
223  return A;
224}
225
226CanonicalForm
227mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G, const Variable& alpha, int m)
228{
229  CanonicalForm A= F;
230  CanonicalForm B= G;
231
232  CanonicalForm denA= bCommonDen (A);
233  CanonicalForm denB= bCommonDen (B);
234
235  A *= denA;
236  B *= denB;
237
238  int degAa= degree (A, alpha);
239  int degBa= degree (B, alpha);
240  int d= degAa + 1 + degBa;
241
242  fmpz_poly_t FLINTA,FLINTB;
243  fmpz_poly_init (FLINTA);
244  fmpz_poly_init (FLINTB);
245  kronSub (FLINTA, A, d);
246  kronSub (FLINTB, B, d);
247
248  int k= d*m;
249  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
250
251  denA *= denB;
252  A= reverseSubst (FLINTA, d, alpha, denA);
253  fmpz_poly_clear (FLINTA);
254  fmpz_poly_clear (FLINTB);
255  return A;
256}
257
258#endif
259
260static
261CFList productsNTL (const CFList& factors, const CanonicalForm& M)
262{
263  zz_p::init (getCharacteristic());
264  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
265  zz_pE::init (NTLMipo);
266  zz_pEX prod;
267  vec_zz_pEX v;
268  v.SetLength (factors.length());
269  int j= 0;
270  for (CFListIterator i= factors; i.hasItem(); i++, j++)
271  {
272    if (i.getItem().inCoeffDomain())
273      v[j]= to_zz_pEX (to_zz_pE (convertFacCF2NTLzzpX (i.getItem())));
274    else
275      v[j]= convertFacCF2NTLzz_pEX (i.getItem(), NTLMipo);
276  }
277  CFList result;
278  Variable x= Variable (1);
279  for (int j= 0; j < factors.length(); j++)
280  {
281    int k= 0;
282    set(prod);
283    for (int i= 0; i < factors.length(); i++, k++)
284    {
285      if (k == j)
286        continue;
287      prod *= v[i];
288    }
289    result.append (convertNTLzz_pEX2CF (prod, x, M.mvar()));
290  }
291  return result;
292}
293
294static
295void tryDiophantine (CFList& result, const CanonicalForm& F,
296                     const CFList& factors, const CanonicalForm& M, bool& fail)
297{
298  ASSERT (M.isUnivariate(), "expected univariate poly");
299
300  CFList bufFactors= factors;
301  bufFactors.removeFirst();
302  bufFactors.insert (factors.getFirst () (0,2));
303  CanonicalForm inv, leadingCoeff= Lc (F);
304  CFListIterator i= bufFactors;
305  if (bufFactors.getFirst().inCoeffDomain())
306  {
307    if (i.hasItem())
308      i++;
309  }
310  for (; i.hasItem(); i++)
311  {
312    tryInvert (Lc (i.getItem()), M, inv ,fail);
313    if (fail)
314      return;
315    i.getItem()= reduce (i.getItem()*inv, M);
316  }
317  bufFactors= productsNTL (bufFactors, M);
318
319  CanonicalForm buf1, buf2, buf3, S, T;
320  i= bufFactors;
321  if (i.hasItem())
322    i++;
323  buf1= bufFactors.getFirst();
324  buf2= i.getItem();
325  tryExtgcd (buf1, buf2, M, buf3, S, T, fail);
326  if (fail)
327    return;
328  result.append (S);
329  result.append (T);
330  if (i.hasItem())
331    i++;
332  for (; i.hasItem(); i++)
333  {
334    buf1= i.getItem();
335    tryExtgcd (buf3, buf1, M, buf3, S, T, fail);
336    if (fail)
337      return;
338    CFListIterator k= factors;
339    for (CFListIterator j= result; j.hasItem(); j++, k++)
340    {
341      j.getItem() *= S;
342      j.getItem()= mod (j.getItem(), k.getItem());
343      j.getItem()= reduce (j.getItem(), M);
344    }
345    result.append (T);
346  }
347}
348
349static inline
350CFList mapinto (const CFList& L)
351{
352  CFList result;
353  for (CFListIterator i= L; i.hasItem(); i++)
354    result.append (mapinto (i.getItem()));
355  return result;
356}
357
358static inline
359int mod (const CFList& L, const CanonicalForm& p)
360{
361  for (CFListIterator i= L; i.hasItem(); i++)
362  {
363    if (mod (i.getItem(), p) == 0)
364      return 0;
365  }
366  return 1;
367}
368
369
370static inline void
371chineseRemainder (const CFList & x1, const CanonicalForm & q1,
372                  const CFList & x2, const CanonicalForm & q2,
373                  CFList & xnew, CanonicalForm & qnew)
374{
375  ASSERT (x1.length() == x2.length(), "expected lists of equal length");
376  CanonicalForm tmp1, tmp2;
377  CFListIterator j= x2;
378  for (CFListIterator i= x1; i.hasItem() && j.hasItem(); i++, j++)
379  {
380    chineseRemainder (i.getItem(), q1, j.getItem(), q2, tmp1, tmp2);
381    xnew.append (tmp1);
382  }
383  qnew= tmp2;
384}
385
386static inline
387CFList Farey (const CFList& L, const CanonicalForm& q)
388{
389  CFList result;
390  for (CFListIterator i= L; i.hasItem(); i++)
391    result.append (Farey (i.getItem(), q));
392  return result;
393}
394
395static inline
396CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
397{
398  CFList result;
399  for (CFListIterator i= L; i.hasItem(); i++)
400    result.append (replacevar (i.getItem(), a, b));
401  return result;
402}
403
404CFList
405modularDiophant (const CanonicalForm& f, const CFList& factors,
406                 const CanonicalForm& M)
407{
408  bool save_rat=!isOn (SW_RATIONAL);
409  On (SW_RATIONAL);
410  CanonicalForm F= f*bCommonDen (f);
411  CFList products= factors;
412  for (CFListIterator i= products; i.hasItem(); i++)
413  {
414    if (products.getFirst().level() == 1)
415      i.getItem() /= Lc (i.getItem());
416    i.getItem() *= bCommonDen (i.getItem());
417  }
418  if (products.getFirst().level() == 1)
419    products.insert (Lc (F));
420  CanonicalForm bound= maxNorm (F);
421  CFList leadingCoeffs;
422  leadingCoeffs.append (lc (F));
423  CanonicalForm dummy;
424  for (CFListIterator i= products; i.hasItem(); i++)
425  {
426    leadingCoeffs.append (lc (i.getItem()));
427    dummy= maxNorm (i.getItem());
428    bound= (dummy > bound) ? dummy : bound;
429  }
430  bound *= maxNorm (Lc (F))*maxNorm (Lc(F))*bound;
431  bound *= bound*bound;
432  bound= power (bound, degree (M));
433  bound *= power (CanonicalForm (2),degree (f));
434  CanonicalForm bufBound= bound;
435  int i = cf_getNumBigPrimes() - 1;
436  int p;
437  CFList resultModP, result, newResult;
438  CanonicalForm q (0), newQ;
439  bool fail= false;
440  Variable a= M.mvar();
441  Variable b= Variable (2);
442  setReduce (M.mvar(), false);
443  CanonicalForm mipo= bCommonDen (M)*M;
444  Off (SW_RATIONAL);
445  CanonicalForm modMipo;
446  leadingCoeffs.append (lc (mipo));
447  CFList tmp1, tmp2;
448  bool equal= false;
449  int count= 0;
450  do
451  {
452    p = cf_getBigPrime( i );
453    i--;
454    while ( i >= 0 && mod( leadingCoeffs, p ) == 0)
455    {
456      p = cf_getBigPrime( i );
457      i--;
458    }
459
460    ASSERT (i >= 0, "ran out of primes"); //sic
461
462    setCharacteristic (p);
463    modMipo= mapinto (mipo);
464    modMipo /= lc (modMipo);
465    resultModP= CFList();
466    tryDiophantine (resultModP, mapinto (F), mapinto (products), modMipo, fail);
467    setCharacteristic (0);
468    if (fail)
469    {
470      fail= false;
471      continue;
472    }
473
474    if ( q.isZero() )
475    {
476      result= replacevar (mapinto(resultModP), a, b);
477      q= p;
478    }
479    else
480    {
481      result= replacevar (result, a, b);
482      newResult= CFList();
483      chineseRemainder( result, q, replacevar (mapinto (resultModP), a, b),
484                        p, newResult, newQ );
485      q= newQ;
486      result= newResult;
487      if (newQ > bound)
488      {
489        count++;
490        tmp1= replacevar (Farey (result, q), b, a);
491        if (tmp2.isEmpty())
492          tmp2= tmp1;
493        else
494        {
495          equal= true;
496          CFListIterator k= tmp1;
497          for (CFListIterator j= tmp2; j.hasItem(); j++, k++)
498          {
499            if (j.getItem() != k.getItem())
500              equal= false;
501          }
502          if (!equal)
503            tmp2= tmp1;
504        }
505        if (count > 2)
506        {
507          bound *= bufBound;
508          equal= false;
509          count= 0;
510        }
511      }
512      if (newQ > bound && equal)
513      {
514        On( SW_RATIONAL );
515        CFList bufResult= result;
516        result= tmp2;
517        setReduce (M.mvar(), true);
518        if (factors.getFirst().level() == 1)
519        {
520          result.removeFirst();
521          CFListIterator j= factors;
522          CanonicalForm denf= bCommonDen (f);
523          for (CFListIterator i= result; i.hasItem(); i++, j++)
524            i.getItem() *= Lc (j.getItem())*denf;
525        }
526        if (factors.getFirst().level() != 1 && 
527            !bCommonDen (factors.getFirst()).isOne())
528        {
529          CanonicalForm denFirst= bCommonDen (factors.getFirst());
530          for (CFListIterator i= result; i.hasItem(); i++)
531            i.getItem() *= denFirst;
532        }
533
534        CanonicalForm test= 0;
535        CFListIterator jj= factors;
536        for (CFListIterator ii= result; ii.hasItem(); ii++, jj++)
537          test += ii.getItem()*(f/jj.getItem());
538        if (!test.isOne())
539        {
540          bound *= bufBound;
541          equal= false;
542          count= 0;
543          setReduce (M.mvar(), false);
544          result= bufResult;
545          Off (SW_RATIONAL);
546        }
547        else
548          break;
549      }
550    }
551  } while (1);
552  if (save_rat) Off(SW_RATIONAL);
553  return result;
554}
555
556CanonicalForm
557mulNTL (const CanonicalForm& F, const CanonicalForm& G)
558{
559  if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
560  {
561    Variable alpha;
562#ifdef HAVE_FLINT
563    if ((!F.inCoeffDomain() && !G.inCoeffDomain()) && (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
564    {
565      CanonicalForm result= mulFLINTQa (F, G, alpha);
566      CanonicalForm test= F*G;
567      if (test != result)
568        printf ("FAILLLLLLLL\n");
569      return result;
570    }
571    else if (!F.inCoeffDomain() && !G.inCoeffDomain())
572    {
573      CanonicalForm result= mulFLINTQ (F,G);
574      CanonicalForm test= F*G;
575      if (test != result)
576        printf ("FAILLLLLLLL2\n");
577      return mulFLINTQ (F, G);
578    }
579#endif
580    return F*G;
581  }
582  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
583  ASSERT (F.level() == G.level(), "expected polys of same level");
584  if (CFFactory::gettype() == GaloisFieldDomain)
585    return F*G;
586  zz_p::init (getCharacteristic());
587  Variable alpha;
588  CanonicalForm result;
589  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
590  {
591    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
592    zz_pE::init (NTLMipo);
593    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
594    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
595    mul (NTLF, NTLF, NTLG);
596    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
597  }
598  else
599  {
600#ifdef HAVE_FLINT
601    nmod_poly_t FLINTF, FLINTG;
602    convertFacCF2nmod_poly_t (FLINTF, F);
603    convertFacCF2nmod_poly_t (FLINTG, G);
604    nmod_poly_mul (FLINTF, FLINTF, FLINTG);
605    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
606    nmod_poly_clear (FLINTF);
607    nmod_poly_clear (FLINTG);
608#else
609    zz_pX NTLF= convertFacCF2NTLzzpX (F);
610    zz_pX NTLG= convertFacCF2NTLzzpX (G);
611    mul (NTLF, NTLF, NTLG);
612    result= convertNTLzzpX2CF(NTLF, F.mvar());
613#endif
614  }
615  return result;
616}
617
618CanonicalForm
619modNTL (const CanonicalForm& F, const CanonicalForm& G)
620{
621  if (F.inCoeffDomain() && G.isUnivariate())
622    return F;
623  else if (F.inCoeffDomain() && G.inCoeffDomain())
624    return mod (F, G);
625  else if (F.isUnivariate() && G.inCoeffDomain())
626    return mod (F,G);
627
628  if (getCharacteristic() == 0)
629  {
630#ifdef HAVE_FLINT
631    Variable alpha;
632    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
633    {
634      CanonicalForm result= modFLINTQ (F,G);
635      CanonicalForm test= mod (F, G);
636      if (test != result)
637        printf ("FAILINMOD\n");
638      return result;
639      //return modFLINTQ (F, G);
640    }
641    else
642     return mod (F, G);
643#else
644    return mod (F, G);
645#endif
646  }
647
648  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
649  ASSERT (F.level() == G.level(), "expected polys of same level");
650  if (CFFactory::gettype() == GaloisFieldDomain)
651    return mod (F, G);
652  zz_p::init (getCharacteristic());
653  Variable alpha;
654  CanonicalForm result;
655  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
656  {
657    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
658    zz_pE::init (NTLMipo);
659    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
660    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
661    rem (NTLF, NTLF, NTLG);
662    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
663  }
664  else
665  {
666#ifdef HAVE_FLINT
667    nmod_poly_t FLINTF, FLINTG;
668    convertFacCF2nmod_poly_t (FLINTF, F);
669    convertFacCF2nmod_poly_t (FLINTG, G);
670    nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
671    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
672    nmod_poly_clear (FLINTF);
673    nmod_poly_clear (FLINTG);
674#else
675    zz_pX NTLF= convertFacCF2NTLzzpX (F);
676    zz_pX NTLG= convertFacCF2NTLzzpX (G);
677    rem (NTLF, NTLF, NTLG);
678    result= convertNTLzzpX2CF(NTLF, F.mvar());
679#endif
680  }
681  return result;
682}
683
684CanonicalForm
685divNTL (const CanonicalForm& F, const CanonicalForm& G)
686{
687  if (F.inCoeffDomain() && G.isUnivariate())
688    return F;
689  else if (F.inCoeffDomain() && G.inCoeffDomain())
690    return div (F, G);
691  else if (F.isUnivariate() && G.inCoeffDomain())
692    return div (F,G);
693
694  if (getCharacteristic() == 0)
695  {
696#ifdef HAVE_FLINT
697    Variable alpha;
698    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
699    {
700      CanonicalForm result= divFLINTQ (F,G);
701      CanonicalForm test= div (F, G);
702      if (test != result)
703        printf ("FAILINDIV\n");
704      return result;
705      //return divFLINTQ (F,G);
706    }
707    else
708      return div (F, G);
709#else
710    return div (F, G);
711#endif
712  }
713
714  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
715  ASSERT (F.level() == G.level(), "expected polys of same level");
716  if (CFFactory::gettype() == GaloisFieldDomain)
717    return div (F, G);
718  zz_p::init (getCharacteristic());
719  Variable alpha;
720  CanonicalForm result;
721  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
722  {
723    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
724    zz_pE::init (NTLMipo);
725    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
726    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
727    div (NTLF, NTLF, NTLG);
728    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
729  }
730  else
731  {
732#ifdef HAVE_FLINT
733    nmod_poly_t FLINTF, FLINTG;
734    convertFacCF2nmod_poly_t (FLINTF, F);
735    convertFacCF2nmod_poly_t (FLINTG, G);
736    nmod_poly_div (FLINTF, FLINTF, FLINTG);
737    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
738    nmod_poly_clear (FLINTF);
739    nmod_poly_clear (FLINTG);
740#else
741    zz_pX NTLF= convertFacCF2NTLzzpX (F);
742    zz_pX NTLG= convertFacCF2NTLzzpX (G);
743    div (NTLF, NTLF, NTLG);
744    result= convertNTLzzpX2CF(NTLF, F.mvar());
745#endif
746  }
747  return result;
748}
749
750/*
751void
752divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
753           CanonicalForm& R)
754{
755  if (F.inCoeffDomain() && G.isUnivariate())
756  {
757    R= F;
758    Q= 0;
759  }
760  else if (F.inCoeffDomain() && G.inCoeffDomain())
761  {
762    divrem (F, G, Q, R);
763    return;
764  }
765  else if (F.isUnivariate() && G.inCoeffDomain())
766  {
767    divrem (F, G, Q, R);
768    return;
769  }
770
771  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
772  ASSERT (F.level() == G.level(), "expected polys of same level");
773  if (CFFactory::gettype() == GaloisFieldDomain)
774  {
775    divrem (F, G, Q, R);
776    return;
777  }
778  zz_p::init (getCharacteristic());
779  Variable alpha;
780  CanonicalForm result;
781  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
782  {
783    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
784    zz_pE::init (NTLMipo);
785    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
786    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
787    zz_pEX NTLQ;
788    zz_pEX NTLR;
789    DivRem (NTLQ, NTLR, NTLF, NTLG);
790    Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
791    R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
792    return;
793  }
794  else
795  {
796    zz_pX NTLF= convertFacCF2NTLzzpX (F);
797    zz_pX NTLG= convertFacCF2NTLzzpX (G);
798    zz_pX NTLQ;
799    zz_pX NTLR;
800    DivRem (NTLQ, NTLR, NTLF, NTLG);
801    Q= convertNTLzzpX2CF(NTLQ, F.mvar());
802    R= convertNTLzzpX2CF(NTLR, F.mvar());
803    return;
804  }
805}*/
806
807CanonicalForm mod (const CanonicalForm& F, const CFList& M)
808{
809  CanonicalForm A= F;
810  for (CFListIterator i= M; i.hasItem(); i++)
811    A= mod (A, i.getItem());
812  return A;
813}
814
815#ifdef HAVE_FLINT
816void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
817{
818  int degAy= degree (A);
819  nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
820
821  nmod_poly_t buf;
822
823  for (CFIterator i= A; i.hasTerms(); i++)
824  {
825    convertFacCF2nmod_poly_t (buf, i.coeff());
826
827    int k= i.exp()*d;
828    int bufRepLength= (int) nmod_poly_length (buf);
829    for (int j= 0; j < bufRepLength; j++)
830      nmod_poly_set_coeff_ui (result, j + k, nmod_poly_get_coeff_ui (buf, j));
831    nmod_poly_clear (buf);
832  }
833  _nmod_poly_normalise (result);
834}
835
836/*void kronSubQ (fmpz_poly_t result, const CanonicalForm& A, int d)
837{
838  int degAy= degree (A);
839  fmpz_poly_init2 (result, d*(degAy + 1));
840  _fmpz_poly_set_length (result, d*(degAy+1));
841
842  CFIterator j;
843  for (CFIterator i= A; i.hasTerms(); i++)
844  {
845    if (i.coeff().inBas
846    convertFacCF2Fmpz_poly_t (buf, i.coeff());
847
848    int k= i.exp()*d;
849    int bufRepLength= (int) fmpz_poly_length (buf);
850    for (int j= 0; j < bufRepLength; j++)
851    {
852      fmpz_poly_get_coeff_fmpz (coeff, buf, j);
853      fmpz_poly_set_coeff_fmpz (result, j + k, coeff);
854    }
855    fmpz_poly_clear (buf);
856  }
857  fmpz_clear (coeff);
858  _fmpz_poly_normalise (result);
859}*/
860
861// A is a bivariate poly over Qa!!!!
862// d2= 2*deg_alpha + 1
863// d1= 2*deg_x*d2+1
864void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
865{
866  int degAy= degree (A);
867  fmpq_poly_init2 (result, d1*(degAy + 1));
868
869  fmpq_poly_t buf;
870  fmpq_t coeff;
871
872  int k, l, bufRepLength;
873  CFIterator j;
874  for (CFIterator i= A; i.hasTerms(); i++)
875  {
876    if (i.coeff().inCoeffDomain())
877    {
878      k= d1*i.exp();
879      convertFacCF2Fmpq_poly_t (buf, i.coeff());
880      bufRepLength= (int) fmpq_poly_length(buf);
881      for (l= 0; l < bufRepLength; l++)
882      {
883        fmpq_poly_get_coeff_fmpq (coeff, buf, l);
884        fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
885      }
886      fmpq_poly_clear (buf);
887    }
888    else
889    {
890      for (j= i.coeff(); j.hasTerms(); j++)
891      {
892        k= d1*i.exp();
893        k += d2*j.exp();
894        convertFacCF2Fmpq_poly_t (buf, j.coeff());
895        bufRepLength= (int) fmpq_poly_length(buf);
896        for (l= 0; l < bufRepLength; l++)
897        {
898          fmpq_poly_get_coeff_fmpq (coeff, buf, l);
899          fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
900        }
901        fmpq_poly_clear (buf);
902      }
903    }
904  }
905  fmpq_clear (coeff);
906  _fmpq_poly_normalise (result);
907}
908#endif
909
910zz_pX kronSubFp (const CanonicalForm& A, int d)
911{
912  int degAy= degree (A);
913  zz_pX result;
914  result.rep.SetLength (d*(degAy + 1));
915
916  zz_p *resultp;
917  resultp= result.rep.elts();
918  zz_pX buf;
919  zz_p *bufp;
920
921  for (CFIterator i= A; i.hasTerms(); i++)
922  {
923    if (i.coeff().inCoeffDomain())
924      buf= convertFacCF2NTLzzpX (i.coeff());
925    else
926      buf= convertFacCF2NTLzzpX (i.coeff());
927
928    int k= i.exp()*d;
929    bufp= buf.rep.elts();
930    int bufRepLength= (int) buf.rep.length();
931    for (int j= 0; j < bufRepLength; j++)
932      resultp [j + k]= bufp [j];
933  }
934  result.normalize();
935
936  return result;
937}
938
939zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
940{
941  int degAy= degree (A);
942  zz_pEX result;
943  result.rep.SetLength (d*(degAy + 1));
944
945  Variable v;
946  zz_pE *resultp;
947  resultp= result.rep.elts();
948  zz_pEX buf1;
949  zz_pE *buf1p;
950  zz_pX buf2;
951  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
952
953  for (CFIterator i= A; i.hasTerms(); i++)
954  {
955    if (i.coeff().inCoeffDomain())
956    {
957      buf2= convertFacCF2NTLzzpX (i.coeff());
958      buf1= to_zz_pEX (to_zz_pE (buf2));
959    }
960    else
961      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
962
963    int k= i.exp()*d;
964    buf1p= buf1.rep.elts();
965    int buf1RepLength= (int) buf1.rep.length();
966    for (int j= 0; j < buf1RepLength; j++)
967      resultp [j + k]= buf1p [j];
968  }
969  result.normalize();
970
971  return result;
972}
973
974void
975kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
976                const Variable& alpha)
977{
978  int degAy= degree (A);
979  subA1.rep.SetLength ((long) d*(degAy + 2));
980  subA2.rep.SetLength ((long) d*(degAy + 2));
981
982  Variable v;
983  zz_pE *subA1p;
984  zz_pE *subA2p;
985  subA1p= subA1.rep.elts();
986  subA2p= subA2.rep.elts();
987  zz_pEX buf;
988  zz_pE *bufp;
989  zz_pX buf2;
990  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
991
992  for (CFIterator i= A; i.hasTerms(); i++)
993  {
994    if (i.coeff().inCoeffDomain())
995    {
996      buf2= convertFacCF2NTLzzpX (i.coeff());
997      buf= to_zz_pEX (to_zz_pE (buf2));
998    }
999    else
1000      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
1001
1002    int k= i.exp()*d;
1003    int kk= (degAy - i.exp())*d;
1004    bufp= buf.rep.elts();
1005    int bufRepLength= (int) buf.rep.length();
1006    for (int j= 0; j < bufRepLength; j++)
1007    {
1008      subA1p [j + k] += bufp [j];
1009      subA2p [j + kk] += bufp [j];
1010    }
1011  }
1012  subA1.normalize();
1013  subA2.normalize();
1014}
1015
1016void
1017kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
1018{
1019  int degAy= degree (A);
1020  subA1.rep.SetLength ((long) d*(degAy + 2));
1021  subA2.rep.SetLength ((long) d*(degAy + 2));
1022
1023  zz_p *subA1p;
1024  zz_p *subA2p;
1025  subA1p= subA1.rep.elts();
1026  subA2p= subA2.rep.elts();
1027  zz_pX buf;
1028  zz_p *bufp;
1029
1030  for (CFIterator i= A; i.hasTerms(); i++)
1031  {
1032    buf= convertFacCF2NTLzzpX (i.coeff());
1033
1034    int k= i.exp()*d;
1035    int kk= (degAy - i.exp())*d;
1036    bufp= buf.rep.elts();
1037    int bufRepLength= (int) buf.rep.length();
1038    for (int j= 0; j < bufRepLength; j++)
1039    {
1040      subA1p [j + k] += bufp [j];
1041      subA2p [j + kk] += bufp [j];
1042    }
1043  }
1044  subA1.normalize();
1045  subA2.normalize();
1046}
1047
1048CanonicalForm
1049reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
1050              const Variable& alpha)
1051{
1052  Variable y= Variable (2);
1053  Variable x= Variable (1);
1054
1055  zz_pEX f= F;
1056  zz_pEX g= G;
1057  int degf= deg(f);
1058  int degg= deg(g);
1059
1060  zz_pEX buf1;
1061  zz_pEX buf2;
1062  zz_pEX buf3;
1063
1064  zz_pE *buf1p;
1065  zz_pE *buf2p;
1066  zz_pE *buf3p;
1067  if (f.rep.length() < (long) d*(k+1)) //zero padding
1068    f.rep.SetLength ((long)d*(k+1));
1069
1070  zz_pE *gp= g.rep.elts();
1071  zz_pE *fp= f.rep.elts();
1072  CanonicalForm result= 0;
1073  int i= 0;
1074  int lf= 0;
1075  int lg= d*k;
1076  int degfSubLf= degf;
1077  int deggSubLg= degg-lg;
1078  int repLengthBuf2;
1079  int repLengthBuf1;
1080  zz_pE zzpEZero= zz_pE();
1081
1082  while (degf >= lf || lg >= 0)
1083  {
1084    if (degfSubLf >= d)
1085      repLengthBuf1= d;
1086    else if (degfSubLf < 0)
1087      repLengthBuf1= 0;
1088    else
1089      repLengthBuf1= degfSubLf + 1;
1090    buf1.rep.SetLength((long) repLengthBuf1);
1091
1092    buf1p= buf1.rep.elts();
1093    for (int ind= 0; ind < repLengthBuf1; ind++)
1094      buf1p [ind]= fp [ind + lf];
1095    buf1.normalize();
1096
1097    repLengthBuf1= buf1.rep.length();
1098
1099    if (deggSubLg >= d - 1)
1100      repLengthBuf2= d - 1;
1101    else if (deggSubLg < 0)
1102      repLengthBuf2= 0;
1103    else
1104      repLengthBuf2= deggSubLg + 1;
1105
1106    buf2.rep.SetLength ((long) repLengthBuf2);
1107    buf2p= buf2.rep.elts();
1108    for (int ind= 0; ind < repLengthBuf2; ind++)
1109    {
1110      buf2p [ind]= gp [ind + lg];
1111    }
1112    buf2.normalize();
1113
1114    repLengthBuf2= buf2.rep.length();
1115
1116    buf3.rep.SetLength((long) repLengthBuf2 + d);
1117    buf3p= buf3.rep.elts();
1118    buf2p= buf2.rep.elts();
1119    buf1p= buf1.rep.elts();
1120    for (int ind= 0; ind < repLengthBuf1; ind++)
1121      buf3p [ind]= buf1p [ind];
1122    for (int ind= repLengthBuf1; ind < d; ind++)
1123      buf3p [ind]= zzpEZero;
1124    for (int ind= 0; ind < repLengthBuf2; ind++)
1125      buf3p [ind + d]= buf2p [ind];
1126    buf3.normalize();
1127
1128    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
1129    i++;
1130
1131
1132    lf= i*d;
1133    degfSubLf= degf - lf;
1134
1135    lg= d*(k-i);
1136    deggSubLg= degg - lg;
1137
1138    buf1p= buf1.rep.elts();
1139
1140    if (lg >= 0 && deggSubLg > 0)
1141    {
1142      if (repLengthBuf2 > degfSubLf + 1)
1143        degfSubLf= repLengthBuf2 - 1;
1144      int tmp= tmin (repLengthBuf1, deggSubLg + 1);
1145      for (int ind= 0; ind < tmp; ind++)
1146        gp [ind + lg] -= buf1p [ind];
1147    }
1148
1149    if (lg < 0)
1150      break;
1151
1152    buf2p= buf2.rep.elts();
1153    if (degfSubLf >= 0)
1154    {
1155      for (int ind= 0; ind < repLengthBuf2; ind++)
1156        fp [ind + lf] -= buf2p [ind];
1157    }
1158  }
1159
1160  return result;
1161}
1162
1163CanonicalForm
1164reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
1165{
1166  Variable y= Variable (2);
1167  Variable x= Variable (1);
1168
1169  zz_pX f= F;
1170  zz_pX g= G;
1171  int degf= deg(f);
1172  int degg= deg(g);
1173
1174  zz_pX buf1;
1175  zz_pX buf2;
1176  zz_pX buf3;
1177
1178  zz_p *buf1p;
1179  zz_p *buf2p;
1180  zz_p *buf3p;
1181
1182  if (f.rep.length() < (long) d*(k+1)) //zero padding
1183    f.rep.SetLength ((long)d*(k+1));
1184
1185  zz_p *gp= g.rep.elts();
1186  zz_p *fp= f.rep.elts();
1187  CanonicalForm result= 0;
1188  int i= 0;
1189  int lf= 0;
1190  int lg= d*k;
1191  int degfSubLf= degf;
1192  int deggSubLg= degg-lg;
1193  int repLengthBuf2;
1194  int repLengthBuf1;
1195  zz_p zzpZero= zz_p();
1196  while (degf >= lf || lg >= 0)
1197  {
1198    if (degfSubLf >= d)
1199      repLengthBuf1= d;
1200    else if (degfSubLf < 0)
1201      repLengthBuf1= 0;
1202    else
1203      repLengthBuf1= degfSubLf + 1;
1204    buf1.rep.SetLength((long) repLengthBuf1);
1205
1206    buf1p= buf1.rep.elts();
1207    for (int ind= 0; ind < repLengthBuf1; ind++)
1208      buf1p [ind]= fp [ind + lf];
1209    buf1.normalize();
1210
1211    repLengthBuf1= buf1.rep.length();
1212
1213    if (deggSubLg >= d - 1)
1214      repLengthBuf2= d - 1;
1215    else if (deggSubLg < 0)
1216      repLengthBuf2= 0;
1217    else
1218      repLengthBuf2= deggSubLg + 1;
1219
1220    buf2.rep.SetLength ((long) repLengthBuf2);
1221    buf2p= buf2.rep.elts();
1222    for (int ind= 0; ind < repLengthBuf2; ind++)
1223      buf2p [ind]= gp [ind + lg];
1224
1225    buf2.normalize();
1226
1227    repLengthBuf2= buf2.rep.length();
1228
1229
1230    buf3.rep.SetLength((long) repLengthBuf2 + d);
1231    buf3p= buf3.rep.elts();
1232    buf2p= buf2.rep.elts();
1233    buf1p= buf1.rep.elts();
1234    for (int ind= 0; ind < repLengthBuf1; ind++)
1235      buf3p [ind]= buf1p [ind];
1236    for (int ind= repLengthBuf1; ind < d; ind++)
1237      buf3p [ind]= zzpZero;
1238    for (int ind= 0; ind < repLengthBuf2; ind++)
1239      buf3p [ind + d]= buf2p [ind];
1240    buf3.normalize();
1241
1242    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
1243    i++;
1244
1245
1246    lf= i*d;
1247    degfSubLf= degf - lf;
1248
1249    lg= d*(k-i);
1250    deggSubLg= degg - lg;
1251
1252    buf1p= buf1.rep.elts();
1253
1254    if (lg >= 0 && deggSubLg > 0)
1255    {
1256      if (repLengthBuf2 > degfSubLf + 1)
1257        degfSubLf= repLengthBuf2 - 1;
1258      int tmp= tmin (repLengthBuf1, deggSubLg + 1);
1259      for (int ind= 0; ind < tmp; ind++)
1260        gp [ind + lg] -= buf1p [ind];
1261    }
1262    if (lg < 0)
1263      break;
1264
1265    buf2p= buf2.rep.elts();
1266    if (degfSubLf >= 0)
1267    {
1268      for (int ind= 0; ind < repLengthBuf2; ind++)
1269        fp [ind + lf] -= buf2p [ind];
1270    }
1271  }
1272
1273  return result;
1274}
1275
1276CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
1277{
1278  Variable y= Variable (2);
1279  Variable x= Variable (1);
1280
1281  zz_pEX f= F;
1282  zz_pE *fp= f.rep.elts();
1283
1284  zz_pEX buf;
1285  zz_pE *bufp;
1286  CanonicalForm result= 0;
1287  int i= 0;
1288  int degf= deg(f);
1289  int k= 0;
1290  int degfSubK;
1291  int repLength;
1292  while (degf >= k)
1293  {
1294    degfSubK= degf - k;
1295    if (degfSubK >= d)
1296      repLength= d;
1297    else
1298      repLength= degfSubK + 1;
1299
1300    buf.rep.SetLength ((long) repLength);
1301    bufp= buf.rep.elts();
1302    for (int j= 0; j < repLength; j++)
1303      bufp [j]= fp [j + k];
1304    buf.normalize();
1305
1306    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
1307    i++;
1308    k= d*i;
1309  }
1310
1311  return result;
1312}
1313
1314CanonicalForm reverseSubstFp (const zz_pX& F, int d)
1315{
1316  Variable y= Variable (2);
1317  Variable x= Variable (1);
1318
1319  zz_pX f= F;
1320  zz_p *fp= f.rep.elts();
1321
1322  zz_pX buf;
1323  zz_p *bufp;
1324  CanonicalForm result= 0;
1325  int i= 0;
1326  int degf= deg(f);
1327  int k= 0;
1328  int degfSubK;
1329  int repLength;
1330  while (degf >= k)
1331  {
1332    degfSubK= degf - k;
1333    if (degfSubK >= d)
1334      repLength= d;
1335    else
1336      repLength= degfSubK + 1;
1337
1338    buf.rep.SetLength ((long) repLength);
1339    bufp= buf.rep.elts();
1340    for (int j= 0; j < repLength; j++)
1341      bufp [j]= fp [j + k];
1342    buf.normalize();
1343
1344    result += convertNTLzzpX2CF (buf, x)*power (y, i);
1345    i++;
1346    k= d*i;
1347  }
1348
1349  return result;
1350}
1351
1352// assumes input to be reduced mod M and to be an element of Fq not Fp
1353CanonicalForm
1354mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
1355                  CanonicalForm& M)
1356{
1357  int d1= degree (F, 1) + degree (G, 1) + 1;
1358  d1 /= 2;
1359  d1 += 1;
1360
1361  zz_pX F1, F2;
1362  kronSubRecipro (F1, F2, F, d1);
1363  zz_pX G1, G2;
1364  kronSubRecipro (G1, G2, G, d1);
1365
1366  int k= d1*degree (M);
1367  MulTrunc (F1, F1, G1, (long) k);
1368
1369  int degtailF= degree (tailcoeff (F), 1);
1370  int degtailG= degree (tailcoeff (G), 1);
1371  int taildegF= taildegree (F);
1372  int taildegG= taildegree (G);
1373  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
1374
1375  reverse (F2, F2);
1376  reverse (G2, G2);
1377  MulTrunc (F2, F2, G2, b + 1);
1378  reverse (F2, F2, b);
1379
1380  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
1381  return reverseSubst (F1, F2, d1, d2);
1382}
1383
1384//Kronecker substitution
1385CanonicalForm
1386mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
1387              CanonicalForm& M)
1388{
1389  CanonicalForm A= F;
1390  CanonicalForm B= G;
1391
1392  int degAx= degree (A, 1);
1393  int degAy= degree (A, 2);
1394  int degBx= degree (B, 1);
1395  int degBy= degree (B, 2);
1396  int d1= degAx + 1 + degBx;
1397  int d2= tmax (degAy, degBy);
1398
1399  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
1400    return mulMod2NTLFpReci (A, B, M);
1401
1402  zz_pX NTLA= kronSubFp (A, d1);
1403  zz_pX NTLB= kronSubFp (B, d1);
1404
1405  int k= d1*degree (M);
1406  MulTrunc (NTLA, NTLA, NTLB, (long) k);
1407
1408  A= reverseSubstFp (NTLA, d1);
1409
1410  return A;
1411}
1412
1413// assumes input to be reduced mod M and to be an element of Fq not Fp
1414CanonicalForm
1415mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
1416                  CanonicalForm& M, const Variable& alpha)
1417{
1418  int d1= degree (F, 1) + degree (G, 1) + 1;
1419  d1 /= 2;
1420  d1 += 1;
1421
1422  zz_pEX F1, F2;
1423  kronSubRecipro (F1, F2, F, d1, alpha);
1424  zz_pEX G1, G2;
1425  kronSubRecipro (G1, G2, G, d1, alpha);
1426
1427  int k= d1*degree (M);
1428  MulTrunc (F1, F1, G1, (long) k);
1429
1430  int degtailF= degree (tailcoeff (F), 1);
1431  int degtailG= degree (tailcoeff (G), 1);
1432  int taildegF= taildegree (F);
1433  int taildegG= taildegree (G);
1434  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
1435
1436  reverse (F2, F2);
1437  reverse (G2, G2);
1438  MulTrunc (F2, F2, G2, b + 1);
1439  reverse (F2, F2, b);
1440
1441  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
1442  return reverseSubst (F1, F2, d1, d2, alpha);
1443}
1444
1445#ifdef HAVE_FLINT
1446CanonicalForm
1447mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
1448                CanonicalForm& M);
1449#endif
1450
1451CanonicalForm
1452mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
1453              CanonicalForm& M)
1454{
1455  Variable alpha;
1456  CanonicalForm A= F;
1457  CanonicalForm B= G;
1458
1459  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
1460  {
1461    int degAx= degree (A, 1);
1462    int degAy= degree (A, 2);
1463    int degBx= degree (B, 1);
1464    int degBy= degree (B, 2);
1465    int d1= degAx + degBx + 1;
1466    int d2= tmax (degAy, degBy);
1467    zz_p::init (getCharacteristic());
1468    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1469    zz_pE::init (NTLMipo);
1470
1471    int degMipo= degree (getMipo (alpha));
1472    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
1473        (2*degAy > degree (M)))
1474      return mulMod2NTLFqReci (A, B, M, alpha);
1475
1476    zz_pEX NTLA= kronSub (A, d1, alpha);
1477    zz_pEX NTLB= kronSub (B, d1, alpha);
1478
1479    int k= d1*degree (M);
1480
1481    MulTrunc (NTLA, NTLA, NTLB, (long) k);
1482
1483    A= reverseSubst (NTLA, d1, alpha);
1484
1485    return A;
1486  }
1487  else
1488#ifdef HAVE_FLINT
1489    return mulMod2FLINTFp (A, B, M);
1490#else
1491    return mulMod2NTLFp (A, B, M);
1492#endif
1493}
1494
1495#ifdef HAVE_FLINT
1496void
1497kronSubRecipro (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A, int d)
1498{
1499  int degAy= degree (A);
1500  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1501  nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
1502  nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
1503
1504  nmod_poly_t buf;
1505
1506  for (CFIterator i= A; i.hasTerms(); i++)
1507  {
1508    convertFacCF2nmod_poly_t (buf, i.coeff());
1509
1510    int k= i.exp()*d;
1511    int kk= (degAy - i.exp())*d;
1512    int bufRepLength= (int) nmod_poly_length (buf);
1513    for (int j= 0; j < bufRepLength; j++)
1514    {
1515      nmod_poly_set_coeff_ui (subA1, j + k, n_addmod (nmod_poly_get_coeff_ui (subA1, j+k), nmod_poly_get_coeff_ui (buf, j), getCharacteristic())); //TODO muß man subA1 coeffs von subA1 vorher auf null setzen??
1516      nmod_poly_set_coeff_ui (subA2, j + kk, n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),nmod_poly_get_coeff_ui (buf, j), getCharacteristic()));
1517    }
1518    nmod_poly_clear (buf);
1519  }
1520  _nmod_poly_normalise (subA1);
1521  _nmod_poly_normalise (subA2);
1522}
1523
1524void
1525kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A, int d)
1526{
1527  int degAy= degree (A);
1528  fmpz_poly_init2 (subA1, d*(degAy + 2));
1529  fmpz_poly_init2 (subA2, d*(degAy + 2));
1530
1531  fmpz_poly_t buf;
1532  fmpz_t coeff1, coeff2;
1533
1534  for (CFIterator i= A; i.hasTerms(); i++)
1535  {
1536    convertFacCF2Fmpz_poly_t (buf, i.coeff());
1537
1538    int k= i.exp()*d;
1539    int kk= (degAy - i.exp())*d;
1540    int bufRepLength= (int) fmpz_poly_length (buf);
1541    for (int j= 0; j < bufRepLength; j++)
1542    {
1543      fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k); //TODO muß man subA1 coeffs von subA1 vorher auf null setzen??
1544                                                     // vielleicht ist es schneller fmpz_poly_get_ptr zu nehmen
1545      fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
1546      fmpz_add (coeff1, coeff1, coeff2);
1547      fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
1548      fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
1549      fmpz_add (coeff1, coeff1, coeff2);
1550      fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
1551    }
1552    fmpz_poly_clear (buf);
1553  }
1554  fmpz_clear (coeff1);
1555  fmpz_clear (coeff2);
1556  _fmpz_poly_normalise (subA1);
1557  _fmpz_poly_normalise (subA2);
1558}
1559
1560CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
1561{
1562  Variable y= Variable (2);
1563  Variable x= Variable (1);
1564
1565  fmpz_poly_t f;
1566  fmpz_poly_init (f);
1567  fmpz_poly_set (f, F);
1568
1569  fmpz_poly_t buf;
1570  CanonicalForm result= 0;
1571  int i= 0;
1572  int degf= fmpz_poly_degree(f);
1573  int k= 0;
1574  int degfSubK;
1575  int repLength;
1576  fmpz_t coeff;
1577  while (degf >= k)
1578  {
1579    degfSubK= degf - k;
1580    if (degfSubK >= d)
1581      repLength= d;
1582    else
1583      repLength= degfSubK + 1;
1584
1585    fmpz_poly_init2 (buf, repLength);
1586    fmpz_init (coeff);
1587    for (int j= 0; j < repLength; j++)
1588    {
1589      fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
1590      fmpz_poly_set_coeff_fmpz (buf, j, coeff);
1591    }
1592    _fmpz_poly_normalise (buf);
1593
1594    result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
1595    i++;
1596    k= d*i;
1597    fmpz_poly_clear (buf);
1598    fmpz_clear (coeff);
1599  }
1600  fmpz_poly_clear (f);
1601
1602  return result;
1603}
1604
1605CanonicalForm
1606reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
1607{
1608  Variable y= Variable (2);
1609  Variable x= Variable (1);
1610
1611  nmod_poly_t f, g;
1612  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1613  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
1614  nmod_poly_init_preinv (g, getCharacteristic(), ninv);
1615  nmod_poly_set (f, F);
1616  nmod_poly_set (g, G);
1617  int degf= nmod_poly_degree(f);
1618  int degg= nmod_poly_degree(g);
1619
1620
1621  nmod_poly_t buf1,buf2, buf3;
1622  /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
1623  nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
1624  nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
1625
1626  if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
1627    nmod_poly_fit_length (f,(long)d*(k+1));
1628
1629  CanonicalForm result= 0;
1630  int i= 0;
1631  int lf= 0;
1632  int lg= d*k;
1633  int degfSubLf= degf;
1634  int deggSubLg= degg-lg;
1635  int repLengthBuf2;
1636  int repLengthBuf1;
1637  //zz_p zzpZero= zz_p();
1638  while (degf >= lf || lg >= 0)
1639  {
1640    if (degfSubLf >= d)
1641      repLengthBuf1= d;
1642    else if (degfSubLf < 0)
1643      repLengthBuf1= 0;
1644    else
1645      repLengthBuf1= degfSubLf + 1;
1646    //cout << "repLengthBuf1= " << repLengthBuf1 << "\n";
1647    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
1648    nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
1649    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
1650
1651    for (int ind= 0; ind < repLengthBuf1; ind++)
1652      nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
1653    _nmod_poly_normalise (buf1);
1654    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
1655
1656    repLengthBuf1= nmod_poly_length (buf1);
1657
1658    if (deggSubLg >= d - 1)
1659      repLengthBuf2= d - 1;
1660    else if (deggSubLg < 0)
1661      repLengthBuf2= 0;
1662    else
1663      repLengthBuf2= deggSubLg + 1;
1664
1665    //cout << "repLengthBuf2= " << repLengthBuf2 << "\n";
1666    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
1667    nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
1668    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
1669    for (int ind= 0; ind < repLengthBuf2; ind++)
1670      nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
1671
1672    _nmod_poly_normalise (buf2);
1673    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
1674    repLengthBuf2= nmod_poly_length (buf2);
1675
1676    nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
1677    for (int ind= 0; ind < repLengthBuf1; ind++)
1678      nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
1679    for (int ind= repLengthBuf1; ind < d; ind++)
1680      nmod_poly_set_coeff_ui (buf3, ind, 0);
1681    for (int ind= 0; ind < repLengthBuf2; ind++)
1682      nmod_poly_set_coeff_ui (buf3, ind + d, nmod_poly_get_coeff_ui (buf2, ind));
1683    _nmod_poly_normalise (buf3);
1684
1685    result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
1686    i++;
1687
1688
1689    lf= i*d;
1690    degfSubLf= degf - lf;
1691
1692    lg= d*(k-i);
1693    deggSubLg= degg - lg;
1694
1695    if (lg >= 0 && deggSubLg > 0)
1696    {
1697      if (repLengthBuf2 > degfSubLf + 1)
1698        degfSubLf= repLengthBuf2 - 1;
1699      int tmp= tmin (repLengthBuf1, deggSubLg + 1);
1700      for (int ind= 0; ind < tmp; ind++)
1701      //{
1702        nmod_poly_set_coeff_ui (g, ind + lg, n_submod (nmod_poly_get_coeff_ui (g, ind + lg), nmod_poly_get_coeff_ui (buf1, ind), getCharacteristic()));
1703        //cout << "nmod_poly_get_coeff_ui (g, ind + lg)= " << nmod_poly_get_coeff_ui (g, ind + lg) << "\n";
1704      //}
1705    }
1706    if (lg < 0)
1707    {
1708      nmod_poly_clear (buf1);
1709      nmod_poly_clear (buf2);
1710      nmod_poly_clear (buf3);
1711      break;
1712    }
1713    if (degfSubLf >= 0)
1714    {
1715      for (int ind= 0; ind < repLengthBuf2; ind++)
1716        nmod_poly_set_coeff_ui (f, ind + lf, n_submod (nmod_poly_get_coeff_ui (f, ind + lf), nmod_poly_get_coeff_ui (buf2, ind), getCharacteristic()));
1717    }
1718    nmod_poly_clear (buf1);
1719    nmod_poly_clear (buf2);
1720    nmod_poly_clear (buf3);
1721    /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
1722    nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
1723    nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
1724  }
1725
1726  /*nmod_poly_clear (buf1);
1727  nmod_poly_clear (buf2);
1728  nmod_poly_clear (buf3);*/
1729  nmod_poly_clear (f);
1730  nmod_poly_clear (g);
1731
1732  return result;
1733}
1734
1735CanonicalForm
1736reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
1737{
1738  Variable y= Variable (2);
1739  Variable x= Variable (1);
1740
1741  fmpz_poly_t f, g;
1742  fmpz_poly_init (f);
1743  fmpz_poly_init (g);
1744  fmpz_poly_set (f, F);
1745  fmpz_poly_set (g, G);
1746  int degf= fmpz_poly_degree(f);
1747  int degg= fmpz_poly_degree(g);
1748
1749
1750  fmpz_poly_t buf1,buf2, buf3;
1751  /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
1752  nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
1753  nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
1754
1755  if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
1756    fmpz_poly_fit_length (f,(long)d*(k+1));
1757
1758  CanonicalForm result= 0;
1759  int i= 0;
1760  int lf= 0;
1761  int lg= d*k;
1762  int degfSubLf= degf;
1763  int deggSubLg= degg-lg;
1764  int repLengthBuf2;
1765  int repLengthBuf1;
1766  //zz_p zzpZero= zz_p();
1767  fmpz_t tmp1, tmp2;
1768  while (degf >= lf || lg >= 0)
1769  {
1770    if (degfSubLf >= d)
1771      repLengthBuf1= d;
1772    else if (degfSubLf < 0)
1773      repLengthBuf1= 0;
1774    else
1775      repLengthBuf1= degfSubLf + 1;
1776    //cout << "repLengthBuf1= " << repLengthBuf1 << "\n";
1777    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
1778    fmpz_poly_init2 (buf1, repLengthBuf1);
1779    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
1780
1781    for (int ind= 0; ind < repLengthBuf1; ind++)
1782    {
1783      fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1784      fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
1785    }
1786    _fmpz_poly_normalise (buf1);
1787    //cout << "nmod_poly_length (buf1)= " << nmod_poly_length (buf1) << "\n";
1788
1789    repLengthBuf1= fmpz_poly_length (buf1);
1790
1791    if (deggSubLg >= d - 1)
1792      repLengthBuf2= d - 1;
1793    else if (deggSubLg < 0)
1794      repLengthBuf2= 0;
1795    else
1796      repLengthBuf2= deggSubLg + 1;
1797
1798    //cout << "repLengthBuf2= " << repLengthBuf2 << "\n";
1799    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
1800    fmpz_poly_init2 (buf2, repLengthBuf2);
1801    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
1802    for (int ind= 0; ind < repLengthBuf2; ind++)
1803    {
1804      fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1805      fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
1806    }
1807
1808    _fmpz_poly_normalise (buf2);
1809    //cout << "nmod_poly_length (buf2)= " << nmod_poly_length (buf2) << "\n";
1810    repLengthBuf2= fmpz_poly_length (buf2);
1811
1812    fmpz_poly_init2 (buf3, repLengthBuf2 + d);
1813    for (int ind= 0; ind < repLengthBuf1; ind++)
1814    {
1815      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))
1816      fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
1817    }
1818    for (int ind= repLengthBuf1; ind < d; ind++)
1819      fmpz_poly_set_coeff_ui (buf3, ind, 0);
1820    for (int ind= 0; ind < repLengthBuf2; ind++)
1821    {
1822      fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
1823      fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
1824    }
1825    _fmpz_poly_normalise (buf3);
1826
1827    result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
1828    i++;
1829
1830
1831    lf= i*d;
1832    degfSubLf= degf - lf;
1833
1834    lg= d*(k-i);
1835    deggSubLg= degg - lg;
1836
1837    if (lg >= 0 && deggSubLg > 0)
1838    {
1839      if (repLengthBuf2 > degfSubLf + 1)
1840        degfSubLf= repLengthBuf2 - 1;
1841      int tmp= tmin (repLengthBuf1, deggSubLg + 1);
1842      for (int ind= 0; ind < tmp; ind++)
1843      {
1844        fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1845        fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
1846        fmpz_sub (tmp1, tmp1, tmp2);
1847        fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
1848        //cout << "nmod_poly_get_coeff_ui (g, ind + lg)= " << nmod_poly_get_coeff_ui (g, ind + lg) << "\n";
1849      }
1850    }
1851    if (lg < 0)
1852    {
1853      fmpz_poly_clear (buf1);
1854      fmpz_poly_clear (buf2);
1855      fmpz_poly_clear (buf3);
1856      break;
1857    }
1858    if (degfSubLf >= 0)
1859    {
1860      for (int ind= 0; ind < repLengthBuf2; ind++)
1861      {
1862        fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1863        fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
1864        fmpz_sub (tmp1, tmp1, tmp2);
1865        fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
1866      }
1867    }
1868    fmpz_poly_clear (buf1);
1869    fmpz_poly_clear (buf2);
1870    fmpz_poly_clear (buf3);
1871    /*nmod_poly_init_preinv (buf1, getCharacteristic(), ninv);
1872    nmod_poly_init_preinv (buf2, getCharacteristic(), ninv);
1873    nmod_poly_init_preinv (buf3, getCharacteristic(), ninv);*/
1874  }
1875
1876  /*nmod_poly_clear (buf1);
1877  nmod_poly_clear (buf2);
1878  nmod_poly_clear (buf3);*/
1879  fmpz_poly_clear (f);
1880  fmpz_poly_clear (g);
1881  fmpz_clear (tmp1);
1882  fmpz_clear (tmp2);
1883
1884  return result;
1885}
1886
1887CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
1888{
1889  Variable y= Variable (2);
1890  Variable x= Variable (1);
1891
1892  nmod_poly_t f;
1893  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1894  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
1895  nmod_poly_set (f, F);
1896
1897  nmod_poly_t buf;
1898  //nmod_poly_init_preinv (buf, getCharacteristic(), ninv);
1899  CanonicalForm result= 0;
1900  int i= 0;
1901  int degf= nmod_poly_degree(f);
1902  int k= 0;
1903  int degfSubK;
1904  int repLength;
1905  while (degf >= k)
1906  {
1907    degfSubK= degf - k;
1908    if (degfSubK >= d)
1909      repLength= d;
1910    else
1911      repLength= degfSubK + 1;
1912
1913    nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
1914    for (int j= 0; j < repLength; j++)
1915      nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
1916    //printf ("after set\n");
1917    _nmod_poly_normalise (buf);
1918    //printf ("after normalise\n");
1919
1920    result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
1921    i++;
1922    k= d*i;
1923    nmod_poly_clear (buf);
1924    //nmod_poly_init_preinv (buf, getCharacteristic(), ninv);
1925  }
1926  //printf ("after while\n");
1927  //nmod_poly_clear (buf);
1928  nmod_poly_clear (f);
1929  //printf ("after clear\n");
1930
1931  return result;
1932}
1933
1934CanonicalForm
1935mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
1936                    CanonicalForm& M)
1937{
1938  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
1939  d1 /= 2;
1940  d1 += 1;
1941
1942  nmod_poly_t F1, F2;
1943  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1944  nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
1945  nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
1946  kronSubRecipro (F1, F2, F, d1);
1947  /*zz_pX TF1, TF2;
1948  kronSubRecipro (TF1, TF2, F, d1);
1949  Variable x= Variable (1);
1950  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";
1951  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
1952
1953  nmod_poly_t G1, G2;
1954  nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
1955  nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
1956  kronSubRecipro (G1, G2, G, d1);
1957  /*zz_pX TG1, TG2;
1958  kronSubRecipro (TG1, TG2, G, d1);
1959  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG1, x) == convertnmod_poly_t2FacCF (G1, x)) << "\n";
1960  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG2, x) == convertnmod_poly_t2FacCF (G2, x)) << "\n";*/
1961
1962
1963  int k= d1*degree (M);
1964  nmod_poly_mullow (F1, F1, G1, (long) k);
1965  /*MulTrunc (TF1, TF1, TG1, (long) k);
1966  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";*/
1967
1968  int degtailF= degree (tailcoeff (F), 1);;
1969  int degtailG= degree (tailcoeff (G), 1);
1970  int taildegF= taildegree (F);
1971  int taildegG= taildegree (G);
1972
1973  int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG + d1*(2+taildegF + taildegG);
1974  nmod_poly_mulhigh (F2, F2, G2, b);
1975  nmod_poly_shift_right (F2, F2, b);
1976  /*mul (TF2, TF2, TG2);
1977  if (deg (TF2) > k)
1978  {
1979    int b= deg (TF2) - k + 2;
1980    TF2 >>= b;
1981  }
1982  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
1983  int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
1984
1985
1986  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
1987  //CanonicalForm NTLres= reverseSubst (TF1, TF2, d1, d2);
1988  //cout << "FINALtestkronSubReci= " << (NTLres == result) << "\n";
1989
1990  nmod_poly_clear (F1);
1991  nmod_poly_clear (F2);
1992  nmod_poly_clear (G1);
1993  nmod_poly_clear (G2);
1994  return result;
1995}
1996
1997CanonicalForm
1998mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
1999                CanonicalForm& M)
2000{
2001  CanonicalForm A= F;
2002  CanonicalForm B= G;
2003
2004  int degAx= degree (A, 1);
2005  int degAy= degree (A, 2);
2006  int degBx= degree (B, 1);
2007  int degBy= degree (B, 2);
2008  int d1= degAx + 1 + degBx;
2009  int d2= tmax (degAy, degBy);
2010
2011  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
2012    return mulMod2FLINTFpReci (A, B, M);
2013
2014  nmod_poly_t FLINTA, FLINTB;
2015  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
2016  nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
2017  nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
2018  kronSubFp (FLINTA, A, d1);
2019  kronSubFp (FLINTB, B, d1);
2020
2021  int k= d1*degree (M);
2022  nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
2023
2024  A= reverseSubstFp (FLINTA, d1);
2025
2026  nmod_poly_clear (FLINTA);
2027  nmod_poly_clear (FLINTB);
2028  return A;
2029}
2030
2031CanonicalForm
2032mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
2033                    CanonicalForm& M)
2034{
2035  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
2036  d1 /= 2;
2037  d1 += 1;
2038
2039  fmpz_poly_t F1, F2;
2040  fmpz_poly_init (F1);
2041  fmpz_poly_init (F2);
2042  kronSubRecipro (F1, F2, F, d1);
2043  /*zz_pX TF1, TF2;
2044  kronSubRecipro (TF1, TF2, F, d1);
2045  Variable x= Variable (1);
2046  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";
2047  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
2048
2049  fmpz_poly_t G1, G2;
2050  fmpz_poly_init (G1);
2051  fmpz_poly_init (G2);
2052  kronSubRecipro (G1, G2, G, d1);
2053  /*zz_pX TG1, TG2;
2054  kronSubRecipro (TG1, TG2, G, d1);
2055  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG1, x) == convertnmod_poly_t2FacCF (G1, x)) << "\n";
2056  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TG2, x) == convertnmod_poly_t2FacCF (G2, x)) << "\n";*/
2057
2058
2059  int k= d1*degree (M);
2060  fmpz_poly_mullow (F1, F1, G1, (long) k);
2061  /*MulTrunc (TF1, TF1, TG1, (long) k);
2062  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF1, x) == convertnmod_poly_t2FacCF (F1, x)) << "\n";*/
2063
2064  int degtailF= degree (tailcoeff (F), 1);;
2065  int degtailG= degree (tailcoeff (G), 1);
2066  int taildegF= taildegree (F);
2067  int taildegG= taildegree (G);
2068
2069  int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG + d1*(2+taildegF + taildegG);
2070  fmpz_poly_mulhigh_n (F2, F2, G2, b);
2071  fmpz_poly_shift_right (F2, F2, b);
2072  /*mul (TF2, TF2, TG2);
2073  if (deg (TF2) > k)
2074  {
2075    int b= deg (TF2) - k + 2;
2076    TF2 >>= b;
2077  }
2078  cout << "testkronSubReci= " << (convertNTLzzpX2CF (TF2, x) == convertnmod_poly_t2FacCF (F2, x)) << "\n";*/
2079  int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
2080
2081
2082  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
2083
2084  //CanonicalForm NTLres= reverseSubst (TF1, TF2, d1, d2);
2085  //cout << "FINALtestkronSubReci= " << (NTLres == result) << "\n";
2086
2087  fmpz_poly_clear (F1);
2088  fmpz_poly_clear (F2);
2089  fmpz_poly_clear (G1);
2090  fmpz_poly_clear (G2);
2091  return result;
2092}
2093
2094CanonicalForm
2095mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
2096                CanonicalForm& M)
2097{
2098  CanonicalForm A= F;
2099  CanonicalForm B= G;
2100
2101  int degAx= degree (A, 1);
2102  //int degAy= degree (A, 2);
2103  int degBx= degree (B, 1);
2104  //int degBy= degree (B, 2);
2105  int d1= degAx + 1 + degBx;
2106  //int d2= tmax (degAy, degBy);
2107
2108  CanonicalForm f= bCommonDen (F);
2109  CanonicalForm g= bCommonDen (G);
2110  A *= f;
2111  B *= g;
2112
2113  fmpz_poly_t FLINTA, FLINTB;
2114  fmpz_poly_init (FLINTA);
2115  fmpz_poly_init (FLINTB);
2116  kronSub (FLINTA, A, d1);
2117  kronSub (FLINTB, B, d1);
2118  int k= d1*degree (M);
2119
2120  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
2121  A= reverseSubstQ (FLINTA, d1);
2122  fmpz_poly_clear (FLINTA);
2123  fmpz_poly_clear (FLINTB);
2124  return A/(f*g);
2125}
2126
2127#endif
2128
2129CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
2130                       const CanonicalForm& M)
2131{
2132  if (A.isZero() || B.isZero())
2133    return 0;
2134
2135  ASSERT (M.isUnivariate(), "M must be univariate");
2136
2137  CanonicalForm F= mod (A, M);
2138  CanonicalForm G= mod (B, M);
2139  if (F.inCoeffDomain() || G.inCoeffDomain())
2140    return F*G;
2141  Variable y= M.mvar();
2142  int degF= degree (F, y);
2143  int degG= degree (G, y);
2144
2145  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
2146      (F.level() == G.level()))
2147  {
2148    CanonicalForm result= mulNTL (F, G);
2149    return mod (result, M);
2150  }
2151  else if (degF <= 1 && degG <= 1)
2152  {
2153    CanonicalForm result= F*G;
2154    return mod (result, M);
2155  }
2156
2157  int sizeF= size (F);
2158  int sizeG= size (G);
2159
2160  int fallBackToNaive= 50;
2161  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
2162    return mod (F*G, M);
2163
2164#ifdef HAVE_FLINT
2165  Variable alpha;
2166  if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
2167  {
2168    CanonicalForm result= mulMod2FLINTQ (F, G, M);
2169    CanonicalForm test= mod (F*G, M);
2170    if (test != result)
2171      printf ("GOTTVERDAMMT\n");
2172    return result;
2173    //return mulMod2FLINTQ (F, G, M);
2174  }
2175#endif
2176
2177  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
2178      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
2179    return mulMod2NTLFq (F, G, M);
2180
2181  int m= (int) ceil (degree (M)/2.0);
2182  if (degF >= m || degG >= m)
2183  {
2184    CanonicalForm MLo= power (y, m);
2185    CanonicalForm MHi= power (y, degree (M) - m);
2186    CanonicalForm F0= mod (F, MLo);
2187    CanonicalForm F1= div (F, MLo);
2188    CanonicalForm G0= mod (G, MLo);
2189    CanonicalForm G1= div (G, MLo);
2190    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
2191    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
2192    CanonicalForm F0G0= mulMod2 (F0, G0, M);
2193    return F0G0 + MLo*(F0G1 + F1G0);
2194  }
2195  else
2196  {
2197    m= (int) ceil (tmax (degF, degG)/2.0);
2198    CanonicalForm yToM= power (y, m);
2199    CanonicalForm F0= mod (F, yToM);
2200    CanonicalForm F1= div (F, yToM);
2201    CanonicalForm G0= mod (G, yToM);
2202    CanonicalForm G1= div (G, yToM);
2203    CanonicalForm H00= mulMod2 (F0, G0, M);
2204    CanonicalForm H11= mulMod2 (F1, G1, M);
2205    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
2206    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
2207  }
2208  DEBOUTLN (cerr, "fatal end in mulMod2");
2209}
2210
2211CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
2212                      const CFList& MOD)
2213{
2214  if (A.isZero() || B.isZero())
2215    return 0;
2216
2217  if (MOD.length() == 1)
2218    return mulMod2 (A, B, MOD.getLast());
2219
2220  CanonicalForm M= MOD.getLast();
2221  CanonicalForm F= mod (A, M);
2222  CanonicalForm G= mod (B, M);
2223  if (F.inCoeffDomain() || G.inCoeffDomain())
2224    return F*G;
2225  Variable y= M.mvar();
2226  int degF= degree (F, y);
2227  int degG= degree (G, y);
2228
2229  if ((degF <= 1 && F.level() <= M.level()) &&
2230      (degG <= 1 && G.level() <= M.level()))
2231  {
2232    CFList buf= MOD;
2233    buf.removeLast();
2234    if (degF == 1 && degG == 1)
2235    {
2236      CanonicalForm F0= mod (F, y);
2237      CanonicalForm F1= div (F, y);
2238      CanonicalForm G0= mod (G, y);
2239      CanonicalForm G1= div (G, y);
2240      if (degree (M) > 2)
2241      {
2242        CanonicalForm H00= mulMod (F0, G0, buf);
2243        CanonicalForm H11= mulMod (F1, G1, buf);
2244        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
2245        return H11*y*y + (H01 - H00 - H11)*y + H00;
2246      }
2247      else //here degree (M) == 2
2248      {
2249        buf.append (y);
2250        CanonicalForm F0G1= mulMod (F0, G1, buf);
2251        CanonicalForm F1G0= mulMod (F1, G0, buf);
2252        CanonicalForm F0G0= mulMod (F0, G0, MOD);
2253        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
2254        return result;
2255      }
2256    }
2257    else if (degF == 1 && degG == 0)
2258      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
2259    else if (degF == 0 && degG == 1)
2260      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
2261    else
2262      return mulMod (F, G, buf);
2263  }
2264  int m= (int) ceil (degree (M)/2.0);
2265  if (degF >= m || degG >= m)
2266  {
2267    CanonicalForm MLo= power (y, m);
2268    CanonicalForm MHi= power (y, degree (M) - m);
2269    CanonicalForm F0= mod (F, MLo);
2270    CanonicalForm F1= div (F, MLo);
2271    CanonicalForm G0= mod (G, MLo);
2272    CanonicalForm G1= div (G, MLo);
2273    CFList buf= MOD;
2274    buf.removeLast();
2275    buf.append (MHi);
2276    CanonicalForm F0G1= mulMod (F0, G1, buf);
2277    CanonicalForm F1G0= mulMod (F1, G0, buf);
2278    CanonicalForm F0G0= mulMod (F0, G0, MOD);
2279    return F0G0 + MLo*(F0G1 + F1G0);
2280  }
2281  else
2282  {
2283    m= (int) ceil (tmax (degF, degG)/2.0);
2284    CanonicalForm yToM= power (y, m);
2285    CanonicalForm F0= mod (F, yToM);
2286    CanonicalForm F1= div (F, yToM);
2287    CanonicalForm G0= mod (G, yToM);
2288    CanonicalForm G1= div (G, yToM);
2289    CanonicalForm H00= mulMod (F0, G0, MOD);
2290    CanonicalForm H11= mulMod (F1, G1, MOD);
2291    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
2292    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
2293  }
2294  DEBOUTLN (cerr, "fatal end in mulMod");
2295}
2296
2297CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
2298{
2299  if (L.isEmpty())
2300    return 1;
2301  int l= L.length();
2302  if (l == 1)
2303    return mod (L.getFirst(), M);
2304  else if (l == 2) {
2305    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
2306    return result;
2307  }
2308  else
2309  {
2310    l /= 2;
2311    CFList tmp1, tmp2;
2312    CFListIterator i= L;
2313    CanonicalForm buf1, buf2;
2314    for (int j= 1; j <= l; j++, i++)
2315      tmp1.append (i.getItem());
2316    tmp2= Difference (L, tmp1);
2317    buf1= prodMod (tmp1, M);
2318    buf2= prodMod (tmp2, M);
2319    CanonicalForm result= mulMod2 (buf1, buf2, M);
2320    return result;
2321  }
2322}
2323
2324CanonicalForm prodMod (const CFList& L, const CFList& M)
2325{
2326  if (L.isEmpty())
2327    return 1;
2328  else if (L.length() == 1)
2329    return L.getFirst();
2330  else if (L.length() == 2)
2331    return mulMod (L.getFirst(), L.getLast(), M);
2332  else
2333  {
2334    int l= L.length()/2;
2335    CFListIterator i= L;
2336    CFList tmp1, tmp2;
2337    CanonicalForm buf1, buf2;
2338    for (int j= 1; j <= l; j++, i++)
2339      tmp1.append (i.getItem());
2340    tmp2= Difference (L, tmp1);
2341    buf1= prodMod (tmp1, M);
2342    buf2= prodMod (tmp2, M);
2343    return mulMod (buf1, buf2, M);
2344  }
2345}
2346
2347
2348CanonicalForm reverse (const CanonicalForm& F, int d)
2349{
2350  if (d == 0)
2351    return F;
2352  CanonicalForm A= F;
2353  Variable y= Variable (2);
2354  Variable x= Variable (1);
2355  if (degree (A, x) > 0)
2356  {
2357    A= swapvar (A, x, y);
2358    CanonicalForm result= 0;
2359    CFIterator i= A;
2360    while (d - i.exp() < 0)
2361      i++;
2362
2363    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
2364      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
2365    return result;
2366  }
2367  else
2368    return A*power (x, d);
2369}
2370
2371CanonicalForm
2372newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
2373{
2374  int l= ilog2(n);
2375
2376  CanonicalForm g= mod (F, M)[0] [0];
2377
2378  ASSERT (!g.isZero(), "expected a unit");
2379
2380  Variable alpha;
2381
2382  if (!g.isOne())
2383    g = 1/g;
2384  Variable x= Variable (1);
2385  CanonicalForm result;
2386  int exp= 0;
2387  if (n & 1)
2388  {
2389    result= g;
2390    exp= 1;
2391  }
2392  CanonicalForm h;
2393
2394  for (int i= 1; i <= l; i++)
2395  {
2396    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
2397    h= mod (h, power (x, (1 << i)) - 1);
2398    h= div (h, power (x, (1 << (i - 1))));
2399    h= mod (h, M);
2400    g -= power (x, (1 << (i - 1)))*
2401         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
2402
2403    if (n & (1 << i))
2404    {
2405      if (exp)
2406      {
2407        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
2408        h= mod (h, power (x, exp + (1 << i)) - 1);
2409        h= div (h, power (x, exp));
2410        h= mod (h, M);
2411        result -= power(x, exp)*mod (mulMod2 (g, h, M),
2412                                       power (x, (1 << i)));
2413        exp += (1 << i);
2414      }
2415      else
2416      {
2417        exp= (1 << i);
2418        result= g;
2419      }
2420    }
2421  }
2422
2423  return result;
2424}
2425
2426CanonicalForm
2427newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
2428           M)
2429{
2430  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
2431  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
2432
2433  CanonicalForm A= mod (F, M);
2434  CanonicalForm B= mod (G, M);
2435
2436  Variable x= Variable (1);
2437  int degA= degree (A, x);
2438  int degB= degree (B, x);
2439  int m= degA - degB;
2440  if (m < 0)
2441    return 0;
2442
2443  Variable v;
2444  CanonicalForm Q;
2445  if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
2446  {
2447    CanonicalForm R;
2448    divrem2 (A, B, Q, R, M);
2449  }
2450  else
2451  {
2452    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
2453    {
2454      CanonicalForm R= reverse (A, degA);
2455      CanonicalForm revB= reverse (B, degB);
2456      revB= newtonInverse (revB, m + 1, M);
2457      Q= mulMod2 (R, revB, M);
2458      Q= mod (Q, power (x, m + 1));
2459      Q= reverse (Q, m);
2460    }
2461    else
2462    {
2463      zz_pX mipo= convertFacCF2NTLzzpX (M);
2464      Variable y= Variable (2);
2465      zz_pEX NTLA, NTLB;
2466      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
2467      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
2468      div (NTLA, NTLA, NTLB);
2469      Q= convertNTLzz_pEX2CF (NTLA, x, y);
2470    }
2471  }
2472
2473  return Q;
2474}
2475
2476void
2477newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2478              CanonicalForm& R, const CanonicalForm& M)
2479{
2480  CanonicalForm A= mod (F, M);
2481  CanonicalForm B= mod (G, M);
2482  Variable x= Variable (1);
2483  int degA= degree (A, x);
2484  int degB= degree (B, x);
2485  int m= degA - degB;
2486
2487  if (m < 0)
2488  {
2489    R= A;
2490    Q= 0;
2491    return;
2492  }
2493
2494  Variable v;
2495  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
2496  {
2497     divrem2 (A, B, Q, R, M);
2498  }
2499  else
2500  {
2501    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
2502    {
2503      R= reverse (A, degA);
2504
2505      CanonicalForm revB= reverse (B, degB);
2506      revB= newtonInverse (revB, m + 1, M);
2507      Q= mulMod2 (R, revB, M);
2508
2509      Q= mod (Q, power (x, m + 1));
2510      Q= reverse (Q, m);
2511
2512      R= A - mulMod2 (Q, B, M);
2513    }
2514    else
2515    {
2516      zz_pX mipo= convertFacCF2NTLzzpX (M);
2517      Variable y= Variable (2);
2518      zz_pEX NTLA, NTLB;
2519      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
2520      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
2521      zz_pEX NTLQ, NTLR;
2522      DivRem (NTLQ, NTLR, NTLA, NTLB);
2523      Q= convertNTLzz_pEX2CF (NTLQ, x, y);
2524      R= convertNTLzz_pEX2CF (NTLR, x, y);
2525    }
2526  }
2527}
2528
2529static inline
2530CFList split (const CanonicalForm& F, const int m, const Variable& x)
2531{
2532  CanonicalForm A= F;
2533  CanonicalForm buf= 0;
2534  bool swap= false;
2535  if (degree (A, x) <= 0)
2536    return CFList(A);
2537  else if (x.level() != A.level())
2538  {
2539    swap= true;
2540    A= swapvar (A, x, A.mvar());
2541  }
2542
2543  int j= (int) floor ((double) degree (A)/ m);
2544  CFList result;
2545  CFIterator i= A;
2546  for (; j >= 0; j--)
2547  {
2548    while (i.hasTerms() && i.exp() - j*m >= 0)
2549    {
2550      if (swap)
2551        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
2552      else
2553        buf += i.coeff()*power (x, i.exp() - j*m);
2554      i++;
2555    }
2556    if (swap)
2557      result.append (swapvar (buf, x, F.mvar()));
2558    else
2559      result.append (buf);
2560    buf= 0;
2561  }
2562  return result;
2563}
2564
2565static inline
2566void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2567               CanonicalForm& R, const CFList& M);
2568
2569static inline
2570void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2571               CanonicalForm& R, const CFList& M)
2572{
2573  CanonicalForm A= mod (F, M);
2574  CanonicalForm B= mod (G, M);
2575  Variable x= Variable (1);
2576  int degB= degree (B, x);
2577  int degA= degree (A, x);
2578  if (degA < degB)
2579  {
2580    Q= 0;
2581    R= A;
2582    return;
2583  }
2584  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
2585  if (degB < 1)
2586  {
2587    divrem (A, B, Q, R);
2588    Q= mod (Q, M);
2589    R= mod (R, M);
2590    return;
2591  }
2592
2593  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
2594  CFList splitA= split (A, m, x);
2595  if (splitA.length() == 3)
2596    splitA.insert (0);
2597  if (splitA.length() == 2)
2598  {
2599    splitA.insert (0);
2600    splitA.insert (0);
2601  }
2602  if (splitA.length() == 1)
2603  {
2604    splitA.insert (0);
2605    splitA.insert (0);
2606    splitA.insert (0);
2607  }
2608
2609  CanonicalForm xToM= power (x, m);
2610
2611  CFListIterator i= splitA;
2612  CanonicalForm H= i.getItem();
2613  i++;
2614  H *= xToM;
2615  H += i.getItem();
2616  i++;
2617  H *= xToM;
2618  H += i.getItem();
2619  i++;
2620
2621  divrem32 (H, B, Q, R, M);
2622
2623  CFList splitR= split (R, m, x);
2624  if (splitR.length() == 1)
2625    splitR.insert (0);
2626
2627  H= splitR.getFirst();
2628  H *= xToM;
2629  H += splitR.getLast();
2630  H *= xToM;
2631  H += i.getItem();
2632
2633  CanonicalForm bufQ;
2634  divrem32 (H, B, bufQ, R, M);
2635
2636  Q *= xToM;
2637  Q += bufQ;
2638  return;
2639}
2640
2641static inline
2642void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2643               CanonicalForm& R, const CFList& M)
2644{
2645  CanonicalForm A= mod (F, M);
2646  CanonicalForm B= mod (G, M);
2647  Variable x= Variable (1);
2648  int degB= degree (B, x);
2649  int degA= degree (A, x);
2650  if (degA < degB)
2651  {
2652    Q= 0;
2653    R= A;
2654    return;
2655  }
2656  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
2657  if (degB < 1)
2658  {
2659    divrem (A, B, Q, R);
2660    Q= mod (Q, M);
2661    R= mod (R, M);
2662    return;
2663  }
2664  int m= (int) ceil ((double) (degB + 1)/ 2.0);
2665
2666  CFList splitA= split (A, m, x);
2667  CFList splitB= split (B, m, x);
2668
2669  if (splitA.length() == 2)
2670  {
2671    splitA.insert (0);
2672  }
2673  if (splitA.length() == 1)
2674  {
2675    splitA.insert (0);
2676    splitA.insert (0);
2677  }
2678  CanonicalForm xToM= power (x, m);
2679
2680  CanonicalForm H;
2681  CFListIterator i= splitA;
2682  i++;
2683
2684  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
2685  {
2686    H= splitA.getFirst()*xToM + i.getItem();
2687    divrem21 (H, splitB.getFirst(), Q, R, M);
2688  }
2689  else
2690  {
2691    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
2692       splitB.getFirst()*xToM;
2693    Q= xToM - 1;
2694  }
2695
2696  H= mulMod (Q, splitB.getLast(), M);
2697
2698  R= R*xToM + splitA.getLast() - H;
2699
2700  while (degree (R, x) >= degB)
2701  {
2702    xToM= power (x, degree (R, x) - degB);
2703    Q += LC (R, x)*xToM;
2704    R -= mulMod (LC (R, x), B, M)*xToM;
2705    Q= mod (Q, M);
2706    R= mod (R, M);
2707  }
2708
2709  return;
2710}
2711
2712void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2713              CanonicalForm& R, const CanonicalForm& M)
2714{
2715  CanonicalForm A= mod (F, M);
2716  CanonicalForm B= mod (G, M);
2717
2718  if (B.inCoeffDomain())
2719  {
2720    divrem (A, B, Q, R);
2721    return;
2722  }
2723  if (A.inCoeffDomain() && !B.inCoeffDomain())
2724  {
2725    Q= 0;
2726    R= A;
2727    return;
2728  }
2729
2730  if (B.level() < A.level())
2731  {
2732    divrem (A, B, Q, R);
2733    return;
2734  }
2735  if (A.level() > B.level())
2736  {
2737    R= A;
2738    Q= 0;
2739    return;
2740  }
2741  if (B.level() == 1 && B.isUnivariate())
2742  {
2743    divrem (A, B, Q, R);
2744    return;
2745  }
2746  if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
2747  {
2748    Q= 0;
2749    R= A;
2750    return;
2751  }
2752
2753  Variable x= Variable (1);
2754  int degB= degree (B, x);
2755  if (degB > degree (A, x))
2756  {
2757    Q= 0;
2758    R= A;
2759    return;
2760  }
2761
2762  CFList splitA= split (A, degB, x);
2763
2764  CanonicalForm xToDegB= power (x, degB);
2765  CanonicalForm H, bufQ;
2766  Q= 0;
2767  CFListIterator i= splitA;
2768  H= i.getItem()*xToDegB;
2769  i++;
2770  H += i.getItem();
2771  CFList buf;
2772  while (i.hasItem())
2773  {
2774    buf= CFList (M);
2775    divrem21 (H, B, bufQ, R, buf);
2776    i++;
2777    if (i.hasItem())
2778      H= R*xToDegB + i.getItem();
2779    Q *= xToDegB;
2780    Q += bufQ;
2781  }
2782  return;
2783}
2784
2785void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2786             CanonicalForm& R, const CFList& MOD)
2787{
2788  CanonicalForm A= mod (F, MOD);
2789  CanonicalForm B= mod (G, MOD);
2790  Variable x= Variable (1);
2791  int degB= degree (B, x);
2792  if (degB > degree (A, x))
2793  {
2794    Q= 0;
2795    R= A;
2796    return;
2797  }
2798
2799  if (degB <= 0)
2800  {
2801    divrem (A, B, Q, R);
2802    Q= mod (Q, MOD);
2803    R= mod (R, MOD);
2804    return;
2805  }
2806  CFList splitA= split (A, degB, x);
2807
2808  CanonicalForm xToDegB= power (x, degB);
2809  CanonicalForm H, bufQ;
2810  Q= 0;
2811  CFListIterator i= splitA;
2812  H= i.getItem()*xToDegB;
2813  i++;
2814  H += i.getItem();
2815  while (i.hasItem())
2816  {
2817    divrem21 (H, B, bufQ, R, MOD);
2818    i++;
2819    if (i.hasItem())
2820      H= R*xToDegB + i.getItem();
2821    Q *= xToDegB;
2822    Q += bufQ;
2823  }
2824  return;
2825}
2826
2827void sortList (CFList& list, const Variable& x)
2828{
2829  int l= 1;
2830  int k= 1;
2831  CanonicalForm buf;
2832  CFListIterator m;
2833  for (CFListIterator i= list; l <= list.length(); i++, l++)
2834  {
2835    for (CFListIterator j= list; k <= list.length() - l; k++)
2836    {
2837      m= j;
2838      m++;
2839      if (degree (j.getItem(), x) > degree (m.getItem(), x))
2840      {
2841        buf= m.getItem();
2842        m.getItem()= j.getItem();
2843        j.getItem()= buf;
2844        j++;
2845        j.getItem()= m.getItem();
2846      }
2847      else
2848        j++;
2849    }
2850    k= 1;
2851  }
2852}
2853
2854static inline
2855CFList diophantine (const CanonicalForm& F, const CFList& factors)
2856{
2857  if (getCharacteristic() == 0)
2858  {
2859    Variable v;
2860    bool hasAlgVar= hasFirstAlgVar (F, v);
2861    for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++)
2862      hasAlgVar= hasFirstAlgVar (i.getItem(), v);
2863    if (hasAlgVar)
2864    {
2865      CFList result= modularDiophant (F, factors, getMipo (v));
2866      return result;
2867    }
2868  }
2869
2870  CanonicalForm buf1, buf2, buf3, S, T;
2871  CFListIterator i= factors;
2872  CFList result;
2873  if (i.hasItem())
2874    i++;
2875  buf1= F/factors.getFirst();
2876  buf2= divNTL (F, i.getItem());
2877  buf3= extgcd (buf1, buf2, S, T);
2878  result.append (S);
2879  result.append (T);
2880  if (i.hasItem())
2881    i++;
2882  for (; i.hasItem(); i++)
2883  {
2884    buf1= divNTL (F, i.getItem());
2885    buf3= extgcd (buf3, buf1, S, T);
2886    CFListIterator k= factors;
2887    for (CFListIterator j= result; j.hasItem(); j++, k++)
2888    {
2889      j.getItem()= mulNTL (j.getItem(), S);
2890      j.getItem()= modNTL (j.getItem(), k.getItem());
2891    }
2892    result.append (T);
2893  }
2894  return result;
2895}
2896
2897void
2898henselStep12 (const CanonicalForm& F, const CFList& factors,
2899              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
2900              CFArray& Pi, int j)
2901{
2902  CanonicalForm E;
2903  CanonicalForm xToJ= power (F.mvar(), j);
2904  Variable x= F.mvar();
2905  // compute the error
2906  if (j == 1)
2907    E= F[j];
2908  else
2909  {
2910    if (degree (Pi [factors.length() - 2], x) > 0)
2911      E= F[j] - Pi [factors.length() - 2] [j];
2912    else
2913      E= F[j];
2914  }
2915
2916  CFArray buf= CFArray (diophant.length());
2917  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
2918  int k= 0;
2919  CanonicalForm remainder;
2920  // actual lifting
2921  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
2922  {
2923    if (degree (bufFactors[k], x) > 0)
2924    {
2925      if (k > 0)
2926        remainder= modNTL (E, bufFactors[k] [0]);
2927      else
2928        remainder= E;
2929    }
2930    else
2931      remainder= modNTL (E, bufFactors[k]);
2932
2933    buf[k]= mulNTL (i.getItem(), remainder);
2934    if (degree (bufFactors[k], x) > 0)
2935      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
2936    else
2937      buf[k]= modNTL (buf[k], bufFactors[k]);
2938  }
2939  for (k= 1; k < factors.length(); k++)
2940    bufFactors[k] += xToJ*buf[k];
2941
2942  // update Pi [0]
2943  int degBuf0= degree (bufFactors[0], x);
2944  int degBuf1= degree (bufFactors[1], x);
2945  if (degBuf0 > 0 && degBuf1 > 0)
2946    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
2947  CanonicalForm uIZeroJ;
2948  if (j == 1)
2949  {
2950    if (degBuf0 > 0 && degBuf1 > 0)
2951      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
2952                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
2953    else if (degBuf0 > 0)
2954      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
2955    else if (degBuf1 > 0)
2956      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
2957    else
2958      uIZeroJ= 0;
2959    Pi [0] += xToJ*uIZeroJ;
2960  }
2961  else
2962  {
2963    if (degBuf0 > 0 && degBuf1 > 0)
2964      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
2965                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
2966    else if (degBuf0 > 0)
2967      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
2968    else if (degBuf1 > 0)
2969      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
2970    else
2971      uIZeroJ= 0;
2972    Pi [0] += xToJ*uIZeroJ;
2973  }
2974  CFArray tmp= CFArray (factors.length() - 1);
2975  for (k= 0; k < factors.length() - 1; k++)
2976    tmp[k]= 0;
2977  CFIterator one, two;
2978  one= bufFactors [0];
2979  two= bufFactors [1];
2980  if (degBuf0 > 0 && degBuf1 > 0)
2981  {
2982    for (k= 1; k <= (int) ceil (j/2.0); k++)
2983    {
2984      if (k != j - k + 1)
2985      {
2986        if ((one.hasTerms() && one.exp() == j - k + 1) && (two.hasTerms() && two.exp() == j - k + 1))
2987        {
2988          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
2989                     two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
2990          one++;
2991          two++;
2992        }
2993        else if (one.hasTerms() && one.exp() == j - k + 1)
2994        {
2995          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) -
2996                     M (k + 1, 1);
2997          one++;
2998        }
2999        else if (two.hasTerms() && two.exp() == j - k + 1)
3000        {
3001          tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) -
3002                    M (k + 1, 1);
3003          two++;
3004        }
3005      }
3006      else
3007      {
3008        tmp[0] += M (k + 1, 1);
3009      }
3010    }
3011  }
3012  Pi [0] += tmp[0]*xToJ*F.mvar();
3013
3014  // update Pi [l]
3015  int degPi, degBuf;
3016  for (int l= 1; l < factors.length() - 1; l++)
3017  {
3018    degPi= degree (Pi [l - 1], x);
3019    degBuf= degree (bufFactors[l + 1], x);
3020    if (degPi > 0 && degBuf > 0)
3021      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
3022    if (j == 1)
3023    {
3024      if (degPi > 0 && degBuf > 0)
3025        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
3026                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
3027                  M (1, l + 1));
3028      else if (degPi > 0)
3029        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
3030      else if (degBuf > 0)
3031        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
3032    }
3033    else
3034    {
3035      if (degPi > 0 && degBuf > 0)
3036      {
3037        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
3038        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
3039      }
3040      else if (degPi > 0)
3041        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
3042      else if (degBuf > 0)
3043      {
3044        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
3045        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
3046      }
3047      Pi[l] += xToJ*uIZeroJ;
3048    }
3049    one= bufFactors [l + 1];
3050    two= Pi [l - 1];
3051    if (two.hasTerms() && two.exp() == j + 1)
3052    {
3053      if (degBuf > 0 && degPi > 0)
3054      {
3055          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
3056          two++;
3057      }
3058      else if (degPi > 0)
3059      {
3060          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
3061          two++;
3062      }
3063    }
3064    if (degBuf > 0 && degPi > 0)
3065    {
3066      for (k= 1; k <= (int) ceil (j/2.0); k++)
3067      {
3068        if (k != j - k + 1)
3069        {
3070          if ((one.hasTerms() && one.exp() == j - k + 1) &&
3071              (two.hasTerms() && two.exp() == j - k + 1))
3072          {
3073            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
3074                       two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
3075            one++;
3076            two++;
3077          }
3078          else if (one.hasTerms() && one.exp() == j - k + 1)
3079          {
3080            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) -
3081                       M (k + 1, l + 1);
3082            one++;
3083          }
3084          else if (two.hasTerms() && two.exp() == j - k + 1)
3085          {
3086            tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) -
3087                      M (k + 1, l + 1);
3088            two++;
3089          }
3090        }
3091        else
3092          tmp[l] += M (k + 1, l + 1);
3093      }
3094    }
3095    Pi[l] += tmp[l]*xToJ*F.mvar();
3096  }
3097  return;
3098}
3099
3100void
3101henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
3102              CFList& diophant, CFMatrix& M, bool sort)
3103{
3104  if (sort)
3105    sortList (factors, Variable (1));
3106  Pi= CFArray (factors.length() - 1);
3107  CFListIterator j= factors;
3108  diophant= diophantine (F[0], factors);
3109  DEBOUTLN (cerr, "diophant= " << diophant);
3110  j++;
3111  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()));
3112  M (1, 1)= Pi [0];
3113  int i= 1;
3114  if (j.hasItem())
3115    j++;
3116  for (; j.hasItem(); j++, i++)
3117  {
3118    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
3119    M (1, i + 1)= Pi [i];
3120  }
3121  CFArray bufFactors= CFArray (factors.length());
3122  i= 0;
3123  for (CFListIterator k= factors; k.hasItem(); i++, k++)
3124  {
3125    if (i == 0)
3126      bufFactors[i]= mod (k.getItem(), F.mvar());
3127    else
3128      bufFactors[i]= k.getItem();
3129  }
3130  for (i= 1; i < l; i++)
3131    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
3132
3133  CFListIterator k= factors;
3134  for (i= 0; i < factors.length (); i++, k++)
3135    k.getItem()= bufFactors[i];
3136  factors.removeFirst();
3137  return;
3138}
3139
3140void
3141henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
3142                    end, CFArray& Pi, const CFList& diophant, CFMatrix& M)
3143{
3144  CFArray bufFactors= CFArray (factors.length());
3145  int i= 0;
3146  CanonicalForm xToStart= power (F.mvar(), start);
3147  for (CFListIterator k= factors; k.hasItem(); k++, i++)
3148  {
3149    if (i == 0)
3150      bufFactors[i]= mod (k.getItem(), xToStart);
3151    else
3152      bufFactors[i]= k.getItem();
3153  }
3154  for (i= start; i < end; i++)
3155    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
3156
3157  CFListIterator k= factors;
3158  for (i= 0; i < factors.length(); k++, i++)
3159    k.getItem()= bufFactors [i];
3160  factors.removeFirst();
3161  return;
3162}
3163
3164static inline
3165CFList
3166biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
3167{
3168  Variable y= F.mvar();
3169  CFList result;
3170  if (y.level() == 1)
3171  {
3172    result= diophantine (F, factors);
3173    return result;
3174  }
3175  else
3176  {
3177    CFList buf= factors;
3178    for (CFListIterator i= buf; i.hasItem(); i++)
3179      i.getItem()= mod (i.getItem(), y);
3180    CanonicalForm A= mod (F, y);
3181    int bufD= 1;
3182    CFList recResult= biDiophantine (A, buf, bufD);
3183    CanonicalForm e= 1;
3184    CFList p;
3185    CFArray bufFactors= CFArray (factors.length());
3186    CanonicalForm yToD= power (y, d);
3187    int k= 0;
3188    for (CFListIterator i= factors; i.hasItem(); i++, k++)
3189    {
3190      bufFactors [k]= i.getItem();
3191    }
3192    CanonicalForm b, quot;
3193    for (k= 0; k < factors.length(); k++) //TODO compute b's faster
3194    {
3195      b= 1;
3196      if (fdivides (bufFactors[k], F, quot))
3197        b= quot;
3198      else
3199      {
3200        for (int l= 0; l < factors.length(); l++)
3201        {
3202          if (l == k)
3203            continue;
3204          else
3205          {
3206            b= mulMod2 (b, bufFactors[l], yToD);
3207          }
3208        }
3209      }
3210      p.append (b);
3211    }
3212
3213    CFListIterator j= p;
3214    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
3215      e -= i.getItem()*j.getItem();
3216
3217    if (e.isZero())
3218      return recResult;
3219    CanonicalForm coeffE;
3220    CFList s;
3221    result= recResult;
3222    CanonicalForm g;
3223    for (int i= 1; i < d; i++)
3224    {
3225      if (degree (e, y) > 0)
3226        coeffE= e[i];
3227      else
3228        coeffE= 0;
3229      if (!coeffE.isZero())
3230      {
3231        CFListIterator k= result;
3232        CFListIterator l= p;
3233        int ii= 0;
3234        j= recResult;
3235        for (; j.hasItem(); j++, k++, l++, ii++)
3236        {
3237          g= coeffE*j.getItem();
3238          if (degree (bufFactors[ii], y) <= 0)
3239            g= mod (g, bufFactors[ii]);
3240          else
3241            g= mod (g, bufFactors[ii][0]);
3242          k.getItem() += g*power (y, i);
3243          e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
3244          DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
3245                    mod (e, power (y, i + 1)));
3246        }
3247      }
3248      if (e.isZero())
3249        break;
3250    }
3251
3252    DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
3253
3254#ifdef DEBUGOUTPUT
3255    CanonicalForm test= 0;
3256    j= p;
3257    for (CFListIterator i= result; i.hasItem(); i++, j++)
3258      test += mod (i.getItem()*j.getItem(), power (y, d));
3259    DEBOUTLN (cerr, "test= " << test);
3260#endif
3261    return result;
3262  }
3263}
3264
3265static inline
3266CFList
3267multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
3268                     const CFList& recResult, const CFList& M, const int d)
3269{
3270  Variable y= F.mvar();
3271  CFList result;
3272  CFListIterator i;
3273  CanonicalForm e= 1;
3274  CFListIterator j= factors;
3275  CFList p;
3276  CFArray bufFactors= CFArray (factors.length());
3277  CanonicalForm yToD= power (y, d);
3278  int k= 0;
3279  for (CFListIterator i= factors; i.hasItem(); i++, k++)
3280    bufFactors [k]= i.getItem();
3281  CanonicalForm b, quot;
3282  CFList buf= M;
3283  buf.removeLast();
3284  buf.append (yToD);
3285  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
3286  {
3287    b= 1;
3288    if (fdivides (bufFactors[k], F, quot))
3289      b= quot;
3290    else
3291    {
3292      for (int l= 0; l < factors.length(); l++)
3293      {
3294        if (l == k)
3295          continue;
3296        else
3297        {
3298          b= mulMod (b, bufFactors[l], buf);
3299        }
3300      }
3301    }
3302    p.append (b);
3303  }
3304  j= p;
3305  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
3306    e -= mulMod (i.getItem(), j.getItem(), M);
3307
3308  if (e.isZero())
3309    return recResult;
3310  CanonicalForm coeffE;
3311  CFList s;
3312  result= recResult;
3313  CanonicalForm g;
3314  for (int i= 1; i < d; i++)
3315  {
3316    if (degree (e, y) > 0)
3317      coeffE= e[i];
3318    else
3319      coeffE= 0;
3320    if (!coeffE.isZero())
3321    {
3322      CFListIterator k= result;
3323      CFListIterator l= p;
3324      j= recResult;
3325      int ii= 0;
3326      CanonicalForm dummy;
3327      for (; j.hasItem(); j++, k++, l++, ii++)
3328      {
3329        g= mulMod (coeffE, j.getItem(), M);
3330        if (degree (bufFactors[ii], y) <= 0)
3331          divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
3332                  g, M);
3333        else
3334          divrem (g, bufFactors[ii][0], dummy, g, M);
3335        k.getItem() += g*power (y, i);
3336        e -= mulMod (g*power (y, i), l.getItem(), M);
3337      }
3338    }
3339
3340    if (e.isZero())
3341      break;
3342  }
3343
3344#ifdef DEBUGOUTPUT
3345    CanonicalForm test= 0;
3346    j= p;
3347    for (CFListIterator i= result; i.hasItem(); i++, j++)
3348      test += mod (i.getItem()*j.getItem(), power (y, d));
3349    DEBOUTLN (cerr, "test= " << test);
3350#endif
3351  return result;
3352}
3353
3354static inline
3355void
3356henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
3357            const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
3358            const CFList& MOD)
3359{
3360  CanonicalForm E;
3361  CanonicalForm xToJ= power (F.mvar(), j);
3362  Variable x= F.mvar();
3363  // compute the error
3364  if (j == 1)
3365  {
3366    E= F[j];
3367#ifdef DEBUGOUTPUT
3368    CanonicalForm test= 1;
3369    for (int i= 0; i < factors.length(); i++)
3370    {
3371      if (i == 0)
3372        test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
3373      else
3374        test= mulMod (test, bufFactors[i], MOD);
3375    }
3376    CanonicalForm test2= mod (F-test, xToJ);
3377
3378    test2= mod (test2, MOD);
3379    DEBOUTLN (cerr, "test= " << test2);
3380#endif
3381  }
3382  else
3383  {
3384#ifdef DEBUGOUTPUT
3385    CanonicalForm test= 1;
3386    for (int i= 0; i < factors.length(); i++)
3387    {
3388      if (i == 0)
3389        test *= mod (bufFactors [i], power (x, j));
3390      else
3391        test *= bufFactors[i];
3392    }
3393    test= mod (test, power (x, j));
3394    test= mod (test, MOD);
3395    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
3396    DEBOUTLN (cerr, "test= " << test2);
3397#endif
3398
3399    if (degree (Pi [factors.length() - 2], x) > 0)
3400      E= F[j] - Pi [factors.length() - 2] [j];
3401    else
3402      E= F[j];
3403  }
3404
3405  CFArray buf= CFArray (diophant.length());
3406  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
3407  int k= 0;
3408  // actual lifting
3409  CanonicalForm dummy, rest1;
3410  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
3411  {
3412    if (degree (bufFactors[k], x) > 0)
3413    {
3414      if (k > 0)
3415        divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
3416      else
3417        rest1= E;
3418    }
3419    else
3420      divrem (E, bufFactors[k], dummy, rest1, MOD);
3421
3422    buf[k]= mulMod (i.getItem(), rest1, MOD);
3423
3424    if (degree (bufFactors[k], x) > 0)
3425      divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
3426    else
3427      divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
3428  }
3429  for (k= 1; k < factors.length(); k++)
3430    bufFactors[k] += xToJ*buf[k];
3431
3432  // update Pi [0]
3433  int degBuf0= degree (bufFactors[0], x);
3434  int degBuf1= degree (bufFactors[1], x);
3435  if (degBuf0 > 0 && degBuf1 > 0)
3436    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
3437  CanonicalForm uIZeroJ;
3438  if (j == 1)
3439  {
3440    if (degBuf0 > 0 && degBuf1 > 0)
3441      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
3442                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
3443    else if (degBuf0 > 0)
3444      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
3445    else if (degBuf1 > 0)
3446      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
3447    else
3448      uIZeroJ= 0;
3449    Pi [0] += xToJ*uIZeroJ;
3450  }
3451  else
3452  {
3453    if (degBuf0 > 0 && degBuf1 > 0)
3454      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
3455                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
3456    else if (degBuf0 > 0)
3457      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
3458    else if (degBuf1 > 0)
3459      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
3460    else
3461      uIZeroJ= 0;
3462    Pi [0] += xToJ*uIZeroJ;
3463  }
3464
3465  CFArray tmp= CFArray (factors.length() - 1);
3466  for (k= 0; k < factors.length() - 1; k++)
3467    tmp[k]= 0;
3468  CFIterator one, two;
3469  one= bufFactors [0];
3470  two= bufFactors [1];
3471  if (degBuf0 > 0 && degBuf1 > 0)
3472  {
3473    for (k= 1; k <= (int) ceil (j/2.0); k++)
3474    {
3475      if (k != j - k + 1)
3476      {
3477        if ((one.hasTerms() && one.exp() == j - k + 1) &&
3478            (two.hasTerms() && two.exp() == j - k + 1))
3479        {
3480          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
3481                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
3482                    M (j - k + 2, 1);
3483          one++;
3484          two++;
3485        }
3486        else if (one.hasTerms() && one.exp() == j - k + 1)
3487        {
3488          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
3489                    bufFactors[1] [k], MOD) - M (k + 1, 1);
3490          one++;
3491        }
3492        else if (two.hasTerms() && two.exp() == j - k + 1)
3493        {
3494          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
3495                    two.coeff()), MOD) - M (k + 1, 1);
3496          two++;
3497        }
3498      }
3499      else
3500      {
3501        tmp[0] += M (k + 1, 1);
3502      }
3503    }
3504  }
3505  Pi [0] += tmp[0]*xToJ*F.mvar();
3506
3507  // update Pi [l]
3508  int degPi, degBuf;
3509  for (int l= 1; l < factors.length() - 1; l++)
3510  {
3511    degPi= degree (Pi [l - 1], x);
3512    degBuf= degree (bufFactors[l + 1], x);
3513    if (degPi > 0 && degBuf > 0)
3514      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
3515    if (j == 1)
3516    {
3517      if (degPi > 0 && degBuf > 0)
3518        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
3519                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
3520                  M (1, l + 1));
3521      else if (degPi > 0)
3522        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
3523      else if (degBuf > 0)
3524        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
3525    }
3526    else
3527    {
3528      if (degPi > 0 && degBuf > 0)
3529      {
3530        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
3531        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
3532      }
3533      else if (degPi > 0)
3534        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
3535      else if (degBuf > 0)
3536      {
3537        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
3538        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
3539      }
3540      Pi[l] += xToJ*uIZeroJ;
3541    }
3542    one= bufFactors [l + 1];
3543    two= Pi [l - 1];
3544    if (two.hasTerms() && two.exp() == j + 1)
3545    {
3546      if (degBuf > 0 && degPi > 0)
3547      {
3548          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
3549          two++;
3550      }
3551      else if (degPi > 0)
3552      {
3553          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
3554          two++;
3555      }
3556    }
3557    if (degBuf > 0 && degPi > 0)
3558    {
3559      for (k= 1; k <= (int) ceil (j/2.0); k++)
3560      {
3561        if (k != j - k + 1)
3562        {
3563          if ((one.hasTerms() && one.exp() == j - k + 1) &&
3564              (two.hasTerms() && two.exp() == j - k + 1))
3565          {
3566            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
3567                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
3568                      M (j - k + 2, l + 1);
3569            one++;
3570            two++;
3571          }
3572          else if (one.hasTerms() && one.exp() == j - k + 1)
3573          {
3574            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
3575                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
3576            one++;
3577          }
3578          else if (two.hasTerms() && two.exp() == j - k + 1)
3579          {
3580            tmp[l] += mulMod (bufFactors[l + 1] [k],
3581                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
3582            two++;
3583          }
3584        }
3585        else
3586          tmp[l] += M (k + 1, l + 1);
3587      }
3588    }
3589    Pi[l] += tmp[l]*xToJ*F.mvar();
3590  }
3591
3592  return;
3593}
3594
3595CFList
3596henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
3597              diophant, CFArray& Pi, CFMatrix& M)
3598{
3599  CFList buf= factors;
3600  int k= 0;
3601  int liftBoundBivar= l[k];
3602  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
3603  CFList MOD;
3604  MOD.append (power (Variable (2), liftBoundBivar));
3605  CFArray bufFactors= CFArray (factors.length());
3606  k= 0;
3607  CFListIterator j= eval;
3608  j++;
3609  buf.removeFirst();
3610  buf.insert (LC (j.getItem(), 1));
3611  for (CFListIterator i= buf; i.hasItem(); i++, k++)
3612    bufFactors[k]= i.getItem();
3613  Pi= CFArray (factors.length() - 1);
3614  CFListIterator i= buf;
3615  i++;
3616  Variable y= j.getItem().mvar();
3617  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
3618  M (1, 1)= Pi [0];
3619  k= 1;
3620  if (i.hasItem())
3621    i++;
3622  for (; i.hasItem(); i++, k++)
3623  {
3624    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
3625    M (1, k + 1)= Pi [k];
3626  }
3627
3628  for (int d= 1; d < l[1]; d++)
3629    henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
3630  CFList result;
3631  for (k= 1; k < factors.length(); k++)
3632    result.append (bufFactors[k]);
3633  return result;
3634}
3635
3636void
3637henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
3638                  CFArray& Pi, const CFList& diophant, CFMatrix& M,
3639                  const CFList& MOD)
3640{
3641  CFArray bufFactors= CFArray (factors.length());
3642  int i= 0;
3643  CanonicalForm xToStart= power (F.mvar(), start);
3644  for (CFListIterator k= factors; k.hasItem(); k++, i++)
3645  {
3646    if (i == 0)
3647      bufFactors[i]= mod (k.getItem(), xToStart);
3648    else
3649      bufFactors[i]= k.getItem();
3650  }
3651  for (i= start; i < end; i++)
3652    henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD);
3653
3654  CFListIterator k= factors;
3655  for (i= 0; i < factors.length(); k++, i++)
3656    k.getItem()= bufFactors [i];
3657  factors.removeFirst();
3658  return;
3659}
3660
3661CFList
3662henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
3663            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
3664            lNew)
3665{
3666  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
3667  int k= 0;
3668  CFArray bufFactors= CFArray (factors.length());
3669  for (CFListIterator i= factors; i.hasItem(); i++, k++)
3670  {
3671    if (k == 0)
3672      bufFactors[k]= LC (F.getLast(), 1);
3673    else
3674      bufFactors[k]= i.getItem();
3675  }
3676  CFList buf= factors;
3677  buf.removeFirst();
3678  buf.insert (LC (F.getLast(), 1));
3679  CFListIterator i= buf;
3680  i++;
3681  Variable y= F.getLast().mvar();
3682  Variable x= F.getFirst().mvar();
3683  CanonicalForm xToLOld= power (x, lOld);
3684  Pi [0]= mod (Pi[0], xToLOld);
3685  M (1, 1)= Pi [0];
3686  k= 1;
3687  if (i.hasItem())
3688    i++;
3689  for (; i.hasItem(); i++, k++)
3690  {
3691    Pi [k]= mod (Pi [k], xToLOld);
3692    M (1, k + 1)= Pi [k];
3693  }
3694
3695  for (int d= 1; d < lNew; d++)
3696    henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
3697  CFList result;
3698  for (k= 1; k < factors.length(); k++)
3699    result.append (bufFactors[k]);
3700  return result;
3701}
3702
3703CFList
3704henselLift (const CFList& eval, const CFList& factors, const int* l, const int
3705            lLength, bool sort)
3706{
3707  CFList diophant;
3708  CFList buf= factors;
3709  buf.insert (LC (eval.getFirst(), 1));
3710  if (sort)
3711    sortList (buf, Variable (1));
3712  CFArray Pi;
3713  CFMatrix M= CFMatrix (l[1], factors.length());
3714  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
3715  if (eval.length() == 2)
3716    return result;
3717  CFList MOD;
3718  for (int i= 0; i < 2; i++)
3719    MOD.append (power (Variable (i + 2), l[i]));
3720  CFListIterator j= eval;
3721  j++;
3722  CFList bufEval;
3723  bufEval.append (j.getItem());
3724  j++;
3725
3726  for (int i= 2; i < lLength && j.hasItem(); i++, j++)
3727  {
3728    result.insert (LC (bufEval.getFirst(), 1));
3729    bufEval.append (j.getItem());
3730    M= CFMatrix (l[i], factors.length());
3731    result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
3732    MOD.append (power (Variable (i + 2), l[i]));
3733    bufEval.removeFirst();
3734  }
3735  return result;
3736}
3737
3738void
3739henselStep122 (const CanonicalForm& F, const CFList& factors,
3740              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
3741              CFArray& Pi, int j, const CFArray& /*LCs*/)
3742{
3743  Variable x= F.mvar();
3744  CanonicalForm xToJ= power (x, j);
3745
3746  CanonicalForm E;
3747  // compute the error
3748  if (degree (Pi [factors.length() - 2], x) > 0)
3749    E= F[j] - Pi [factors.length() - 2] [j];
3750  else
3751    E= F[j];
3752
3753  CFArray buf= CFArray (diophant.length());
3754
3755  int k= 0;
3756  CanonicalForm remainder;
3757  // actual lifting
3758  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
3759  {
3760    if (degree (bufFactors[k], x) > 0)
3761      remainder= modNTL (E, bufFactors[k] [0]);
3762    else
3763      remainder= modNTL (E, bufFactors[k]);
3764    buf[k]= mulNTL (i.getItem(), remainder);
3765    if (degree (bufFactors[k], x) > 0)
3766      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
3767    else
3768      buf[k]= modNTL (buf[k], bufFactors[k]);
3769  }
3770
3771  for (k= 0; k < factors.length(); k++)
3772    bufFactors[k] += xToJ*buf[k];
3773
3774  // update Pi [0]
3775  int degBuf0= degree (bufFactors[0], x);
3776  int degBuf1= degree (bufFactors[1], x);
3777  if (degBuf0 > 0 && degBuf1 > 0)
3778  {
3779    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
3780    if (j + 2 <= M.rows())
3781      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
3782  }
3783  CanonicalForm uIZeroJ;
3784  if (degBuf0 > 0 && degBuf1 > 0)
3785    uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]);
3786  else if (degBuf0 > 0)
3787    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
3788  else if (degBuf1 > 0)
3789    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
3790  else
3791    uIZeroJ= 0;
3792  Pi [0] += xToJ*uIZeroJ;
3793
3794  CFArray tmp= CFArray (factors.length() - 1);
3795  for (k= 0; k < factors.length() - 1; k++)
3796    tmp[k]= 0;
3797  CFIterator one, two;
3798  one= bufFactors [0];
3799  two= bufFactors [1];
3800  if (degBuf0 > 0 && degBuf1 > 0)
3801  {
3802    while (one.hasTerms() && one.exp() > j) one++;
3803    while (two.hasTerms() && two.exp() > j) two++;
3804    for (k= 1; k <= (int) ceil (j/2.0); k++)
3805    {
3806      if (one.hasTerms() && two.hasTerms())
3807      {
3808        if (k != j - k + 1)
3809        {
3810          if ((one.hasTerms() && one.exp() == j - k + 1) && +
3811              (two.hasTerms() && two.exp() == j - k + 1))
3812        {
3813          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
3814                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
3815          one++;
3816          two++;
3817        }
3818        else if (one.hasTerms() && one.exp() == j - k + 1)
3819        {
3820          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
3821                      M (k + 1, 1);
3822          one++;
3823        }
3824        else if (two.hasTerms() && two.exp() == j - k + 1)
3825        {
3826          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
3827                    M (k + 1, 1);
3828          two++;
3829        }
3830      }
3831      else
3832        tmp[0] += M (k + 1, 1);
3833      }
3834    }
3835  }
3836
3837  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
3838  {
3839    if (j + 2 <= M.rows())
3840      tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
3841                         (bufFactors [1] [j + 1] + bufFactors [1] [0]))
3842                         - M(1,1) - M (j + 2,1);
3843  }
3844  else if (degBuf0 >= j + 1)
3845  {
3846    if (degBuf1 > 0)
3847      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
3848    else
3849      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]);
3850  }
3851  else if (degBuf1 >= j + 1)
3852  {
3853    if (degBuf0 > 0)
3854      tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
3855    else
3856      tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]);
3857  }
3858
3859  Pi [0] += tmp[0]*xToJ*F.mvar();
3860
3861  int degPi, degBuf;
3862  for (int l= 1; l < factors.length() - 1; l++)
3863  {
3864    degPi= degree (Pi [l - 1], x);
3865    degBuf= degree (bufFactors[l + 1], x);
3866    if (degPi > 0 && degBuf > 0)
3867    {
3868      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
3869      if (j + 2 <= M.rows())
3870        M (j + 2, l + 1)= mulNTL (Pi [l - 1][j + 1], bufFactors[l + 1] [j + 1]);
3871    }
3872
3873    if (degPi > 0 && degBuf > 0)
3874      uIZeroJ= mulNTL (Pi[l -1] [0], buf[l + 1]) +
3875               mulNTL (uIZeroJ, bufFactors[l+1] [0]);
3876    else if (degPi > 0)
3877      uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]);
3878    else if (degBuf > 0)
3879      uIZeroJ= mulNTL (Pi[l - 1], buf[1]);
3880    else
3881      uIZeroJ= 0;
3882
3883    Pi [l] += xToJ*uIZeroJ;
3884
3885    one= bufFactors [l + 1];
3886    two= Pi [l - 1];
3887    if (degBuf > 0 && degPi > 0)
3888    {
3889      while (one.hasTerms() && one.exp() > j) one++;
3890      while (two.hasTerms() && two.exp() > j) two++;
3891      for (k= 1; k <= (int) ceil (j/2.0); k++)
3892      {
3893        if (k != j - k + 1)
3894        {
3895          if ((one.hasTerms() && one.exp() == j - k + 1) &&
3896              (two.hasTerms() && two.exp() == j - k + 1))
3897          {
3898            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
3899                      (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1) -
3900                      M (j - k + 2, l + 1);
3901            one++;
3902            two++;
3903          }
3904          else if (one.hasTerms() && one.exp() == j - k + 1)
3905          {
3906            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
3907                               Pi[l - 1] [k]) - M (k + 1, l + 1);
3908            one++;
3909          }
3910          else if (two.hasTerms() && two.exp() == j - k + 1)
3911          {
3912            tmp[l] += mulNTL (bufFactors[l + 1] [k],
3913                      (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1);
3914            two++;
3915           }
3916        }
3917        else
3918          tmp[l] += M (k + 1, l + 1);
3919      }
3920    }
3921
3922    if (degPi >= j + 1 && degBuf >= j + 1)
3923    {
3924      if (j + 2 <= M.rows())
3925        tmp [l] += mulNTL ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
3926                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
3927                          ) - M(1,l+1) - M (j + 2,l+1);
3928    }
3929    else if (degPi >= j + 1)
3930    {
3931      if (degBuf > 0)
3932        tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1] [0]);
3933      else
3934        tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1]);
3935    }
3936    else if (degBuf >= j + 1)
3937    {
3938      if (degPi > 0)
3939        tmp[l] += mulNTL (Pi [l - 1] [0], bufFactors [l + 1] [j + 1]);
3940      else
3941        tmp[l] += mulNTL (Pi [l - 1], bufFactors [l + 1] [j + 1]);
3942    }
3943
3944    Pi[l] += tmp[l]*xToJ*F.mvar();
3945  }
3946  return;
3947}
3948
3949void
3950henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
3951              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
3952{
3953  if (sort)
3954    sortList (factors, Variable (1));
3955  Pi= CFArray (factors.length() - 2);
3956  CFList bufFactors2= factors;
3957  bufFactors2.removeFirst();
3958  diophant= diophantine (F[0], bufFactors2);
3959  DEBOUTLN (cerr, "diophant= " << diophant);
3960
3961  CFArray bufFactors= CFArray (bufFactors2.length());
3962  int i= 0;
3963  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
3964    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
3965
3966  Variable x= F.mvar();
3967  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
3968  {
3969    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
3970    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
3971                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
3972  }
3973  else if (degree (bufFactors[0], x) > 0)
3974  {
3975    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
3976    Pi [0]= M (1, 1) +
3977            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
3978  }
3979  else if (degree (bufFactors[1], x) > 0)
3980  {
3981    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
3982    Pi [0]= M (1, 1) +
3983            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
3984  }
3985  else
3986  {
3987    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
3988    Pi [0]= M (1, 1);
3989  }
3990
3991  for (i= 1; i < Pi.size(); i++)
3992  {
3993    if (degree (Pi[i-1], x) > 0 && degree (bufFactors [i+1], x) > 0)
3994    {
3995      M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors[i+1] [0]);
3996      Pi [i]= M (1,i+1) + (mulNTL (Pi[i-1] [1], bufFactors[i+1] [0]) +
3997                       mulNTL (Pi[i-1] [0], bufFactors [i+1] [1]))*x;
3998    }
3999    else if (degree (Pi[i-1], x) > 0)
4000    {
4001      M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors [i+1]);
4002      Pi [i]=  M(1,i+1) + mulNTL (Pi[i-1] [1], bufFactors[i+1])*x;
4003    }
4004    else if (degree (bufFactors[i+1], x) > 0)
4005    {
4006      M (1,i+1)= mulNTL (Pi[i-1], bufFactors [i+1] [0]);
4007      Pi [i]= M (1,i+1) + mulNTL (Pi[i-1], bufFactors[i+1] [1])*x;
4008    }
4009    else
4010    {
4011      M (1,i+1)= mulNTL (Pi [i-1], bufFactors [i+1]);
4012      Pi [i]= M (1,i+1);
4013    }
4014  }
4015
4016  for (i= 1; i < l; i++)
4017    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
4018
4019  factors= CFList();
4020  for (i= 0; i < bufFactors.size(); i++)
4021    factors.append (bufFactors[i]);
4022  return;
4023}
4024
4025
4026/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
4027/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
4028static inline
4029CFList
4030diophantine (const CFList& recResult, const CFList& factors,
4031             const CFList& products, const CFList& M, const CanonicalForm& E,
4032             bool& bad)
4033{
4034  if (M.isEmpty())
4035  {
4036    CFList result;
4037    CFListIterator j= factors;
4038    CanonicalForm buf;
4039    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
4040    {
4041      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
4042              "constant or univariate poly expected");
4043      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
4044              "constant or univariate poly expected");
4045      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
4046              "constant or univariate poly expected");
4047      buf= mulNTL (E, i.getItem());
4048      result.append (modNTL (buf, j.getItem()));
4049    }
4050    return result;
4051  }
4052  Variable y= M.getLast().mvar();
4053  CFList bufFactors= factors;
4054  for (CFListIterator i= bufFactors; i.hasItem(); i++)
4055    i.getItem()= mod (i.getItem(), y);
4056  CFList bufProducts= products;
4057  for (CFListIterator i= bufProducts; i.hasItem(); i++)
4058    i.getItem()= mod (i.getItem(), y);
4059  CFList buf= M;
4060  buf.removeLast();
4061  CanonicalForm bufE= mod (E, y);
4062  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
4063                                      bufE, bad);
4064
4065  if (bad)
4066    return CFList();
4067
4068  CanonicalForm e= E;
4069  CFListIterator j= products;
4070  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
4071    e -= i.getItem()*j.getItem();
4072
4073  CFList result= recDiophantine;
4074  int d= degree (M.getLast());
4075  CanonicalForm coeffE;
4076  for (int i= 1; i < d; i++)
4077  {
4078    if (degree (e, y) > 0)
4079      coeffE= e[i];
4080    else
4081      coeffE= 0;
4082    if (!coeffE.isZero())
4083    {
4084      CFListIterator k= result;
4085      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
4086                                   coeffE, bad);
4087      if (bad)
4088        return CFList();
4089      CFListIterator l= products;
4090      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
4091      {
4092        k.getItem() += j.getItem()*power (y, i);
4093        e -= j.getItem()*power (y, i)*l.getItem();
4094      }
4095    }
4096    if (e.isZero())
4097      break;
4098  }
4099  if (!e.isZero())
4100  {
4101    bad= true;
4102    return CFList();
4103  }
4104  return result;
4105}
4106
4107static inline
4108void
4109henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
4110            const CFList& diophant, CFMatrix& M, CFArray& Pi,
4111            const CFList& products, int j, const CFList& MOD, bool& noOneToOne)
4112{
4113  CanonicalForm E;
4114  CanonicalForm xToJ= power (F.mvar(), j);
4115  Variable x= F.mvar();
4116
4117  // compute the error
4118#ifdef DEBUGOUTPUT
4119    CanonicalForm test= 1;
4120    for (int i= 0; i < factors.length(); i++)
4121    {
4122      if (i == 0)
4123        test *= mod (bufFactors [i], power (x, j));
4124      else
4125        test *= bufFactors[i];
4126    }
4127    test= mod (test, power (x, j));
4128    test= mod (test, MOD);
4129    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
4130    DEBOUTLN (cerr, "test= " << test2);
4131#endif
4132
4133  if (degree (Pi [factors.length() - 2], x) > 0)
4134    E= F[j] - Pi [factors.length() - 2] [j];
4135  else
4136    E= F[j];
4137
4138  CFArray buf= CFArray (diophant.length());
4139
4140  // actual lifting
4141  CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
4142                                    noOneToOne);
4143
4144  if (noOneToOne)
4145    return;
4146
4147  int k= 0;
4148  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
4149  {
4150    buf[k]= i.getItem();
4151    bufFactors[k] += xToJ*i.getItem();
4152  }
4153
4154  // update Pi [0]
4155  int degBuf0= degree (bufFactors[0], x);
4156  int degBuf1= degree (bufFactors[1], x);
4157  if (degBuf0 > 0 && degBuf1 > 0)
4158  {
4159    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
4160    if (j + 2 <= M.rows())
4161      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
4162  }
4163  CanonicalForm uIZeroJ;
4164
4165  if (degBuf0 > 0 && degBuf1 > 0)
4166    uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
4167             mulMod (bufFactors[1] [0], buf[0], MOD);
4168  else if (degBuf0 > 0)
4169    uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
4170  else if (degBuf1 > 0)
4171    uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
4172  else
4173    uIZeroJ= 0;
4174  Pi [0] += xToJ*uIZeroJ;
4175
4176  CFArray tmp= CFArray (factors.length() - 1);
4177  for (k= 0; k < factors.length() - 1; k++)
4178    tmp[k]= 0;
4179  CFIterator one, two;
4180  one= bufFactors [0];
4181  two= bufFactors [1];
4182  if (degBuf0 > 0 && degBuf1 > 0)
4183  {
4184    while (one.hasTerms() && one.exp() > j) one++;
4185    while (two.hasTerms() && two.exp() > j) two++;
4186    for (k= 1; k <= (int) ceil (j/2.0); k++)
4187    {
4188      if (k != j - k + 1)
4189      {
4190        if ((one.hasTerms() && one.exp() == j - k + 1) &&
4191            (two.hasTerms() && two.exp() == j - k + 1))
4192        {
4193          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
4194                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
4195                    M (j - k + 2, 1);
4196          one++;
4197          two++;
4198        }
4199        else if (one.hasTerms() && one.exp() == j - k + 1)
4200        {
4201          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
4202                    bufFactors[1] [k], MOD) - M (k + 1, 1);
4203          one++;
4204        }
4205        else if (two.hasTerms() && two.exp() == j - k + 1)
4206        {
4207          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
4208                    two.coeff()), MOD) - M (k + 1, 1);
4209          two++;
4210        }
4211      }
4212      else
4213      {
4214        tmp[0] += M (k + 1, 1);
4215      }
4216    }
4217  }
4218
4219  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
4220  {
4221    if (j + 2 <= M.rows())
4222      tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
4223                         (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
4224                         - M(1,1) - M (j + 2,1);
4225  }
4226  else if (degBuf0 >= j + 1)
4227  {
4228    if (degBuf1 > 0)
4229      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
4230    else
4231      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
4232  }
4233  else if (degBuf1 >= j + 1)
4234  {
4235    if (degBuf0 > 0)
4236      tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
4237    else
4238      tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
4239  }
4240  Pi [0] += tmp[0]*xToJ*F.mvar();
4241
4242  // update Pi [l]
4243  int degPi, degBuf;
4244  for (int l= 1; l < factors.length() - 1; l++)
4245  {
4246    degPi= degree (Pi [l - 1], x);
4247    degBuf= degree (bufFactors[l + 1], x);
4248    if (degPi > 0 && degBuf > 0)
4249    {
4250      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
4251      if (j + 2 <= M.rows())
4252        M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
4253                                  MOD);
4254    }
4255
4256    if (degPi > 0 && degBuf > 0)
4257      uIZeroJ= mulMod (Pi[l -1] [0], buf[l + 1], MOD) +
4258               mulMod (uIZeroJ, bufFactors[l+1] [0], MOD);
4259    else if (degPi > 0)
4260      uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD);
4261    else if (degBuf > 0)
4262      uIZeroJ= mulMod (Pi[l - 1], buf[1], MOD);
4263    else
4264      uIZeroJ= 0;
4265
4266    Pi [l] += xToJ*uIZeroJ;
4267
4268    one= bufFactors [l + 1];
4269    two= Pi [l - 1];
4270    if (degBuf > 0 && degPi > 0)
4271    {
4272      while (one.hasTerms() && one.exp() > j) one++;
4273      while (two.hasTerms() && two.exp() > j) two++;
4274      for (k= 1; k <= (int) ceil (j/2.0); k++)
4275      {
4276        if (k != j - k + 1)
4277        {
4278          if ((one.hasTerms() && one.exp() == j - k + 1) &&
4279              (two.hasTerms() && two.exp() == j - k + 1))
4280          {
4281            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
4282                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
4283                      M (j - k + 2, l + 1);
4284            one++;
4285            two++;
4286          }
4287          else if (one.hasTerms() && one.exp() == j - k + 1)
4288          {
4289            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
4290                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
4291            one++;
4292          }
4293          else if (two.hasTerms() && two.exp() == j - k + 1)
4294          {
4295            tmp[l] += mulMod (bufFactors[l + 1] [k],
4296                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
4297            two++;
4298           }
4299        }
4300        else
4301          tmp[l] += M (k + 1, l + 1);
4302      }
4303    }
4304
4305    if (degPi >= j + 1 && degBuf >= j + 1)
4306    {
4307      if (j + 2 <= M.rows())
4308        tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
4309                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
4310                            , MOD) - M(1,l+1) - M (j + 2,l+1);
4311    }
4312    else if (degPi >= j + 1)
4313    {
4314      if (degBuf > 0)
4315        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
4316      else
4317        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
4318    }
4319    else if (degBuf >= j + 1)
4320    {
4321      if (degPi > 0)
4322        tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
4323      else
4324        tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
4325    }
4326
4327    Pi[l] += tmp[l]*xToJ*F.mvar();
4328  }
4329  return;
4330}
4331
4332// wrt. Variable (1)
4333CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
4334{
4335  if (degree (F, 1) <= 0)
4336   return c;
4337  else
4338  {
4339    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
4340    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
4341              - LC (result))*power (result.mvar(), degree (result));
4342    return swapvar (result, Variable (F.level() + 1), Variable (1));
4343  }
4344}
4345
4346CFList
4347henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
4348              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
4349              const CFList& LCs2, bool& bad)
4350{
4351  CFList buf= factors;
4352  int k= 0;
4353  int liftBoundBivar= l[k];
4354  CFList bufbuf= factors;
4355  Variable v= Variable (2);
4356
4357  CFList MOD;
4358  MOD.append (power (Variable (2), liftBoundBivar));
4359  CFArray bufFactors= CFArray (factors.length());
4360  k= 0;
4361  CFListIterator j= eval;
4362  j++;
4363  CFListIterator iter1= LCs1;
4364  CFListIterator iter2= LCs2;
4365  iter1++;
4366  iter2++;
4367  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
4368  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
4369
4370  CFListIterator i= buf;
4371  i++;
4372  Variable y= j.getItem().mvar();
4373  if (y.level() != 3)
4374    y= Variable (3);
4375
4376  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
4377  M (1, 1)= Pi[0];
4378  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
4379    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
4380                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
4381  else if (degree (bufFactors[0], y) > 0)
4382    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
4383  else if (degree (bufFactors[1], y) > 0)
4384    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
4385
4386  CFList products;
4387  for (int i= 0; i < bufFactors.size(); i++)
4388  {
4389    if (degree (bufFactors[i], y) > 0)
4390      products.append (eval.getFirst()/bufFactors[i] [0]);
4391    else
4392      products.append (eval.getFirst()/bufFactors[i]);
4393  }
4394
4395  for (int d= 1; d < l[1]; d++)
4396  {
4397    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
4398    if (bad)
4399      return CFList();
4400  }
4401  CFList result;
4402  for (k= 0; k < factors.length(); k++)
4403    result.append (bufFactors[k]);
4404  return result;
4405}
4406
4407
4408CFList
4409henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
4410            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
4411            lNew, const CFList& LCs1, const CFList& LCs2, bool& bad)
4412{
4413  int k= 0;
4414  CFArray bufFactors= CFArray (factors.length());
4415  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
4416  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
4417  CFList buf= factors;
4418  Variable y= F.getLast().mvar();
4419  Variable x= F.getFirst().mvar();
4420  CanonicalForm xToLOld= power (x, lOld);
4421  Pi [0]= mod (Pi[0], xToLOld);
4422  M (1, 1)= Pi [0];
4423
4424  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
4425    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
4426                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
4427  else if (degree (bufFactors[0], y) > 0)
4428    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
4429  else if (degree (bufFactors[1], y) > 0)
4430    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
4431
4432  CFList products;
4433  CanonicalForm quot;
4434  for (int i= 0; i < bufFactors.size(); i++)
4435  {
4436    if (degree (bufFactors[i], y) > 0)
4437    {
4438      if (!fdivides (bufFactors[i] [0], F.getFirst(), quot))
4439      {
4440        bad= true;
4441        return CFList();
4442      }
4443      products.append (quot);
4444    }
4445    else
4446    {
4447      if (!fdivides (bufFactors[i], F.getFirst(), quot))
4448      {
4449        bad= true;
4450        return CFList();
4451      }
4452      products.append (quot);
4453    }
4454  }
4455
4456  for (int d= 1; d < lNew; d++)
4457  {
4458    henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
4459    if (bad)
4460      return CFList();
4461  }
4462
4463  CFList result;
4464  for (k= 0; k < factors.length(); k++)
4465    result.append (bufFactors[k]);
4466  return result;
4467}
4468
4469CFList
4470henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
4471             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
4472             const CFArray& Pi, const CFList& diophant, bool& bad)
4473{
4474  CFList bufDiophant= diophant;
4475  CFList buf= factors;
4476  if (sort)
4477    sortList (buf, Variable (1));
4478  CFArray bufPi= Pi;
4479  CFMatrix M= CFMatrix (l[1], factors.length());
4480  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2,
4481                               bad);
4482  if (bad)
4483    return CFList();
4484
4485  if (eval.length() == 2)
4486    return result;
4487  CFList MOD;
4488  for (int i= 0; i < 2; i++)
4489    MOD.append (power (Variable (i + 2), l[i]));
4490  CFListIterator j= eval;
4491  j++;
4492  CFList bufEval;
4493  bufEval.append (j.getItem());
4494  j++;
4495  CFListIterator jj= LCs1;
4496  CFListIterator jjj= LCs2;
4497  CFList bufLCs1, bufLCs2;
4498  jj++, jjj++;
4499  bufLCs1.append (jj.getItem());
4500  bufLCs2.append (jjj.getItem());
4501  jj++, jjj++;
4502
4503  for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
4504  {
4505    bufEval.append (j.getItem());
4506    bufLCs1.append (jj.getItem());
4507    bufLCs2.append (jjj.getItem());
4508    M= CFMatrix (l[i], factors.length());
4509    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
4510                         l[i], bufLCs1, bufLCs2, bad);
4511    if (bad)
4512      return CFList();
4513    MOD.append (power (Variable (i + 2), l[i]));
4514    bufEval.removeFirst();
4515    bufLCs1.removeFirst();
4516    bufLCs2.removeFirst();
4517  }
4518  return result;
4519}
4520
4521CFList
4522nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
4523                      CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
4524                      int bivarLiftBound, bool& bad)
4525{
4526  CFList bufFactors2= factors;
4527
4528  Variable y= Variable (2);
4529  for (CFListIterator i= bufFactors2; i.hasItem(); i++)
4530    i.getItem()= mod (i.getItem(), y);
4531
4532  CanonicalForm bufF= F;
4533  bufF= mod (bufF, y);
4534  bufF= mod (bufF, Variable (3));
4535
4536  diophant= diophantine (bufF, bufFactors2);
4537
4538  CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
4539
4540  Pi= CFArray (bufFactors2.length() - 1);
4541
4542  CFArray bufFactors= CFArray (bufFactors2.length());
4543  CFListIterator j= LCs;
4544  int i= 0;
4545  for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
4546    bufFactors[i]= replaceLC (k.getItem(), j.getItem());
4547
4548  //initialise Pi
4549  Variable v= Variable (3);
4550  CanonicalForm yToL= power (y, bivarLiftBound);
4551  if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
4552  {
4553    M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
4554    Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
4555                       mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
4556  }
4557  else if (degree (bufFactors[0], v) > 0)
4558  {
4559    M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
4560    Pi [0]=  M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
4561  }
4562  else if (degree (bufFactors[1], v) > 0)
4563  {
4564    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
4565    Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
4566  }
4567  else
4568  {
4569    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
4570    Pi [0]= M (1,1);
4571  }
4572
4573  for (i= 1; i < Pi.size(); i++)
4574  {
4575    if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
4576    {
4577      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
4578      Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
4579                       mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
4580    }
4581    else if (degree (Pi[i-1], v) > 0)
4582    {
4583      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
4584      Pi [i]=  M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
4585    }
4586    else if (degree (bufFactors[i+1], v) > 0)
4587    {
4588      M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
4589      Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
4590    }
4591    else
4592    {
4593      M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
4594      Pi [i]= M (1,i+1);
4595    }
4596  }
4597
4598  CFList products;
4599  bufF= mod (F, Variable (3));
4600  for (CFListIterator k= factors; k.hasItem(); k++)
4601    products.append (bufF/k.getItem());
4602
4603  CFList MOD= CFList (power (v, liftBound));
4604  MOD.insert (yToL);
4605  for (int d= 1; d < liftBound; d++)
4606  {
4607    henselStep2 (F, factors, bufFactors, diophant, M, Pi, products, d, MOD, bad);
4608    if (bad)
4609      return CFList();
4610  }
4611
4612  CFList result;
4613  for (i= 0; i < factors.length(); i++)
4614    result.append (bufFactors[i]);
4615  return result;
4616}
4617
4618CFList
4619nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
4620                    CFList& diophant, CFArray& Pi, CFMatrix& M, const int lOld,
4621                    int& lNew, const CFList& MOD, bool& noOneToOne
4622                   )
4623{
4624
4625  int k= 0;
4626  CFArray bufFactors= CFArray (factors.length());
4627  CFListIterator j= LCs;
4628  for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
4629    bufFactors [k]= replaceLC (i.getItem(), j.getItem());
4630
4631  Variable y= F.getLast().mvar();
4632  Variable x= F.getFirst().mvar();
4633  CanonicalForm xToLOld= power (x, lOld);
4634
4635  Pi [0]= mod (Pi[0], xToLOld);
4636  M (1, 1)= Pi [0];
4637
4638  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
4639    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
4640                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
4641  else if (degree (bufFactors[0], y) > 0)
4642    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
4643  else if (degree (bufFactors[1], y) > 0)
4644    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
4645
4646  for (int i= 1; i < Pi.size(); i++)
4647  {
4648    Pi [i]= mod (Pi [i], xToLOld);
4649    M (1, i + 1)= Pi [i];
4650
4651    if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
4652      Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
4653                 mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
4654    else if (degree (Pi[i-1], y) > 0)
4655      Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
4656    else if (degree (bufFactors[i+1], y) > 0)
4657      Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
4658  }
4659
4660  CFList products;
4661  CanonicalForm quot, bufF= F.getFirst();
4662
4663  for (int i= 0; i < bufFactors.size(); i++)
4664  {
4665    if (degree (bufFactors[i], y) > 0)
4666    {
4667      if (!fdivides (bufFactors[i] [0], bufF, quot))
4668      {
4669        noOneToOne= true;
4670        return factors;
4671      }
4672      products.append (quot);
4673    }
4674    else
4675    {
4676      if (!fdivides (bufFactors[i], bufF, quot))
4677      {
4678        noOneToOne= true;
4679        return factors;
4680      }
4681      products.append (quot);
4682    }
4683  }
4684
4685  for (int d= 1; d < lNew; d++)
4686  {
4687    henselStep2 (F.getLast(), factors, bufFactors, diophant, M, Pi, products, d,
4688                 MOD, noOneToOne);
4689    if (noOneToOne)
4690      return CFList();
4691  }
4692
4693  CFList result;
4694  for (k= 0; k < factors.length(); k++)
4695    result.append (bufFactors[k]);
4696  return result;
4697}
4698
4699CFList
4700nonMonicHenselLift (const CFList& eval, const CFList& factors,
4701                    CFList* const& LCs, CFList& diophant, CFArray& Pi,
4702                    int* liftBound, int length, bool& noOneToOne
4703                   )
4704{
4705  CFList bufDiophant= diophant;
4706  CFList buf= factors;
4707  CFArray bufPi= Pi;
4708  CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
4709  int k= 0;
4710
4711  CFList result=
4712  nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
4713                        liftBound[1], liftBound[0], noOneToOne);
4714
4715  if (noOneToOne)
4716    return CFList();
4717
4718  if (eval.length() == 1)
4719    return result;
4720
4721  k++;
4722  CFList MOD;
4723  for (int i= 0; i < 2; i++)
4724    MOD.append (power (Variable (i + 2), liftBound[i]));
4725
4726  CFListIterator j= eval;
4727  CFList bufEval;
4728  bufEval.append (j.getItem());
4729  j++;
4730
4731  for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
4732  {
4733    bufEval.append (j.getItem());
4734    M= CFMatrix (liftBound[i], factors.length() - 1);
4735    result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
4736                                liftBound[i-1], liftBound[i], MOD, noOneToOne);
4737    if (noOneToOne)
4738      return result;
4739    MOD.append (power (Variable (i + 2), liftBound[i]));
4740    bufEval.removeFirst();
4741  }
4742
4743  return result;
4744}
4745
4746#endif
4747/* HAVE_NTL */
4748
Note: See TracBrowser for help on using the repository browser.