source: git/factory/facHensel.cc @ 428b38

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