source: git/factory/facHensel.cc @ 49ba09

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