source: git/factory/facHensel.cc @ 81d96c

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