source: git/factory/facHensel.cc @ 2537fa0

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