source: git/factory/facHensel.cc @ c9aa52

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