source: git/factory/facHensel.cc @ a37b34

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