source: git/factory/facHensel.cc @ 73c222

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