source: git/factory/facHensel.cc @ 82e0a7

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