source: git/factory/facHensel.cc @ 27ab36

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