source: git/factory/facHensel.cc @ 90c9e3

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