source: git/factory/facHensel.cc @ 86faff

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