source: git/factory/facHensel.cc @ 72bfc8

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