source: git/factory/facHensel.cc @ 9ebec2

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