source: git/factory/facHensel.cc @ c9733f

spielwiese
Last change on this file since c9733f was c9733f, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: add new function diophantineQa
  • Property mode set to 100644
File size: 74.6 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facHensel.cc
5 *
6 * This file implements functions to lift factors via Hensel lifting.
7 *
8 * ABSTRACT: Hensel lifting is described in "Efficient Multivariate
9 * Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with
10 * remainder is described in "Fast Recursive Division" by C. Burnikel and
11 * J. Ziegler. Karatsuba multiplication is described in "Modern Computer
12 * Algebra" by J. von zur Gathen and J. Gerhard.
13 *
14 * @author Martin Lee
15 *
16 **/
17/*****************************************************************************/
18
19#include "config.h"
20
21#include "cf_assert.h"
22#include "debug.h"
23#include "timing.h"
24
25#include "algext.h"
26#include "facHensel.h"
27#include "facMul.h"
28#include "fac_util.h"
29#include "cf_algorithm.h"
30#include "cf_primes.h"
31#include "facBivar.h"
32#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);
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  if (j == 1)
921  {
922    if (degBuf0 > 0 && degBuf1 > 0)
923      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
924                  (bufFactors[1] [0] + buf[1]), b) - M(1, 1) - M(j + 1, 1);
925    else if (degBuf0 > 0)
926      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1], b);
927    else if (degBuf1 > 0)
928      uIZeroJ= mulNTL (bufFactors[0], buf[1], b);
929    else
930      uIZeroJ= 0;
931    if (b.getp() != 0)
932      uIZeroJ= b (uIZeroJ);
933    Pi [0] += xToJ*uIZeroJ;
934    if (b.getp() != 0)
935      Pi [0]= b (Pi[0]);
936  }
937  else
938  {
939    if (degBuf0 > 0 && degBuf1 > 0)
940      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
941                  (bufFactors[1] [0] + buf[1]), b) - M(1, 1) - M(j + 1, 1);
942    else if (degBuf0 > 0)
943      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1], b);
944    else if (degBuf1 > 0)
945      uIZeroJ= mulNTL (bufFactors[0], buf[1], b);
946    else
947      uIZeroJ= 0;
948    if (b.getp() != 0)
949      uIZeroJ= b (uIZeroJ);
950    Pi [0] += xToJ*uIZeroJ;
951    if (b.getp() != 0)
952      Pi [0]= b (Pi[0]);
953  }
954  CFArray tmp= CFArray (factors.length() - 1);
955  for (k= 0; k < factors.length() - 1; k++)
956    tmp[k]= 0;
957  CFIterator one, two;
958  one= bufFactors [0];
959  two= bufFactors [1];
960  if (degBuf0 > 0 && degBuf1 > 0)
961  {
962    for (k= 1; k <= (int) ceil (j/2.0); k++)
963    {
964      if (k != j - k + 1)
965      {
966        if ((one.hasTerms() && one.exp() == j - k + 1)
967             && (two.hasTerms() && two.exp() == j - k + 1))
968        {
969          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), (bufFactors[1][k]+
970                     two.coeff()), b) - M (k + 1, 1) - M (j - k + 2, 1);
971          one++;
972          two++;
973        }
974        else if (one.hasTerms() && one.exp() == j - k + 1)
975        {
976          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1][k], b)
977                    - M (k + 1, 1);
978          one++;
979        }
980        else if (two.hasTerms() && two.exp() == j - k + 1)
981        {
982          tmp[0] += mulNTL (bufFactors[0][k], (bufFactors[1][k]+two.coeff()), b)
983                    - M (k + 1, 1);
984          two++;
985        }
986      }
987      else
988      {
989        tmp[0] += M (k + 1, 1);
990      }
991    }
992  }
993  if (b.getp() != 0)
994    tmp[0]= b (tmp[0]);
995  Pi [0] += tmp[0]*xToJ*F.mvar();
996
997  // update Pi [l]
998  int degPi, degBuf;
999  for (int l= 1; l < factors.length() - 1; l++)
1000  {
1001    degPi= degree (Pi [l - 1], x);
1002    degBuf= degree (bufFactors[l + 1], x);
1003    if (degPi > 0 && degBuf > 0)
1004      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j], b);
1005    if (j == 1)
1006    {
1007      if (degPi > 0 && degBuf > 0)
1008        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
1009                  bufFactors[l + 1] [0] + buf[l + 1], b) - M (j + 1, l +1) -
1010                  M (1, l + 1));
1011      else if (degPi > 0)
1012        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1], b));
1013      else if (degBuf > 0)
1014        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1], b));
1015    }
1016    else
1017    {
1018      if (degPi > 0 && degBuf > 0)
1019      {
1020        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0], b);
1021        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1], b);
1022      }
1023      else if (degPi > 0)
1024        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1], b);
1025      else if (degBuf > 0)
1026      {
1027        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0], b);
1028        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1], b);
1029      }
1030      Pi[l] += xToJ*uIZeroJ;
1031    }
1032    one= bufFactors [l + 1];
1033    two= Pi [l - 1];
1034    if (two.hasTerms() && two.exp() == j + 1)
1035    {
1036      if (degBuf > 0 && degPi > 0)
1037      {
1038          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0], b);
1039          two++;
1040      }
1041      else if (degPi > 0)
1042      {
1043          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1], b);
1044          two++;
1045      }
1046    }
1047    if (degBuf > 0 && degPi > 0)
1048    {
1049      for (k= 1; k <= (int) ceil (j/2.0); k++)
1050      {
1051        if (k != j - k + 1)
1052        {
1053          if ((one.hasTerms() && one.exp() == j - k + 1) &&
1054              (two.hasTerms() && two.exp() == j - k + 1))
1055          {
1056            tmp[l] += mulNTL ((bufFactors[l+1][k] + one.coeff()), (Pi[l-1][k] +
1057                      two.coeff()),b) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
1058            one++;
1059            two++;
1060          }
1061          else if (one.hasTerms() && one.exp() == j - k + 1)
1062          {
1063            tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()), Pi[l-1][k], b) -
1064                       M (k + 1, l + 1);
1065            one++;
1066          }
1067          else if (two.hasTerms() && two.exp() == j - k + 1)
1068          {
1069            tmp[l] += mulNTL (bufFactors[l+1][k], (Pi[l-1][k] + two.coeff()), b)
1070                      - M (k + 1, l + 1);
1071            two++;
1072          }
1073        }
1074        else
1075          tmp[l] += M (k + 1, l + 1);
1076      }
1077    }
1078    if (b.getp() != 0)
1079      tmp[l]= b (tmp[l]);
1080    Pi[l] += tmp[l]*xToJ*F.mvar();
1081  }
1082}
1083
1084void
1085henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1086              CFList& diophant, CFMatrix& M, modpk& b, bool sort)
1087{
1088  if (sort)
1089    sortList (factors, Variable (1));
1090  Pi= CFArray (factors.length() - 1);
1091  CFListIterator j= factors;
1092  diophant= diophantine (F[0], F, factors, b);
1093  CanonicalForm bufF= F;
1094  if (getCharacteristic() == 0 && b.getp() != 0)
1095  {
1096    Variable v;
1097    bool hasAlgVar= hasFirstAlgVar (F, v);
1098    for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++)
1099      hasAlgVar= hasFirstAlgVar (i.getItem(), v);
1100    Variable w;
1101    bool hasAlgVar2= false;
1102    for (CFListIterator i= diophant; i.hasItem() && !hasAlgVar2; i++)
1103      hasAlgVar2= hasFirstAlgVar (i.getItem(), w);
1104    if (hasAlgVar && hasAlgVar2 && v!=w)
1105    {
1106      bufF= replacevar (bufF, v, w);
1107      for (CFListIterator i= factors; i.hasItem(); i++)
1108        i.getItem()= replacevar (i.getItem(), v, w);
1109    }
1110  }
1111
1112  DEBOUTLN (cerr, "diophant= " << diophant);
1113  j++;
1114  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()), b);
1115  M (1, 1)= Pi [0];
1116  int i= 1;
1117  if (j.hasItem())
1118    j++;
1119  for (; j.hasItem(); j++, i++)
1120  {
1121    Pi [i]= mulNTL (Pi [i - 1], j.getItem(), b);
1122    M (1, i + 1)= Pi [i];
1123  }
1124  CFArray bufFactors= CFArray (factors.length());
1125  i= 0;
1126  for (CFListIterator k= factors; k.hasItem(); i++, k++)
1127  {
1128    if (i == 0)
1129      bufFactors[i]= mod (k.getItem(), F.mvar());
1130    else
1131      bufFactors[i]= k.getItem();
1132  }
1133  for (i= 1; i < l; i++)
1134    henselStep12 (bufF, factors, bufFactors, diophant, M, Pi, i, b);
1135
1136  CFListIterator k= factors;
1137  for (i= 0; i < factors.length (); i++, k++)
1138    k.getItem()= bufFactors[i];
1139  factors.removeFirst();
1140}
1141
1142void
1143henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1144              CFList& diophant, CFMatrix& M, bool sort)
1145{
1146  modpk dummy= modpk();
1147  henselLift12 (F, factors, l, Pi, diophant, M, dummy, sort);
1148}
1149
1150void
1151henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
1152                    end, CFArray& Pi, const CFList& diophant, CFMatrix& M,
1153                    const modpk& b)
1154{
1155  CFArray bufFactors= CFArray (factors.length());
1156  int i= 0;
1157  CanonicalForm xToStart= power (F.mvar(), start);
1158  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1159  {
1160    if (i == 0)
1161      bufFactors[i]= mod (k.getItem(), xToStart);
1162    else
1163      bufFactors[i]= k.getItem();
1164  }
1165  for (i= start; i < end; i++)
1166    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i, b);
1167
1168  CFListIterator k= factors;
1169  for (i= 0; i < factors.length(); k++, i++)
1170    k.getItem()= bufFactors [i];
1171  factors.removeFirst();
1172  return;
1173}
1174
1175CFList
1176biDiophantine (const CanonicalForm& F, const CFList& factors, int d)
1177{
1178  Variable y= F.mvar();
1179  CFList result;
1180  if (y.level() == 1)
1181  {
1182    result= diophantine (F, factors);
1183    return result;
1184  }
1185  else
1186  {
1187    CFList buf= factors;
1188    for (CFListIterator i= buf; i.hasItem(); i++)
1189      i.getItem()= mod (i.getItem(), y);
1190    CanonicalForm A= mod (F, y);
1191    int bufD= 1;
1192    CFList recResult= biDiophantine (A, buf, bufD);
1193    CanonicalForm e= 1;
1194    CFList p;
1195    CFArray bufFactors= CFArray (factors.length());
1196    CanonicalForm yToD= power (y, d);
1197    int k= 0;
1198    for (CFListIterator i= factors; i.hasItem(); i++, k++)
1199    {
1200      bufFactors [k]= i.getItem();
1201    }
1202    CanonicalForm b, quot;
1203    for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1204    {
1205      b= 1;
1206      if (fdivides (bufFactors[k], F, quot))
1207        b= quot;
1208      else
1209      {
1210        for (int l= 0; l < factors.length(); l++)
1211        {
1212          if (l == k)
1213            continue;
1214          else
1215          {
1216            b= mulMod2 (b, bufFactors[l], yToD);
1217          }
1218        }
1219      }
1220      p.append (b);
1221    }
1222
1223    CFListIterator j= p;
1224    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1225      e -= i.getItem()*j.getItem();
1226
1227    if (e.isZero())
1228      return recResult;
1229    CanonicalForm coeffE;
1230    CFList s;
1231    result= recResult;
1232    CanonicalForm g;
1233    for (int i= 1; i < d; i++)
1234    {
1235      if (degree (e, y) > 0)
1236        coeffE= e[i];
1237      else
1238        coeffE= 0;
1239      if (!coeffE.isZero())
1240      {
1241        CFListIterator k= result;
1242        CFListIterator l= p;
1243        int ii= 0;
1244        j= recResult;
1245        for (; j.hasItem(); j++, k++, l++, ii++)
1246        {
1247          g= coeffE*j.getItem();
1248          if (degree (bufFactors[ii], y) <= 0)
1249            g= mod (g, bufFactors[ii]);
1250          else
1251            g= mod (g, bufFactors[ii][0]);
1252          k.getItem() += g*power (y, i);
1253          e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
1254          DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
1255                    mod (e, power (y, i + 1)));
1256        }
1257      }
1258      if (e.isZero())
1259        break;
1260    }
1261
1262    DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
1263
1264#ifdef DEBUGOUTPUT
1265    CanonicalForm test= 0;
1266    j= p;
1267    for (CFListIterator i= result; i.hasItem(); i++, j++)
1268      test += mod (i.getItem()*j.getItem(), power (y, d));
1269    DEBOUTLN (cerr, "test= " << test);
1270#endif
1271    return result;
1272  }
1273}
1274
1275CFList
1276multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
1277                     const CFList& recResult, const CFList& M, int d)
1278{
1279  Variable y= F.mvar();
1280  CFList result;
1281  CFListIterator i;
1282  CanonicalForm e= 1;
1283  CFListIterator j= factors;
1284  CFList p;
1285  CFArray bufFactors= CFArray (factors.length());
1286  CanonicalForm yToD= power (y, d);
1287  int k= 0;
1288  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1289    bufFactors [k]= i.getItem();
1290  CanonicalForm b, quot;
1291  CFList buf= M;
1292  buf.removeLast();
1293  buf.append (yToD);
1294  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1295  {
1296    b= 1;
1297    if (fdivides (bufFactors[k], F, quot))
1298      b= quot;
1299    else
1300    {
1301      for (int l= 0; l < factors.length(); l++)
1302      {
1303        if (l == k)
1304          continue;
1305        else
1306        {
1307          b= mulMod (b, bufFactors[l], buf);
1308        }
1309      }
1310    }
1311    p.append (b);
1312  }
1313  j= p;
1314  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1315    e -= mulMod (i.getItem(), j.getItem(), M);
1316
1317  if (e.isZero())
1318    return recResult;
1319  CanonicalForm coeffE;
1320  CFList s;
1321  result= recResult;
1322  CanonicalForm g;
1323  for (int i= 1; i < d; i++)
1324  {
1325    if (degree (e, y) > 0)
1326      coeffE= e[i];
1327    else
1328      coeffE= 0;
1329    if (!coeffE.isZero())
1330    {
1331      CFListIterator k= result;
1332      CFListIterator l= p;
1333      j= recResult;
1334      int ii= 0;
1335      CanonicalForm dummy;
1336      for (; j.hasItem(); j++, k++, l++, ii++)
1337      {
1338        g= mulMod (coeffE, j.getItem(), M);
1339        if (degree (bufFactors[ii], y) <= 0)
1340          divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
1341                  g, M);
1342        else
1343          divrem (g, bufFactors[ii][0], dummy, g, M);
1344        k.getItem() += g*power (y, i);
1345        e -= mulMod (g*power (y, i), l.getItem(), M);
1346      }
1347    }
1348
1349    if (e.isZero())
1350      break;
1351  }
1352
1353#ifdef DEBUGOUTPUT
1354    CanonicalForm test= 0;
1355    j= p;
1356    for (CFListIterator i= result; i.hasItem(); i++, j++)
1357      test += mod (i.getItem()*j.getItem(), power (y, d));
1358    DEBOUTLN (cerr, "test in multiRecDiophantine= " << test);
1359#endif
1360  return result;
1361}
1362
1363static
1364void
1365henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
1366            const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
1367            const CFList& MOD)
1368{
1369  CanonicalForm E;
1370  CanonicalForm xToJ= power (F.mvar(), j);
1371  Variable x= F.mvar();
1372  // compute the error
1373  if (j == 1)
1374  {
1375    E= F[j];
1376#ifdef DEBUGOUTPUT
1377    CanonicalForm test= 1;
1378    for (int i= 0; i < factors.length(); i++)
1379    {
1380      if (i == 0)
1381        test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
1382      else
1383        test= mulMod (test, bufFactors[i], MOD);
1384    }
1385    CanonicalForm test2= mod (F-test, xToJ);
1386
1387    test2= mod (test2, MOD);
1388    DEBOUTLN (cerr, "test in henselStep= " << test2);
1389#endif
1390  }
1391  else
1392  {
1393#ifdef DEBUGOUTPUT
1394    CanonicalForm test= 1;
1395    for (int i= 0; i < factors.length(); i++)
1396    {
1397      if (i == 0)
1398        test *= mod (bufFactors [i], power (x, j));
1399      else
1400        test *= bufFactors[i];
1401    }
1402    test= mod (test, power (x, j));
1403    test= mod (test, MOD);
1404    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
1405    DEBOUTLN (cerr, "test in henselStep= " << test2);
1406#endif
1407
1408    if (degree (Pi [factors.length() - 2], x) > 0)
1409      E= F[j] - Pi [factors.length() - 2] [j];
1410    else
1411      E= F[j];
1412  }
1413
1414  CFArray buf= CFArray (diophant.length());
1415  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1416  int k= 0;
1417  // actual lifting
1418  CanonicalForm dummy, rest1;
1419  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1420  {
1421    if (degree (bufFactors[k], x) > 0)
1422    {
1423      if (k > 0)
1424        divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
1425      else
1426        rest1= E;
1427    }
1428    else
1429      divrem (E, bufFactors[k], dummy, rest1, MOD);
1430
1431    buf[k]= mulMod (i.getItem(), rest1, MOD);
1432
1433    if (degree (bufFactors[k], x) > 0)
1434      divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
1435    else
1436      divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
1437  }
1438  for (k= 1; k < factors.length(); k++)
1439    bufFactors[k] += xToJ*buf[k];
1440
1441  // update Pi [0]
1442  int degBuf0= degree (bufFactors[0], x);
1443  int degBuf1= degree (bufFactors[1], x);
1444  if (degBuf0 > 0 && degBuf1 > 0)
1445    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
1446  CanonicalForm uIZeroJ;
1447  if (j == 1)
1448  {
1449    if (degBuf0 > 0 && degBuf1 > 0)
1450      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
1451                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
1452    else if (degBuf0 > 0)
1453      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
1454    else if (degBuf1 > 0)
1455      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
1456    else
1457      uIZeroJ= 0;
1458    Pi [0] += xToJ*uIZeroJ;
1459  }
1460  else
1461  {
1462    if (degBuf0 > 0 && degBuf1 > 0)
1463      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
1464                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
1465    else if (degBuf0 > 0)
1466      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
1467    else if (degBuf1 > 0)
1468      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
1469    else
1470      uIZeroJ= 0;
1471    Pi [0] += xToJ*uIZeroJ;
1472  }
1473
1474  CFArray tmp= CFArray (factors.length() - 1);
1475  for (k= 0; k < factors.length() - 1; k++)
1476    tmp[k]= 0;
1477  CFIterator one, two;
1478  one= bufFactors [0];
1479  two= bufFactors [1];
1480  if (degBuf0 > 0 && degBuf1 > 0)
1481  {
1482    for (k= 1; k <= (int) ceil (j/2.0); k++)
1483    {
1484      if (k != j - k + 1)
1485      {
1486        if ((one.hasTerms() && one.exp() == j - k + 1) &&
1487            (two.hasTerms() && two.exp() == j - k + 1))
1488        {
1489          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
1490                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
1491                    M (j - k + 2, 1);
1492          one++;
1493          two++;
1494        }
1495        else if (one.hasTerms() && one.exp() == j - k + 1)
1496        {
1497          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
1498                    bufFactors[1] [k], MOD) - M (k + 1, 1);
1499          one++;
1500        }
1501        else if (two.hasTerms() && two.exp() == j - k + 1)
1502        {
1503          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
1504                    two.coeff()), MOD) - M (k + 1, 1);
1505          two++;
1506        }
1507      }
1508      else
1509      {
1510        tmp[0] += M (k + 1, 1);
1511      }
1512    }
1513  }
1514  Pi [0] += tmp[0]*xToJ*F.mvar();
1515
1516  // update Pi [l]
1517  int degPi, degBuf;
1518  for (int l= 1; l < factors.length() - 1; l++)
1519  {
1520    degPi= degree (Pi [l - 1], x);
1521    degBuf= degree (bufFactors[l + 1], x);
1522    if (degPi > 0 && degBuf > 0)
1523      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
1524    if (j == 1)
1525    {
1526      if (degPi > 0 && degBuf > 0)
1527        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
1528                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
1529                  M (1, l + 1));
1530      else if (degPi > 0)
1531        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
1532      else if (degBuf > 0)
1533        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
1534    }
1535    else
1536    {
1537      if (degPi > 0 && degBuf > 0)
1538      {
1539        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
1540        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
1541      }
1542      else if (degPi > 0)
1543        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
1544      else if (degBuf > 0)
1545      {
1546        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
1547        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
1548      }
1549      Pi[l] += xToJ*uIZeroJ;
1550    }
1551    one= bufFactors [l + 1];
1552    two= Pi [l - 1];
1553    if (two.hasTerms() && two.exp() == j + 1)
1554    {
1555      if (degBuf > 0 && degPi > 0)
1556      {
1557          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
1558          two++;
1559      }
1560      else if (degPi > 0)
1561      {
1562          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
1563          two++;
1564      }
1565    }
1566    if (degBuf > 0 && degPi > 0)
1567    {
1568      for (k= 1; k <= (int) ceil (j/2.0); k++)
1569      {
1570        if (k != j - k + 1)
1571        {
1572          if ((one.hasTerms() && one.exp() == j - k + 1) &&
1573              (two.hasTerms() && two.exp() == j - k + 1))
1574          {
1575            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1576                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
1577                      M (j - k + 2, l + 1);
1578            one++;
1579            two++;
1580          }
1581          else if (one.hasTerms() && one.exp() == j - k + 1)
1582          {
1583            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1584                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
1585            one++;
1586          }
1587          else if (two.hasTerms() && two.exp() == j - k + 1)
1588          {
1589            tmp[l] += mulMod (bufFactors[l + 1] [k],
1590                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
1591            two++;
1592          }
1593        }
1594        else
1595          tmp[l] += M (k + 1, l + 1);
1596      }
1597    }
1598    Pi[l] += tmp[l]*xToJ*F.mvar();
1599  }
1600
1601  return;
1602}
1603
1604CFList
1605henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList&
1606              diophant, CFArray& Pi, CFMatrix& M)
1607{
1608  CFList buf= factors;
1609  int k= 0;
1610  int liftBoundBivar= l[k];
1611  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
1612  CFList MOD;
1613  MOD.append (power (Variable (2), liftBoundBivar));
1614  CFArray bufFactors= CFArray (factors.length());
1615  k= 0;
1616  CFListIterator j= eval;
1617  j++;
1618  buf.removeFirst();
1619  buf.insert (LC (j.getItem(), 1));
1620  for (CFListIterator i= buf; i.hasItem(); i++, k++)
1621    bufFactors[k]= i.getItem();
1622  Pi= CFArray (factors.length() - 1);
1623  CFListIterator i= buf;
1624  i++;
1625  Variable y= j.getItem().mvar();
1626  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
1627  M (1, 1)= Pi [0];
1628  k= 1;
1629  if (i.hasItem())
1630    i++;
1631  for (; i.hasItem(); i++, k++)
1632  {
1633    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
1634    M (1, k + 1)= Pi [k];
1635  }
1636
1637  for (int d= 1; d < l[1]; d++)
1638    henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
1639  CFList result;
1640  for (k= 1; k < factors.length(); k++)
1641    result.append (bufFactors[k]);
1642  return result;
1643}
1644
1645void
1646henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
1647                  CFArray& Pi, const CFList& diophant, CFMatrix& M,
1648                  const CFList& MOD)
1649{
1650  CFArray bufFactors= CFArray (factors.length());
1651  int i= 0;
1652  CanonicalForm xToStart= power (F.mvar(), start);
1653  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1654  {
1655    if (i == 0)
1656      bufFactors[i]= mod (k.getItem(), xToStart);
1657    else
1658      bufFactors[i]= k.getItem();
1659  }
1660  for (i= start; i < end; i++)
1661    henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD);
1662
1663  CFListIterator k= factors;
1664  for (i= 0; i < factors.length(); k++, i++)
1665    k.getItem()= bufFactors [i];
1666  factors.removeFirst();
1667  return;
1668}
1669
1670CFList
1671henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
1672            diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew)
1673{
1674  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
1675  int k= 0;
1676  CFArray bufFactors= CFArray (factors.length());
1677  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1678  {
1679    if (k == 0)
1680      bufFactors[k]= LC (F.getLast(), 1);
1681    else
1682      bufFactors[k]= i.getItem();
1683  }
1684  CFList buf= factors;
1685  buf.removeFirst();
1686  buf.insert (LC (F.getLast(), 1));
1687  CFListIterator i= buf;
1688  i++;
1689  Variable y= F.getLast().mvar();
1690  Variable x= F.getFirst().mvar();
1691  CanonicalForm xToLOld= power (x, lOld);
1692  Pi [0]= mod (Pi[0], xToLOld);
1693  M (1, 1)= Pi [0];
1694  k= 1;
1695  if (i.hasItem())
1696    i++;
1697  for (; i.hasItem(); i++, k++)
1698  {
1699    Pi [k]= mod (Pi [k], xToLOld);
1700    M (1, k + 1)= Pi [k];
1701  }
1702
1703  for (int d= 1; d < lNew; d++)
1704    henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
1705  CFList result;
1706  for (k= 1; k < factors.length(); k++)
1707    result.append (bufFactors[k]);
1708  return result;
1709}
1710
1711CFList
1712henselLift (const CFList& eval, const CFList& factors, int* l, int lLength,
1713            bool sort)
1714{
1715  CFList diophant;
1716  CFList buf= factors;
1717  buf.insert (LC (eval.getFirst(), 1));
1718  if (sort)
1719    sortList (buf, Variable (1));
1720  CFArray Pi;
1721  CFMatrix M= CFMatrix (l[1], factors.length());
1722  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
1723  if (eval.length() == 2)
1724    return result;
1725  CFList MOD;
1726  for (int i= 0; i < 2; i++)
1727    MOD.append (power (Variable (i + 2), l[i]));
1728  CFListIterator j= eval;
1729  j++;
1730  CFList bufEval;
1731  bufEval.append (j.getItem());
1732  j++;
1733
1734  for (int i= 2; i < lLength && j.hasItem(); i++, j++)
1735  {
1736    result.insert (LC (bufEval.getFirst(), 1));
1737    bufEval.append (j.getItem());
1738    M= CFMatrix (l[i], factors.length());
1739    result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
1740    MOD.append (power (Variable (i + 2), l[i]));
1741    bufEval.removeFirst();
1742  }
1743  return result;
1744}
1745
1746// nonmonic
1747
1748void
1749nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors,
1750                      CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1751                      CFArray& Pi, int j, const CFArray& /*LCs*/)
1752{
1753  Variable x= F.mvar();
1754  CanonicalForm xToJ= power (x, j);
1755
1756  CanonicalForm E;
1757  // compute the error
1758  if (degree (Pi [factors.length() - 2], x) > 0)
1759    E= F[j] - Pi [factors.length() - 2] [j];
1760  else
1761    E= F[j];
1762
1763  CFArray buf= CFArray (diophant.length());
1764
1765  int k= 0;
1766  CanonicalForm remainder;
1767  // actual lifting
1768  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1769  {
1770    if (degree (bufFactors[k], x) > 0)
1771      remainder= modNTL (E, bufFactors[k] [0]);
1772    else
1773      remainder= modNTL (E, bufFactors[k]);
1774    buf[k]= mulNTL (i.getItem(), remainder);
1775    if (degree (bufFactors[k], x) > 0)
1776      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
1777    else
1778      buf[k]= modNTL (buf[k], bufFactors[k]);
1779  }
1780
1781  for (k= 0; k < factors.length(); k++)
1782    bufFactors[k] += xToJ*buf[k];
1783
1784  // update Pi [0]
1785  int degBuf0= degree (bufFactors[0], x);
1786  int degBuf1= degree (bufFactors[1], x);
1787  if (degBuf0 > 0 && degBuf1 > 0)
1788  {
1789    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
1790    if (j + 2 <= M.rows())
1791      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
1792  }
1793  CanonicalForm uIZeroJ;
1794  if (degBuf0 > 0 && degBuf1 > 0)
1795    uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]);
1796  else if (degBuf0 > 0)
1797    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
1798  else if (degBuf1 > 0)
1799    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
1800  else
1801    uIZeroJ= 0;
1802  Pi [0] += xToJ*uIZeroJ;
1803
1804  CFArray tmp= CFArray (factors.length() - 1);
1805  for (k= 0; k < factors.length() - 1; k++)
1806    tmp[k]= 0;
1807  CFIterator one, two;
1808  one= bufFactors [0];
1809  two= bufFactors [1];
1810  if (degBuf0 > 0 && degBuf1 > 0)
1811  {
1812    while (one.hasTerms() && one.exp() > j) one++;
1813    while (two.hasTerms() && two.exp() > j) two++;
1814    for (k= 1; k <= (int) ceil (j/2.0); k++)
1815    {
1816      if (one.hasTerms() && two.hasTerms())
1817      {
1818        if (k != j - k + 1)
1819        {
1820          if ((one.hasTerms() && one.exp() == j - k + 1) && +
1821              (two.hasTerms() && two.exp() == j - k + 1))
1822        {
1823          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
1824                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
1825          one++;
1826          two++;
1827        }
1828        else if (one.hasTerms() && one.exp() == j - k + 1)
1829        {
1830          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
1831                      M (k + 1, 1);
1832          one++;
1833        }
1834        else if (two.hasTerms() && two.exp() == j - k + 1)
1835        {
1836          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
1837                    M (k + 1, 1);
1838          two++;
1839        }
1840      }
1841      else
1842        tmp[0] += M (k + 1, 1);
1843      }
1844    }
1845  }
1846
1847  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
1848  {
1849    if (j + 2 <= M.rows())
1850      tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
1851                         (bufFactors [1] [j + 1] + bufFactors [1] [0]))
1852                         - M(1,1) - M (j + 2,1);
1853  }
1854  else if (degBuf0 >= j + 1)
1855  {
1856    if (degBuf1 > 0)
1857      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
1858    else
1859      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]);
1860  }
1861  else if (degBuf1 >= j + 1)
1862  {
1863    if (degBuf0 > 0)
1864      tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
1865    else
1866      tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]);
1867  }
1868
1869  Pi [0] += tmp[0]*xToJ*F.mvar();
1870
1871  int degPi, degBuf;
1872  for (int l= 1; l < factors.length() - 1; l++)
1873  {
1874    degPi= degree (Pi [l - 1], x);
1875    degBuf= degree (bufFactors[l + 1], x);
1876    if (degPi > 0 && degBuf > 0)
1877    {
1878      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
1879      if (j + 2 <= M.rows())
1880        M (j + 2, l + 1)= mulNTL (Pi [l - 1][j + 1], bufFactors[l + 1] [j + 1]);
1881    }
1882
1883    if (degPi > 0 && degBuf > 0)
1884      uIZeroJ= mulNTL (Pi[l -1] [0], buf[l + 1]) +
1885               mulNTL (uIZeroJ, bufFactors[l+1] [0]);
1886    else if (degPi > 0)
1887      uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]);
1888    else if (degBuf > 0)
1889      uIZeroJ= mulNTL (Pi[l - 1], buf[1]);
1890    else
1891      uIZeroJ= 0;
1892
1893    Pi [l] += xToJ*uIZeroJ;
1894
1895    one= bufFactors [l + 1];
1896    two= Pi [l - 1];
1897    if (degBuf > 0 && degPi > 0)
1898    {
1899      while (one.hasTerms() && one.exp() > j) one++;
1900      while (two.hasTerms() && two.exp() > j) two++;
1901      for (k= 1; k <= (int) ceil (j/2.0); k++)
1902      {
1903        if (k != j - k + 1)
1904        {
1905          if ((one.hasTerms() && one.exp() == j - k + 1) &&
1906              (two.hasTerms() && two.exp() == j - k + 1))
1907          {
1908            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
1909                      (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1) -
1910                      M (j - k + 2, l + 1);
1911            one++;
1912            two++;
1913          }
1914          else if (one.hasTerms() && one.exp() == j - k + 1)
1915          {
1916            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
1917                               Pi[l - 1] [k]) - M (k + 1, l + 1);
1918            one++;
1919          }
1920          else if (two.hasTerms() && two.exp() == j - k + 1)
1921          {
1922            tmp[l] += mulNTL (bufFactors[l + 1] [k],
1923                      (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1);
1924            two++;
1925           }
1926        }
1927        else
1928          tmp[l] += M (k + 1, l + 1);
1929      }
1930    }
1931
1932    if (degPi >= j + 1 && degBuf >= j + 1)
1933    {
1934      if (j + 2 <= M.rows())
1935        tmp [l] += mulNTL ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
1936                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
1937                          ) - M(1,l+1) - M (j + 2,l+1);
1938    }
1939    else if (degPi >= j + 1)
1940    {
1941      if (degBuf > 0)
1942        tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1] [0]);
1943      else
1944        tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1]);
1945    }
1946    else if (degBuf >= j + 1)
1947    {
1948      if (degPi > 0)
1949        tmp[l] += mulNTL (Pi [l - 1] [0], bufFactors [l + 1] [j + 1]);
1950      else
1951        tmp[l] += mulNTL (Pi [l - 1], bufFactors [l + 1] [j + 1]);
1952    }
1953
1954    Pi[l] += tmp[l]*xToJ*F.mvar();
1955  }
1956  return;
1957}
1958
1959void
1960nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l,
1961                      CFArray& Pi, CFList& diophant, CFMatrix& M,
1962                      const CFArray& LCs, bool sort)
1963{
1964  if (sort)
1965    sortList (factors, Variable (1));
1966  Pi= CFArray (factors.length() - 2);
1967  CFList bufFactors2= factors;
1968  bufFactors2.removeFirst();
1969  diophant= diophantine (F[0], bufFactors2);
1970  DEBOUTLN (cerr, "diophant= " << diophant);
1971
1972  CFArray bufFactors= CFArray (bufFactors2.length());
1973  int i= 0;
1974  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
1975    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
1976
1977  Variable x= F.mvar();
1978  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
1979  {
1980    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
1981    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
1982                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
1983  }
1984  else if (degree (bufFactors[0], x) > 0)
1985  {
1986    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
1987    Pi [0]= M (1, 1) +
1988            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
1989  }
1990  else if (degree (bufFactors[1], x) > 0)
1991  {
1992    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
1993    Pi [0]= M (1, 1) +
1994            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
1995  }
1996  else
1997  {
1998    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
1999    Pi [0]= M (1, 1);
2000  }
2001
2002  for (i= 1; i < Pi.size(); i++)
2003  {
2004    if (degree (Pi[i-1], x) > 0 && degree (bufFactors [i+1], x) > 0)
2005    {
2006      M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors[i+1] [0]);
2007      Pi [i]= M (1,i+1) + (mulNTL (Pi[i-1] [1], bufFactors[i+1] [0]) +
2008                       mulNTL (Pi[i-1] [0], bufFactors [i+1] [1]))*x;
2009    }
2010    else if (degree (Pi[i-1], x) > 0)
2011    {
2012      M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors [i+1]);
2013      Pi [i]=  M(1,i+1) + mulNTL (Pi[i-1] [1], bufFactors[i+1])*x;
2014    }
2015    else if (degree (bufFactors[i+1], x) > 0)
2016    {
2017      M (1,i+1)= mulNTL (Pi[i-1], bufFactors [i+1] [0]);
2018      Pi [i]= M (1,i+1) + mulNTL (Pi[i-1], bufFactors[i+1] [1])*x;
2019    }
2020    else
2021    {
2022      M (1,i+1)= mulNTL (Pi [i-1], bufFactors [i+1]);
2023      Pi [i]= M (1,i+1);
2024    }
2025  }
2026
2027  for (i= 1; i < l; i++)
2028    nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
2029
2030  factors= CFList();
2031  for (i= 0; i < bufFactors.size(); i++)
2032    factors.append (bufFactors[i]);
2033  return;
2034}
2035
2036
2037/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
2038/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
2039CFList
2040diophantine (const CFList& recResult, const CFList& factors,
2041             const CFList& products, const CFList& M, const CanonicalForm& E,
2042             bool& bad)
2043{
2044  if (M.isEmpty())
2045  {
2046    CFList result;
2047    CFListIterator j= factors;
2048    CanonicalForm buf;
2049    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2050    {
2051      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
2052              "constant or univariate poly expected");
2053      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
2054              "constant or univariate poly expected");
2055      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
2056              "constant or univariate poly expected");
2057      buf= mulNTL (E, i.getItem());
2058      result.append (modNTL (buf, j.getItem()));
2059    }
2060    return result;
2061  }
2062  Variable y= M.getLast().mvar();
2063  CFList bufFactors= factors;
2064  for (CFListIterator i= bufFactors; i.hasItem(); i++)
2065    i.getItem()= mod (i.getItem(), y);
2066  CFList bufProducts= products;
2067  for (CFListIterator i= bufProducts; i.hasItem(); i++)
2068    i.getItem()= mod (i.getItem(), y);
2069  CFList buf= M;
2070  buf.removeLast();
2071  CanonicalForm bufE= mod (E, y);
2072  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2073                                      bufE, bad);
2074
2075  if (bad)
2076    return CFList();
2077
2078  CanonicalForm e= E;
2079  CFListIterator j= products;
2080  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
2081    e -= i.getItem()*j.getItem();
2082
2083  CFList result= recDiophantine;
2084  int d= degree (M.getLast());
2085  CanonicalForm coeffE;
2086  for (int i= 1; i < d; i++)
2087  {
2088    if (degree (e, y) > 0)
2089      coeffE= e[i];
2090    else
2091      coeffE= 0;
2092    if (!coeffE.isZero())
2093    {
2094      CFListIterator k= result;
2095      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2096                                   coeffE, bad);
2097      if (bad)
2098        return CFList();
2099      CFListIterator l= products;
2100      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
2101      {
2102        k.getItem() += j.getItem()*power (y, i);
2103        e -= j.getItem()*power (y, i)*l.getItem();
2104      }
2105    }
2106    if (e.isZero())
2107      break;
2108  }
2109  if (!e.isZero())
2110  {
2111    bad= true;
2112    return CFList();
2113  }
2114  return result;
2115}
2116
2117void
2118nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
2119                    CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
2120                    CFArray& Pi, const CFList& products, int j,
2121                    const CFList& MOD, bool& noOneToOne)
2122{
2123  CanonicalForm E;
2124  CanonicalForm xToJ= power (F.mvar(), j);
2125  Variable x= F.mvar();
2126
2127  // compute the error
2128#ifdef DEBUGOUTPUT
2129    CanonicalForm test= 1;
2130    for (int i= 0; i < factors.length(); i++)
2131    {
2132      if (i == 0)
2133        test *= mod (bufFactors [i], power (x, j));
2134      else
2135        test *= bufFactors[i];
2136    }
2137    test= mod (test, power (x, j));
2138    test= mod (test, MOD);
2139    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
2140    DEBOUTLN (cerr, "test in nonMonicHenselStep= " << test2);
2141#endif
2142
2143  if (degree (Pi [factors.length() - 2], x) > 0)
2144    E= F[j] - Pi [factors.length() - 2] [j];
2145  else
2146    E= F[j];
2147
2148  CFArray buf= CFArray (diophant.length());
2149
2150  // actual lifting
2151  CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
2152                                    noOneToOne);
2153
2154  if (noOneToOne)
2155    return;
2156
2157  int k= 0;
2158  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
2159  {
2160    buf[k]= i.getItem();
2161    bufFactors[k] += xToJ*i.getItem();
2162  }
2163
2164  // update Pi [0]
2165  int degBuf0= degree (bufFactors[0], x);
2166  int degBuf1= degree (bufFactors[1], x);
2167  if (degBuf0 > 0 && degBuf1 > 0)
2168  {
2169    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
2170    if (j + 2 <= M.rows())
2171      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
2172  }
2173  CanonicalForm uIZeroJ;
2174
2175  if (degBuf0 > 0 && degBuf1 > 0)
2176    uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
2177             mulMod (bufFactors[1] [0], buf[0], MOD);
2178  else if (degBuf0 > 0)
2179    uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
2180  else if (degBuf1 > 0)
2181    uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
2182  else
2183    uIZeroJ= 0;
2184  Pi [0] += xToJ*uIZeroJ;
2185
2186  CFArray tmp= CFArray (factors.length() - 1);
2187  for (k= 0; k < factors.length() - 1; k++)
2188    tmp[k]= 0;
2189  CFIterator one, two;
2190  one= bufFactors [0];
2191  two= bufFactors [1];
2192  if (degBuf0 > 0 && degBuf1 > 0)
2193  {
2194    while (one.hasTerms() && one.exp() > j) one++;
2195    while (two.hasTerms() && two.exp() > j) two++;
2196    for (k= 1; k <= (int) ceil (j/2.0); k++)
2197    {
2198      if (k != j - k + 1)
2199      {
2200        if ((one.hasTerms() && one.exp() == j - k + 1) &&
2201            (two.hasTerms() && two.exp() == j - k + 1))
2202        {
2203          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2204                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
2205                    M (j - k + 2, 1);
2206          one++;
2207          two++;
2208        }
2209        else if (one.hasTerms() && one.exp() == j - k + 1)
2210        {
2211          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2212                    bufFactors[1] [k], MOD) - M (k + 1, 1);
2213          one++;
2214        }
2215        else if (two.hasTerms() && two.exp() == j - k + 1)
2216        {
2217          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
2218                    two.coeff()), MOD) - M (k + 1, 1);
2219          two++;
2220        }
2221      }
2222      else
2223      {
2224        tmp[0] += M (k + 1, 1);
2225      }
2226    }
2227  }
2228
2229  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
2230  {
2231    if (j + 2 <= M.rows())
2232      tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
2233                         (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
2234                         - M(1,1) - M (j + 2,1);
2235  }
2236  else if (degBuf0 >= j + 1)
2237  {
2238    if (degBuf1 > 0)
2239      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
2240    else
2241      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
2242  }
2243  else if (degBuf1 >= j + 1)
2244  {
2245    if (degBuf0 > 0)
2246      tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
2247    else
2248      tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
2249  }
2250  Pi [0] += tmp[0]*xToJ*F.mvar();
2251
2252  // update Pi [l]
2253  int degPi, degBuf;
2254  for (int l= 1; l < factors.length() - 1; l++)
2255  {
2256    degPi= degree (Pi [l - 1], x);
2257    degBuf= degree (bufFactors[l + 1], x);
2258    if (degPi > 0 && degBuf > 0)
2259    {
2260      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
2261      if (j + 2 <= M.rows())
2262        M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
2263                                  MOD);
2264    }
2265
2266    if (degPi > 0 && degBuf > 0)
2267      uIZeroJ= mulMod (Pi[l -1] [0], buf[l + 1], MOD) +
2268               mulMod (uIZeroJ, bufFactors[l+1] [0], MOD);
2269    else if (degPi > 0)
2270      uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD);
2271    else if (degBuf > 0)
2272      uIZeroJ= mulMod (Pi[l - 1], buf[1], MOD);
2273    else
2274      uIZeroJ= 0;
2275
2276    Pi [l] += xToJ*uIZeroJ;
2277
2278    one= bufFactors [l + 1];
2279    two= Pi [l - 1];
2280    if (degBuf > 0 && degPi > 0)
2281    {
2282      while (one.hasTerms() && one.exp() > j) one++;
2283      while (two.hasTerms() && two.exp() > j) two++;
2284      for (k= 1; k <= (int) ceil (j/2.0); k++)
2285      {
2286        if (k != j - k + 1)
2287        {
2288          if ((one.hasTerms() && one.exp() == j - k + 1) &&
2289              (two.hasTerms() && two.exp() == j - k + 1))
2290          {
2291            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2292                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
2293                      M (j - k + 2, l + 1);
2294            one++;
2295            two++;
2296          }
2297          else if (one.hasTerms() && one.exp() == j - k + 1)
2298          {
2299            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2300                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
2301            one++;
2302          }
2303          else if (two.hasTerms() && two.exp() == j - k + 1)
2304          {
2305            tmp[l] += mulMod (bufFactors[l + 1] [k],
2306                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
2307            two++;
2308           }
2309        }
2310        else
2311          tmp[l] += M (k + 1, l + 1);
2312      }
2313    }
2314
2315    if (degPi >= j + 1 && degBuf >= j + 1)
2316    {
2317      if (j + 2 <= M.rows())
2318        tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
2319                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
2320                            , MOD) - M(1,l+1) - M (j + 2,l+1);
2321    }
2322    else if (degPi >= j + 1)
2323    {
2324      if (degBuf > 0)
2325        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
2326      else
2327        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
2328    }
2329    else if (degBuf >= j + 1)
2330    {
2331      if (degPi > 0)
2332        tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
2333      else
2334        tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
2335    }
2336
2337    Pi[l] += tmp[l]*xToJ*F.mvar();
2338  }
2339  return;
2340}
2341
2342// wrt. Variable (1)
2343CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
2344{
2345  if (degree (F, 1) <= 0)
2346   return c;
2347  else
2348  {
2349    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
2350    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
2351              - LC (result))*power (result.mvar(), degree (result));
2352    return swapvar (result, Variable (F.level() + 1), Variable (1));
2353  }
2354}
2355
2356CFList
2357nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
2358                      diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
2359                      const CFList& LCs2, bool& bad)
2360{
2361  CFList buf= factors;
2362  int k= 0;
2363  int liftBoundBivar= l[k];
2364  CFList bufbuf= factors;
2365  Variable v= Variable (2);
2366
2367  CFList MOD;
2368  MOD.append (power (Variable (2), liftBoundBivar));
2369  CFArray bufFactors= CFArray (factors.length());
2370  k= 0;
2371  CFListIterator j= eval;
2372  j++;
2373  CFListIterator iter1= LCs1;
2374  CFListIterator iter2= LCs2;
2375  iter1++;
2376  iter2++;
2377  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
2378  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
2379
2380  CFListIterator i= buf;
2381  i++;
2382  Variable y= j.getItem().mvar();
2383  if (y.level() != 3)
2384    y= Variable (3);
2385
2386  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
2387  M (1, 1)= Pi[0];
2388  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2389    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2390                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2391  else if (degree (bufFactors[0], y) > 0)
2392    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2393  else if (degree (bufFactors[1], y) > 0)
2394    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2395
2396  CFList products;
2397  for (int i= 0; i < bufFactors.size(); i++)
2398  {
2399    if (degree (bufFactors[i], y) > 0)
2400      products.append (eval.getFirst()/bufFactors[i] [0]);
2401    else
2402      products.append (eval.getFirst()/bufFactors[i]);
2403  }
2404
2405  for (int d= 1; d < l[1]; d++)
2406  {
2407    nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
2408                        d, MOD, bad);
2409    if (bad)
2410      return CFList();
2411  }
2412  CFList result;
2413  for (k= 0; k < factors.length(); k++)
2414    result.append (bufFactors[k]);
2415  return result;
2416}
2417
2418
2419CFList
2420nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
2421                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2422                    int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
2423                   )
2424{
2425  int k= 0;
2426  CFArray bufFactors= CFArray (factors.length());
2427  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
2428  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
2429  CFList buf= factors;
2430  Variable y= F.getLast().mvar();
2431  Variable x= F.getFirst().mvar();
2432  CanonicalForm xToLOld= power (x, lOld);
2433  Pi [0]= mod (Pi[0], xToLOld);
2434  M (1, 1)= Pi [0];
2435
2436  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2437    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2438                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2439  else if (degree (bufFactors[0], y) > 0)
2440    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2441  else if (degree (bufFactors[1], y) > 0)
2442    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2443
2444  CFList products;
2445  CanonicalForm quot;
2446  for (int i= 0; i < bufFactors.size(); i++)
2447  {
2448    if (degree (bufFactors[i], y) > 0)
2449    {
2450      if (!fdivides (bufFactors[i] [0], F.getFirst(), quot))
2451      {
2452        bad= true;
2453        return CFList();
2454      }
2455      products.append (quot);
2456    }
2457    else
2458    {
2459      if (!fdivides (bufFactors[i], F.getFirst(), quot))
2460      {
2461        bad= true;
2462        return CFList();
2463      }
2464      products.append (quot);
2465    }
2466  }
2467
2468  for (int d= 1; d < lNew; d++)
2469  {
2470    nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
2471                        d, MOD, bad);
2472    if (bad)
2473      return CFList();
2474  }
2475
2476  CFList result;
2477  for (k= 0; k < factors.length(); k++)
2478    result.append (bufFactors[k]);
2479  return result;
2480}
2481
2482CFList
2483nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
2484                    lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
2485                    const CFArray& Pi, const CFList& diophant, bool& bad)
2486{
2487  CFList bufDiophant= diophant;
2488  CFList buf= factors;
2489  if (sort)
2490    sortList (buf, Variable (1));
2491  CFArray bufPi= Pi;
2492  CFMatrix M= CFMatrix (l[1], factors.length());
2493  CFList result= 
2494    nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
2495  if (bad)
2496    return CFList();
2497
2498  if (eval.length() == 2)
2499    return result;
2500  CFList MOD;
2501  for (int i= 0; i < 2; i++)
2502    MOD.append (power (Variable (i + 2), l[i]));
2503  CFListIterator j= eval;
2504  j++;
2505  CFList bufEval;
2506  bufEval.append (j.getItem());
2507  j++;
2508  CFListIterator jj= LCs1;
2509  CFListIterator jjj= LCs2;
2510  CFList bufLCs1, bufLCs2;
2511  jj++, jjj++;
2512  bufLCs1.append (jj.getItem());
2513  bufLCs2.append (jjj.getItem());
2514  jj++, jjj++;
2515
2516  for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
2517  {
2518    bufEval.append (j.getItem());
2519    bufLCs1.append (jj.getItem());
2520    bufLCs2.append (jjj.getItem());
2521    M= CFMatrix (l[i], factors.length());
2522    result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
2523                                 l[i - 1], l[i], bufLCs1, bufLCs2, bad);
2524    if (bad)
2525      return CFList();
2526    MOD.append (power (Variable (i + 2), l[i]));
2527    bufEval.removeFirst();
2528    bufLCs1.removeFirst();
2529    bufLCs2.removeFirst();
2530  }
2531  return result;
2532}
2533
2534CFList
2535nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
2536                      CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
2537                      int bivarLiftBound, bool& bad)
2538{
2539  CFList bufFactors2= factors;
2540
2541  Variable y= Variable (2);
2542  for (CFListIterator i= bufFactors2; i.hasItem(); i++)
2543    i.getItem()= mod (i.getItem(), y);
2544
2545  CanonicalForm bufF= F;
2546  bufF= mod (bufF, y);
2547  bufF= mod (bufF, Variable (3));
2548
2549  diophant= diophantine (bufF, bufFactors2);
2550
2551  CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
2552
2553  Pi= CFArray (bufFactors2.length() - 1);
2554
2555  CFArray bufFactors= CFArray (bufFactors2.length());
2556  CFListIterator j= LCs;
2557  int i= 0;
2558  for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
2559    bufFactors[i]= replaceLC (k.getItem(), j.getItem());
2560
2561  //initialise Pi
2562  Variable v= Variable (3);
2563  CanonicalForm yToL= power (y, bivarLiftBound);
2564  if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
2565  {
2566    M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
2567    Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
2568                       mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
2569  }
2570  else if (degree (bufFactors[0], v) > 0)
2571  {
2572    M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
2573    Pi [0]=  M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
2574  }
2575  else if (degree (bufFactors[1], v) > 0)
2576  {
2577    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
2578    Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
2579  }
2580  else
2581  {
2582    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
2583    Pi [0]= M (1,1);
2584  }
2585
2586  for (i= 1; i < Pi.size(); i++)
2587  {
2588    if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
2589    {
2590      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
2591      Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
2592                       mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
2593    }
2594    else if (degree (Pi[i-1], v) > 0)
2595    {
2596      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
2597      Pi [i]=  M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
2598    }
2599    else if (degree (bufFactors[i+1], v) > 0)
2600    {
2601      M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
2602      Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
2603    }
2604    else
2605    {
2606      M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
2607      Pi [i]= M (1,i+1);
2608    }
2609  }
2610
2611  CFList products;
2612  bufF= mod (F, Variable (3));
2613  for (CFListIterator k= factors; k.hasItem(); k++)
2614    products.append (bufF/k.getItem());
2615
2616  CFList MOD= CFList (power (v, liftBound));
2617  MOD.insert (yToL);
2618  for (int d= 1; d < liftBound; d++)
2619  {
2620    nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
2621                        MOD, bad);
2622    if (bad)
2623      return CFList();
2624  }
2625
2626  CFList result;
2627  for (i= 0; i < factors.length(); i++)
2628    result.append (bufFactors[i]);
2629  return result;
2630}
2631
2632CFList
2633nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
2634                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2635                    int& lNew, const CFList& MOD, bool& noOneToOne
2636                   )
2637{
2638
2639  int k= 0;
2640  CFArray bufFactors= CFArray (factors.length());
2641  CFListIterator j= LCs;
2642  for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
2643    bufFactors [k]= replaceLC (i.getItem(), j.getItem());
2644
2645  Variable y= F.getLast().mvar();
2646  Variable x= F.getFirst().mvar();
2647  CanonicalForm xToLOld= power (x, lOld);
2648
2649  Pi [0]= mod (Pi[0], xToLOld);
2650  M (1, 1)= Pi [0];
2651
2652  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2653    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2654                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2655  else if (degree (bufFactors[0], y) > 0)
2656    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2657  else if (degree (bufFactors[1], y) > 0)
2658    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2659
2660  for (int i= 1; i < Pi.size(); i++)
2661  {
2662    Pi [i]= mod (Pi [i], xToLOld);
2663    M (1, i + 1)= Pi [i];
2664
2665    if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
2666      Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
2667                 mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
2668    else if (degree (Pi[i-1], y) > 0)
2669      Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
2670    else if (degree (bufFactors[i+1], y) > 0)
2671      Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
2672  }
2673
2674  CFList products;
2675  CanonicalForm quot, bufF= F.getFirst();
2676
2677  for (int i= 0; i < bufFactors.size(); i++)
2678  {
2679    if (degree (bufFactors[i], y) > 0)
2680    {
2681      if (!fdivides (bufFactors[i] [0], bufF, quot))
2682      {
2683        noOneToOne= true;
2684        return factors;
2685      }
2686      products.append (quot);
2687    }
2688    else
2689    {
2690      if (!fdivides (bufFactors[i], bufF, quot))
2691      {
2692        noOneToOne= true;
2693        return factors;
2694      }
2695      products.append (quot);
2696    }
2697  }
2698
2699  for (int d= 1; d < lNew; d++)
2700  {
2701    nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
2702                        products, d, MOD, noOneToOne);
2703    if (noOneToOne)
2704      return CFList();
2705  }
2706
2707  CFList result;
2708  for (k= 0; k < factors.length(); k++)
2709    result.append (bufFactors[k]);
2710  return result;
2711}
2712
2713CFList
2714nonMonicHenselLift (const CFList& eval, const CFList& factors,
2715                    CFList* const& LCs, CFList& diophant, CFArray& Pi,
2716                    int* liftBound, int length, bool& noOneToOne
2717                   )
2718{
2719  CFList bufDiophant= diophant;
2720  CFList buf= factors;
2721  CFArray bufPi= Pi;
2722  CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
2723  int k= 0;
2724
2725  CFList result=
2726  nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
2727                        liftBound[1], liftBound[0], noOneToOne);
2728
2729  if (noOneToOne)
2730    return CFList();
2731
2732  if (eval.length() == 1)
2733    return result;
2734
2735  k++;
2736  CFList MOD;
2737  for (int i= 0; i < 2; i++)
2738    MOD.append (power (Variable (i + 2), liftBound[i]));
2739
2740  CFListIterator j= eval;
2741  CFList bufEval;
2742  bufEval.append (j.getItem());
2743  j++;
2744
2745  for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
2746  {
2747    bufEval.append (j.getItem());
2748    M= CFMatrix (liftBound[i], factors.length() - 1);
2749    result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
2750                                liftBound[i-1], liftBound[i], MOD, noOneToOne);
2751    if (noOneToOne)
2752      return result;
2753    MOD.append (power (Variable (i + 2), liftBound[i]));
2754    bufEval.removeFirst();
2755  }
2756
2757  return result;
2758}
2759
2760#endif
2761/* HAVE_NTL */
2762
Note: See TracBrowser for help on using the repository browser.