source: git/factory/facHensel.cc @ 8e51ca

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