source: git/factory/facHensel.cc @ e016ba

fieker-DuValspielwiese
Last change on this file since e016ba was e016ba, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: code formatting
  • 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)
493             && (two.hasTerms() && two.exp() == j - k + 1))
494        {
495          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), (bufFactors[1][k]+
496                     two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
497          one++;
498          two++;
499        }
500        else if (one.hasTerms() && one.exp() == j - k + 1)
501        {
502          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1][k])-
503                     M (k + 1, 1);
504          one++;
505        }
506        else if (two.hasTerms() && two.exp() == j - k + 1)
507        {
508          tmp[0] += mulNTL (bufFactors[0][k], (bufFactors[1][k]+two.coeff()))-
509                    M (k + 1, 1);
510          two++;
511        }
512      }
513      else
514      {
515        tmp[0] += M (k + 1, 1);
516      }
517    }
518  }
519  Pi [0] += tmp[0]*xToJ*F.mvar();
520
521  // update Pi [l]
522  int degPi, degBuf;
523  for (int l= 1; l < factors.length() - 1; l++)
524  {
525    degPi= degree (Pi [l - 1], x);
526    degBuf= degree (bufFactors[l + 1], x);
527    if (degPi > 0 && degBuf > 0)
528      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
529    if (j == 1)
530    {
531      if (degPi > 0 && degBuf > 0)
532        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
533                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
534                  M (1, l + 1));
535      else if (degPi > 0)
536        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
537      else if (degBuf > 0)
538        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
539    }
540    else
541    {
542      if (degPi > 0 && degBuf > 0)
543      {
544        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
545        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
546      }
547      else if (degPi > 0)
548        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
549      else if (degBuf > 0)
550      {
551        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
552        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
553      }
554      Pi[l] += xToJ*uIZeroJ;
555    }
556    one= bufFactors [l + 1];
557    two= Pi [l - 1];
558    if (two.hasTerms() && two.exp() == j + 1)
559    {
560      if (degBuf > 0 && degPi > 0)
561      {
562          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
563          two++;
564      }
565      else if (degPi > 0)
566      {
567          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
568          two++;
569      }
570    }
571    if (degBuf > 0 && degPi > 0)
572    {
573      for (k= 1; k <= (int) ceil (j/2.0); k++)
574      {
575        if (k != j - k + 1)
576        {
577          if ((one.hasTerms() && one.exp() == j - k + 1) &&
578              (two.hasTerms() && two.exp() == j - k + 1))
579          {
580            tmp[l] += mulNTL ((bufFactors[l+1][k] + one.coeff()), (Pi[l-1][k] +
581                       two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
582            one++;
583            two++;
584          }
585          else if (one.hasTerms() && one.exp() == j - k + 1)
586          {
587            tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()), Pi[l-1][k]) -
588                       M (k + 1, l + 1);
589            one++;
590          }
591          else if (two.hasTerms() && two.exp() == j - k + 1)
592          {
593            tmp[l] += mulNTL (bufFactors[l+1][k], (Pi[l-1][k] + two.coeff())) -
594                      M (k + 1, l + 1);
595            two++;
596          }
597        }
598        else
599          tmp[l] += M (k + 1, l + 1);
600      }
601    }
602    Pi[l] += tmp[l]*xToJ*F.mvar();
603  }
604  return;
605}
606
607void
608henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
609              CFList& diophant, CFMatrix& M, bool sort)
610{
611  if (sort)
612    sortList (factors, Variable (1));
613  Pi= CFArray (factors.length() - 1);
614  CFListIterator j= factors;
615  diophant= diophantine (F[0], factors);
616  DEBOUTLN (cerr, "diophant= " << diophant);
617  j++;
618  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()));
619  M (1, 1)= Pi [0];
620  int i= 1;
621  if (j.hasItem())
622    j++;
623  for (; j.hasItem(); j++, i++)
624  {
625    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
626    M (1, i + 1)= Pi [i];
627  }
628  CFArray bufFactors= CFArray (factors.length());
629  i= 0;
630  for (CFListIterator k= factors; k.hasItem(); i++, k++)
631  {
632    if (i == 0)
633      bufFactors[i]= mod (k.getItem(), F.mvar());
634    else
635      bufFactors[i]= k.getItem();
636  }
637  for (i= 1; i < l; i++)
638    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
639
640  CFListIterator k= factors;
641  for (i= 0; i < factors.length (); i++, k++)
642    k.getItem()= bufFactors[i];
643  factors.removeFirst();
644  return;
645}
646
647void
648henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
649                    end, CFArray& Pi, const CFList& diophant, CFMatrix& M)
650{
651  CFArray bufFactors= CFArray (factors.length());
652  int i= 0;
653  CanonicalForm xToStart= power (F.mvar(), start);
654  for (CFListIterator k= factors; k.hasItem(); k++, i++)
655  {
656    if (i == 0)
657      bufFactors[i]= mod (k.getItem(), xToStart);
658    else
659      bufFactors[i]= k.getItem();
660  }
661  for (i= start; i < end; i++)
662    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
663
664  CFListIterator k= factors;
665  for (i= 0; i < factors.length(); k++, i++)
666    k.getItem()= bufFactors [i];
667  factors.removeFirst();
668  return;
669}
670
671CFList
672biDiophantine (const CanonicalForm& F, const CFList& factors, int d)
673{
674  Variable y= F.mvar();
675  CFList result;
676  if (y.level() == 1)
677  {
678    result= diophantine (F, factors);
679    return result;
680  }
681  else
682  {
683    CFList buf= factors;
684    for (CFListIterator i= buf; i.hasItem(); i++)
685      i.getItem()= mod (i.getItem(), y);
686    CanonicalForm A= mod (F, y);
687    int bufD= 1;
688    CFList recResult= biDiophantine (A, buf, bufD);
689    CanonicalForm e= 1;
690    CFList p;
691    CFArray bufFactors= CFArray (factors.length());
692    CanonicalForm yToD= power (y, d);
693    int k= 0;
694    for (CFListIterator i= factors; i.hasItem(); i++, k++)
695    {
696      bufFactors [k]= i.getItem();
697    }
698    CanonicalForm b, quot;
699    for (k= 0; k < factors.length(); k++) //TODO compute b's faster
700    {
701      b= 1;
702      if (fdivides (bufFactors[k], F, quot))
703        b= quot;
704      else
705      {
706        for (int l= 0; l < factors.length(); l++)
707        {
708          if (l == k)
709            continue;
710          else
711          {
712            b= mulMod2 (b, bufFactors[l], yToD);
713          }
714        }
715      }
716      p.append (b);
717    }
718
719    CFListIterator j= p;
720    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
721      e -= i.getItem()*j.getItem();
722
723    if (e.isZero())
724      return recResult;
725    CanonicalForm coeffE;
726    CFList s;
727    result= recResult;
728    CanonicalForm g;
729    for (int i= 1; i < d; i++)
730    {
731      if (degree (e, y) > 0)
732        coeffE= e[i];
733      else
734        coeffE= 0;
735      if (!coeffE.isZero())
736      {
737        CFListIterator k= result;
738        CFListIterator l= p;
739        int ii= 0;
740        j= recResult;
741        for (; j.hasItem(); j++, k++, l++, ii++)
742        {
743          g= coeffE*j.getItem();
744          if (degree (bufFactors[ii], y) <= 0)
745            g= mod (g, bufFactors[ii]);
746          else
747            g= mod (g, bufFactors[ii][0]);
748          k.getItem() += g*power (y, i);
749          e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
750          DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
751                    mod (e, power (y, i + 1)));
752        }
753      }
754      if (e.isZero())
755        break;
756    }
757
758    DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
759
760#ifdef DEBUGOUTPUT
761    CanonicalForm test= 0;
762    j= p;
763    for (CFListIterator i= result; i.hasItem(); i++, j++)
764      test += mod (i.getItem()*j.getItem(), power (y, d));
765    DEBOUTLN (cerr, "test= " << test);
766#endif
767    return result;
768  }
769}
770
771CFList
772multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
773                     const CFList& recResult, const CFList& M, int d)
774{
775  Variable y= F.mvar();
776  CFList result;
777  CFListIterator i;
778  CanonicalForm e= 1;
779  CFListIterator j= factors;
780  CFList p;
781  CFArray bufFactors= CFArray (factors.length());
782  CanonicalForm yToD= power (y, d);
783  int k= 0;
784  for (CFListIterator i= factors; i.hasItem(); i++, k++)
785    bufFactors [k]= i.getItem();
786  CanonicalForm b, quot;
787  CFList buf= M;
788  buf.removeLast();
789  buf.append (yToD);
790  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
791  {
792    b= 1;
793    if (fdivides (bufFactors[k], F, quot))
794      b= quot;
795    else
796    {
797      for (int l= 0; l < factors.length(); l++)
798      {
799        if (l == k)
800          continue;
801        else
802        {
803          b= mulMod (b, bufFactors[l], buf);
804        }
805      }
806    }
807    p.append (b);
808  }
809  j= p;
810  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
811    e -= mulMod (i.getItem(), j.getItem(), M);
812
813  if (e.isZero())
814    return recResult;
815  CanonicalForm coeffE;
816  CFList s;
817  result= recResult;
818  CanonicalForm g;
819  for (int i= 1; i < d; i++)
820  {
821    if (degree (e, y) > 0)
822      coeffE= e[i];
823    else
824      coeffE= 0;
825    if (!coeffE.isZero())
826    {
827      CFListIterator k= result;
828      CFListIterator l= p;
829      j= recResult;
830      int ii= 0;
831      CanonicalForm dummy;
832      for (; j.hasItem(); j++, k++, l++, ii++)
833      {
834        g= mulMod (coeffE, j.getItem(), M);
835        if (degree (bufFactors[ii], y) <= 0)
836          divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
837                  g, M);
838        else
839          divrem (g, bufFactors[ii][0], dummy, g, M);
840        k.getItem() += g*power (y, i);
841        e -= mulMod (g*power (y, i), l.getItem(), M);
842      }
843    }
844
845    if (e.isZero())
846      break;
847  }
848
849#ifdef DEBUGOUTPUT
850    CanonicalForm test= 0;
851    j= p;
852    for (CFListIterator i= result; i.hasItem(); i++, j++)
853      test += mod (i.getItem()*j.getItem(), power (y, d));
854    DEBOUTLN (cerr, "test= " << test);
855#endif
856  return result;
857}
858
859static
860void
861henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
862            const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
863            const CFList& MOD)
864{
865  CanonicalForm E;
866  CanonicalForm xToJ= power (F.mvar(), j);
867  Variable x= F.mvar();
868  // compute the error
869  if (j == 1)
870  {
871    E= F[j];
872#ifdef DEBUGOUTPUT
873    CanonicalForm test= 1;
874    for (int i= 0; i < factors.length(); i++)
875    {
876      if (i == 0)
877        test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
878      else
879        test= mulMod (test, bufFactors[i], MOD);
880    }
881    CanonicalForm test2= mod (F-test, xToJ);
882
883    test2= mod (test2, MOD);
884    DEBOUTLN (cerr, "test= " << test2);
885#endif
886  }
887  else
888  {
889#ifdef DEBUGOUTPUT
890    CanonicalForm test= 1;
891    for (int i= 0; i < factors.length(); i++)
892    {
893      if (i == 0)
894        test *= mod (bufFactors [i], power (x, j));
895      else
896        test *= bufFactors[i];
897    }
898    test= mod (test, power (x, j));
899    test= mod (test, MOD);
900    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
901    DEBOUTLN (cerr, "test= " << test2);
902#endif
903
904    if (degree (Pi [factors.length() - 2], x) > 0)
905      E= F[j] - Pi [factors.length() - 2] [j];
906    else
907      E= F[j];
908  }
909
910  CFArray buf= CFArray (diophant.length());
911  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
912  int k= 0;
913  // actual lifting
914  CanonicalForm dummy, rest1;
915  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
916  {
917    if (degree (bufFactors[k], x) > 0)
918    {
919      if (k > 0)
920        divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
921      else
922        rest1= E;
923    }
924    else
925      divrem (E, bufFactors[k], dummy, rest1, MOD);
926
927    buf[k]= mulMod (i.getItem(), rest1, MOD);
928
929    if (degree (bufFactors[k], x) > 0)
930      divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
931    else
932      divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
933  }
934  for (k= 1; k < factors.length(); k++)
935    bufFactors[k] += xToJ*buf[k];
936
937  // update Pi [0]
938  int degBuf0= degree (bufFactors[0], x);
939  int degBuf1= degree (bufFactors[1], x);
940  if (degBuf0 > 0 && degBuf1 > 0)
941    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
942  CanonicalForm uIZeroJ;
943  if (j == 1)
944  {
945    if (degBuf0 > 0 && degBuf1 > 0)
946      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
947                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
948    else if (degBuf0 > 0)
949      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
950    else if (degBuf1 > 0)
951      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
952    else
953      uIZeroJ= 0;
954    Pi [0] += xToJ*uIZeroJ;
955  }
956  else
957  {
958    if (degBuf0 > 0 && degBuf1 > 0)
959      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
960                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
961    else if (degBuf0 > 0)
962      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
963    else if (degBuf1 > 0)
964      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
965    else
966      uIZeroJ= 0;
967    Pi [0] += xToJ*uIZeroJ;
968  }
969
970  CFArray tmp= CFArray (factors.length() - 1);
971  for (k= 0; k < factors.length() - 1; k++)
972    tmp[k]= 0;
973  CFIterator one, two;
974  one= bufFactors [0];
975  two= bufFactors [1];
976  if (degBuf0 > 0 && degBuf1 > 0)
977  {
978    for (k= 1; k <= (int) ceil (j/2.0); k++)
979    {
980      if (k != j - k + 1)
981      {
982        if ((one.hasTerms() && one.exp() == j - k + 1) &&
983            (two.hasTerms() && two.exp() == j - k + 1))
984        {
985          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
986                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
987                    M (j - k + 2, 1);
988          one++;
989          two++;
990        }
991        else if (one.hasTerms() && one.exp() == j - k + 1)
992        {
993          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
994                    bufFactors[1] [k], MOD) - M (k + 1, 1);
995          one++;
996        }
997        else if (two.hasTerms() && two.exp() == j - k + 1)
998        {
999          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
1000                    two.coeff()), MOD) - M (k + 1, 1);
1001          two++;
1002        }
1003      }
1004      else
1005      {
1006        tmp[0] += M (k + 1, 1);
1007      }
1008    }
1009  }
1010  Pi [0] += tmp[0]*xToJ*F.mvar();
1011
1012  // update Pi [l]
1013  int degPi, degBuf;
1014  for (int l= 1; l < factors.length() - 1; l++)
1015  {
1016    degPi= degree (Pi [l - 1], x);
1017    degBuf= degree (bufFactors[l + 1], x);
1018    if (degPi > 0 && degBuf > 0)
1019      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
1020    if (j == 1)
1021    {
1022      if (degPi > 0 && degBuf > 0)
1023        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
1024                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
1025                  M (1, l + 1));
1026      else if (degPi > 0)
1027        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
1028      else if (degBuf > 0)
1029        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
1030    }
1031    else
1032    {
1033      if (degPi > 0 && degBuf > 0)
1034      {
1035        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
1036        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
1037      }
1038      else if (degPi > 0)
1039        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
1040      else if (degBuf > 0)
1041      {
1042        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
1043        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
1044      }
1045      Pi[l] += xToJ*uIZeroJ;
1046    }
1047    one= bufFactors [l + 1];
1048    two= Pi [l - 1];
1049    if (two.hasTerms() && two.exp() == j + 1)
1050    {
1051      if (degBuf > 0 && degPi > 0)
1052      {
1053          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
1054          two++;
1055      }
1056      else if (degPi > 0)
1057      {
1058          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
1059          two++;
1060      }
1061    }
1062    if (degBuf > 0 && degPi > 0)
1063    {
1064      for (k= 1; k <= (int) ceil (j/2.0); k++)
1065      {
1066        if (k != j - k + 1)
1067        {
1068          if ((one.hasTerms() && one.exp() == j - k + 1) &&
1069              (two.hasTerms() && two.exp() == j - k + 1))
1070          {
1071            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1072                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
1073                      M (j - k + 2, l + 1);
1074            one++;
1075            two++;
1076          }
1077          else if (one.hasTerms() && one.exp() == j - k + 1)
1078          {
1079            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1080                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
1081            one++;
1082          }
1083          else if (two.hasTerms() && two.exp() == j - k + 1)
1084          {
1085            tmp[l] += mulMod (bufFactors[l + 1] [k],
1086                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
1087            two++;
1088          }
1089        }
1090        else
1091          tmp[l] += M (k + 1, l + 1);
1092      }
1093    }
1094    Pi[l] += tmp[l]*xToJ*F.mvar();
1095  }
1096
1097  return;
1098}
1099
1100CFList
1101henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList&
1102              diophant, CFArray& Pi, CFMatrix& M)
1103{
1104  CFList buf= factors;
1105  int k= 0;
1106  int liftBoundBivar= l[k];
1107  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
1108  CFList MOD;
1109  MOD.append (power (Variable (2), liftBoundBivar));
1110  CFArray bufFactors= CFArray (factors.length());
1111  k= 0;
1112  CFListIterator j= eval;
1113  j++;
1114  buf.removeFirst();
1115  buf.insert (LC (j.getItem(), 1));
1116  for (CFListIterator i= buf; i.hasItem(); i++, k++)
1117    bufFactors[k]= i.getItem();
1118  Pi= CFArray (factors.length() - 1);
1119  CFListIterator i= buf;
1120  i++;
1121  Variable y= j.getItem().mvar();
1122  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
1123  M (1, 1)= Pi [0];
1124  k= 1;
1125  if (i.hasItem())
1126    i++;
1127  for (; i.hasItem(); i++, k++)
1128  {
1129    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
1130    M (1, k + 1)= Pi [k];
1131  }
1132
1133  for (int d= 1; d < l[1]; d++)
1134    henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
1135  CFList result;
1136  for (k= 1; k < factors.length(); k++)
1137    result.append (bufFactors[k]);
1138  return result;
1139}
1140
1141void
1142henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
1143                  CFArray& Pi, const CFList& diophant, CFMatrix& M,
1144                  const CFList& MOD)
1145{
1146  CFArray bufFactors= CFArray (factors.length());
1147  int i= 0;
1148  CanonicalForm xToStart= power (F.mvar(), start);
1149  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1150  {
1151    if (i == 0)
1152      bufFactors[i]= mod (k.getItem(), xToStart);
1153    else
1154      bufFactors[i]= k.getItem();
1155  }
1156  for (i= start; i < end; i++)
1157    henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD);
1158
1159  CFListIterator k= factors;
1160  for (i= 0; i < factors.length(); k++, i++)
1161    k.getItem()= bufFactors [i];
1162  factors.removeFirst();
1163  return;
1164}
1165
1166CFList
1167henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
1168            diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew)
1169{
1170  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
1171  int k= 0;
1172  CFArray bufFactors= CFArray (factors.length());
1173  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1174  {
1175    if (k == 0)
1176      bufFactors[k]= LC (F.getLast(), 1);
1177    else
1178      bufFactors[k]= i.getItem();
1179  }
1180  CFList buf= factors;
1181  buf.removeFirst();
1182  buf.insert (LC (F.getLast(), 1));
1183  CFListIterator i= buf;
1184  i++;
1185  Variable y= F.getLast().mvar();
1186  Variable x= F.getFirst().mvar();
1187  CanonicalForm xToLOld= power (x, lOld);
1188  Pi [0]= mod (Pi[0], xToLOld);
1189  M (1, 1)= Pi [0];
1190  k= 1;
1191  if (i.hasItem())
1192    i++;
1193  for (; i.hasItem(); i++, k++)
1194  {
1195    Pi [k]= mod (Pi [k], xToLOld);
1196    M (1, k + 1)= Pi [k];
1197  }
1198
1199  for (int d= 1; d < lNew; d++)
1200    henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
1201  CFList result;
1202  for (k= 1; k < factors.length(); k++)
1203    result.append (bufFactors[k]);
1204  return result;
1205}
1206
1207CFList
1208henselLift (const CFList& eval, const CFList& factors, int* l, int lLength,
1209            bool sort)
1210{
1211  CFList diophant;
1212  CFList buf= factors;
1213  buf.insert (LC (eval.getFirst(), 1));
1214  if (sort)
1215    sortList (buf, Variable (1));
1216  CFArray Pi;
1217  CFMatrix M= CFMatrix (l[1], factors.length());
1218  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
1219  if (eval.length() == 2)
1220    return result;
1221  CFList MOD;
1222  for (int i= 0; i < 2; i++)
1223    MOD.append (power (Variable (i + 2), l[i]));
1224  CFListIterator j= eval;
1225  j++;
1226  CFList bufEval;
1227  bufEval.append (j.getItem());
1228  j++;
1229
1230  for (int i= 2; i < lLength && j.hasItem(); i++, j++)
1231  {
1232    result.insert (LC (bufEval.getFirst(), 1));
1233    bufEval.append (j.getItem());
1234    M= CFMatrix (l[i], factors.length());
1235    result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
1236    MOD.append (power (Variable (i + 2), l[i]));
1237    bufEval.removeFirst();
1238  }
1239  return result;
1240}
1241
1242// nonmonic
1243
1244void
1245nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors,
1246                      CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1247                      CFArray& Pi, int j, const CFArray& /*LCs*/)
1248{
1249  Variable x= F.mvar();
1250  CanonicalForm xToJ= power (x, j);
1251
1252  CanonicalForm E;
1253  // compute the error
1254  if (degree (Pi [factors.length() - 2], x) > 0)
1255    E= F[j] - Pi [factors.length() - 2] [j];
1256  else
1257    E= F[j];
1258
1259  CFArray buf= CFArray (diophant.length());
1260
1261  int k= 0;
1262  CanonicalForm remainder;
1263  // actual lifting
1264  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1265  {
1266    if (degree (bufFactors[k], x) > 0)
1267      remainder= modNTL (E, bufFactors[k] [0]);
1268    else
1269      remainder= modNTL (E, bufFactors[k]);
1270    buf[k]= mulNTL (i.getItem(), remainder);
1271    if (degree (bufFactors[k], x) > 0)
1272      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
1273    else
1274      buf[k]= modNTL (buf[k], bufFactors[k]);
1275  }
1276
1277  for (k= 0; k < factors.length(); k++)
1278    bufFactors[k] += xToJ*buf[k];
1279
1280  // update Pi [0]
1281  int degBuf0= degree (bufFactors[0], x);
1282  int degBuf1= degree (bufFactors[1], x);
1283  if (degBuf0 > 0 && degBuf1 > 0)
1284  {
1285    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
1286    if (j + 2 <= M.rows())
1287      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
1288  }
1289  CanonicalForm uIZeroJ;
1290  if (degBuf0 > 0 && degBuf1 > 0)
1291    uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]);
1292  else if (degBuf0 > 0)
1293    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
1294  else if (degBuf1 > 0)
1295    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
1296  else
1297    uIZeroJ= 0;
1298  Pi [0] += xToJ*uIZeroJ;
1299
1300  CFArray tmp= CFArray (factors.length() - 1);
1301  for (k= 0; k < factors.length() - 1; k++)
1302    tmp[k]= 0;
1303  CFIterator one, two;
1304  one= bufFactors [0];
1305  two= bufFactors [1];
1306  if (degBuf0 > 0 && degBuf1 > 0)
1307  {
1308    while (one.hasTerms() && one.exp() > j) one++;
1309    while (two.hasTerms() && two.exp() > j) two++;
1310    for (k= 1; k <= (int) ceil (j/2.0); k++)
1311    {
1312      if (one.hasTerms() && two.hasTerms())
1313      {
1314        if (k != j - k + 1)
1315        {
1316          if ((one.hasTerms() && one.exp() == j - k + 1) && +
1317              (two.hasTerms() && two.exp() == j - k + 1))
1318        {
1319          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
1320                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
1321          one++;
1322          two++;
1323        }
1324        else if (one.hasTerms() && one.exp() == j - k + 1)
1325        {
1326          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
1327                      M (k + 1, 1);
1328          one++;
1329        }
1330        else if (two.hasTerms() && two.exp() == j - k + 1)
1331        {
1332          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
1333                    M (k + 1, 1);
1334          two++;
1335        }
1336      }
1337      else
1338        tmp[0] += M (k + 1, 1);
1339      }
1340    }
1341  }
1342
1343  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
1344  {
1345    if (j + 2 <= M.rows())
1346      tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
1347                         (bufFactors [1] [j + 1] + bufFactors [1] [0]))
1348                         - M(1,1) - M (j + 2,1);
1349  }
1350  else if (degBuf0 >= j + 1)
1351  {
1352    if (degBuf1 > 0)
1353      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
1354    else
1355      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]);
1356  }
1357  else if (degBuf1 >= j + 1)
1358  {
1359    if (degBuf0 > 0)
1360      tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
1361    else
1362      tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]);
1363  }
1364
1365  Pi [0] += tmp[0]*xToJ*F.mvar();
1366
1367  int degPi, degBuf;
1368  for (int l= 1; l < factors.length() - 1; l++)
1369  {
1370    degPi= degree (Pi [l - 1], x);
1371    degBuf= degree (bufFactors[l + 1], x);
1372    if (degPi > 0 && degBuf > 0)
1373    {
1374      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
1375      if (j + 2 <= M.rows())
1376        M (j + 2, l + 1)= mulNTL (Pi [l - 1][j + 1], bufFactors[l + 1] [j + 1]);
1377    }
1378
1379    if (degPi > 0 && degBuf > 0)
1380      uIZeroJ= mulNTL (Pi[l -1] [0], buf[l + 1]) +
1381               mulNTL (uIZeroJ, bufFactors[l+1] [0]);
1382    else if (degPi > 0)
1383      uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]);
1384    else if (degBuf > 0)
1385      uIZeroJ= mulNTL (Pi[l - 1], buf[1]);
1386    else
1387      uIZeroJ= 0;
1388
1389    Pi [l] += xToJ*uIZeroJ;
1390
1391    one= bufFactors [l + 1];
1392    two= Pi [l - 1];
1393    if (degBuf > 0 && degPi > 0)
1394    {
1395      while (one.hasTerms() && one.exp() > j) one++;
1396      while (two.hasTerms() && two.exp() > j) two++;
1397      for (k= 1; k <= (int) ceil (j/2.0); k++)
1398      {
1399        if (k != j - k + 1)
1400        {
1401          if ((one.hasTerms() && one.exp() == j - k + 1) &&
1402              (two.hasTerms() && two.exp() == j - k + 1))
1403          {
1404            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
1405                      (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1) -
1406                      M (j - k + 2, l + 1);
1407            one++;
1408            two++;
1409          }
1410          else if (one.hasTerms() && one.exp() == j - k + 1)
1411          {
1412            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
1413                               Pi[l - 1] [k]) - M (k + 1, l + 1);
1414            one++;
1415          }
1416          else if (two.hasTerms() && two.exp() == j - k + 1)
1417          {
1418            tmp[l] += mulNTL (bufFactors[l + 1] [k],
1419                      (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1);
1420            two++;
1421           }
1422        }
1423        else
1424          tmp[l] += M (k + 1, l + 1);
1425      }
1426    }
1427
1428    if (degPi >= j + 1 && degBuf >= j + 1)
1429    {
1430      if (j + 2 <= M.rows())
1431        tmp [l] += mulNTL ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
1432                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
1433                          ) - M(1,l+1) - M (j + 2,l+1);
1434    }
1435    else if (degPi >= j + 1)
1436    {
1437      if (degBuf > 0)
1438        tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1] [0]);
1439      else
1440        tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1]);
1441    }
1442    else if (degBuf >= j + 1)
1443    {
1444      if (degPi > 0)
1445        tmp[l] += mulNTL (Pi [l - 1] [0], bufFactors [l + 1] [j + 1]);
1446      else
1447        tmp[l] += mulNTL (Pi [l - 1], bufFactors [l + 1] [j + 1]);
1448    }
1449
1450    Pi[l] += tmp[l]*xToJ*F.mvar();
1451  }
1452  return;
1453}
1454
1455void
1456nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l,
1457                      CFArray& Pi, CFList& diophant, CFMatrix& M,
1458                      const CFArray& LCs, bool sort)
1459{
1460  if (sort)
1461    sortList (factors, Variable (1));
1462  Pi= CFArray (factors.length() - 2);
1463  CFList bufFactors2= factors;
1464  bufFactors2.removeFirst();
1465  diophant= diophantine (F[0], bufFactors2);
1466  DEBOUTLN (cerr, "diophant= " << diophant);
1467
1468  CFArray bufFactors= CFArray (bufFactors2.length());
1469  int i= 0;
1470  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
1471    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
1472
1473  Variable x= F.mvar();
1474  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
1475  {
1476    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
1477    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
1478                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
1479  }
1480  else if (degree (bufFactors[0], x) > 0)
1481  {
1482    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
1483    Pi [0]= M (1, 1) +
1484            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
1485  }
1486  else if (degree (bufFactors[1], x) > 0)
1487  {
1488    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
1489    Pi [0]= M (1, 1) +
1490            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
1491  }
1492  else
1493  {
1494    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
1495    Pi [0]= M (1, 1);
1496  }
1497
1498  for (i= 1; i < Pi.size(); i++)
1499  {
1500    if (degree (Pi[i-1], x) > 0 && degree (bufFactors [i+1], x) > 0)
1501    {
1502      M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors[i+1] [0]);
1503      Pi [i]= M (1,i+1) + (mulNTL (Pi[i-1] [1], bufFactors[i+1] [0]) +
1504                       mulNTL (Pi[i-1] [0], bufFactors [i+1] [1]))*x;
1505    }
1506    else if (degree (Pi[i-1], x) > 0)
1507    {
1508      M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors [i+1]);
1509      Pi [i]=  M(1,i+1) + mulNTL (Pi[i-1] [1], bufFactors[i+1])*x;
1510    }
1511    else if (degree (bufFactors[i+1], x) > 0)
1512    {
1513      M (1,i+1)= mulNTL (Pi[i-1], bufFactors [i+1] [0]);
1514      Pi [i]= M (1,i+1) + mulNTL (Pi[i-1], bufFactors[i+1] [1])*x;
1515    }
1516    else
1517    {
1518      M (1,i+1)= mulNTL (Pi [i-1], bufFactors [i+1]);
1519      Pi [i]= M (1,i+1);
1520    }
1521  }
1522
1523  for (i= 1; i < l; i++)
1524    nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
1525
1526  factors= CFList();
1527  for (i= 0; i < bufFactors.size(); i++)
1528    factors.append (bufFactors[i]);
1529  return;
1530}
1531
1532
1533/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
1534/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
1535CFList
1536diophantine (const CFList& recResult, const CFList& factors,
1537             const CFList& products, const CFList& M, const CanonicalForm& E,
1538             bool& bad)
1539{
1540  if (M.isEmpty())
1541  {
1542    CFList result;
1543    CFListIterator j= factors;
1544    CanonicalForm buf;
1545    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1546    {
1547      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
1548              "constant or univariate poly expected");
1549      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
1550              "constant or univariate poly expected");
1551      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
1552              "constant or univariate poly expected");
1553      buf= mulNTL (E, i.getItem());
1554      result.append (modNTL (buf, j.getItem()));
1555    }
1556    return result;
1557  }
1558  Variable y= M.getLast().mvar();
1559  CFList bufFactors= factors;
1560  for (CFListIterator i= bufFactors; i.hasItem(); i++)
1561    i.getItem()= mod (i.getItem(), y);
1562  CFList bufProducts= products;
1563  for (CFListIterator i= bufProducts; i.hasItem(); i++)
1564    i.getItem()= mod (i.getItem(), y);
1565  CFList buf= M;
1566  buf.removeLast();
1567  CanonicalForm bufE= mod (E, y);
1568  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
1569                                      bufE, bad);
1570
1571  if (bad)
1572    return CFList();
1573
1574  CanonicalForm e= E;
1575  CFListIterator j= products;
1576  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
1577    e -= i.getItem()*j.getItem();
1578
1579  CFList result= recDiophantine;
1580  int d= degree (M.getLast());
1581  CanonicalForm coeffE;
1582  for (int i= 1; i < d; i++)
1583  {
1584    if (degree (e, y) > 0)
1585      coeffE= e[i];
1586    else
1587      coeffE= 0;
1588    if (!coeffE.isZero())
1589    {
1590      CFListIterator k= result;
1591      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
1592                                   coeffE, bad);
1593      if (bad)
1594        return CFList();
1595      CFListIterator l= products;
1596      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
1597      {
1598        k.getItem() += j.getItem()*power (y, i);
1599        e -= j.getItem()*power (y, i)*l.getItem();
1600      }
1601    }
1602    if (e.isZero())
1603      break;
1604  }
1605  if (!e.isZero())
1606  {
1607    bad= true;
1608    return CFList();
1609  }
1610  return result;
1611}
1612
1613void
1614nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
1615                    CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1616                    CFArray& Pi, const CFList& products, int j,
1617                    const CFList& MOD, bool& noOneToOne)
1618{
1619  CanonicalForm E;
1620  CanonicalForm xToJ= power (F.mvar(), j);
1621  Variable x= F.mvar();
1622
1623  // compute the error
1624#ifdef DEBUGOUTPUT
1625    CanonicalForm test= 1;
1626    for (int i= 0; i < factors.length(); i++)
1627    {
1628      if (i == 0)
1629        test *= mod (bufFactors [i], power (x, j));
1630      else
1631        test *= bufFactors[i];
1632    }
1633    test= mod (test, power (x, j));
1634    test= mod (test, MOD);
1635    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
1636    DEBOUTLN (cerr, "test= " << test2);
1637#endif
1638
1639  if (degree (Pi [factors.length() - 2], x) > 0)
1640    E= F[j] - Pi [factors.length() - 2] [j];
1641  else
1642    E= F[j];
1643
1644  CFArray buf= CFArray (diophant.length());
1645
1646  // actual lifting
1647  CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
1648                                    noOneToOne);
1649
1650  if (noOneToOne)
1651    return;
1652
1653  int k= 0;
1654  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
1655  {
1656    buf[k]= i.getItem();
1657    bufFactors[k] += xToJ*i.getItem();
1658  }
1659
1660  // update Pi [0]
1661  int degBuf0= degree (bufFactors[0], x);
1662  int degBuf1= degree (bufFactors[1], x);
1663  if (degBuf0 > 0 && degBuf1 > 0)
1664  {
1665    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
1666    if (j + 2 <= M.rows())
1667      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
1668  }
1669  CanonicalForm uIZeroJ;
1670
1671  if (degBuf0 > 0 && degBuf1 > 0)
1672    uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
1673             mulMod (bufFactors[1] [0], buf[0], MOD);
1674  else if (degBuf0 > 0)
1675    uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
1676  else if (degBuf1 > 0)
1677    uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
1678  else
1679    uIZeroJ= 0;
1680  Pi [0] += xToJ*uIZeroJ;
1681
1682  CFArray tmp= CFArray (factors.length() - 1);
1683  for (k= 0; k < factors.length() - 1; k++)
1684    tmp[k]= 0;
1685  CFIterator one, two;
1686  one= bufFactors [0];
1687  two= bufFactors [1];
1688  if (degBuf0 > 0 && degBuf1 > 0)
1689  {
1690    while (one.hasTerms() && one.exp() > j) one++;
1691    while (two.hasTerms() && two.exp() > j) two++;
1692    for (k= 1; k <= (int) ceil (j/2.0); k++)
1693    {
1694      if (k != j - k + 1)
1695      {
1696        if ((one.hasTerms() && one.exp() == j - k + 1) &&
1697            (two.hasTerms() && two.exp() == j - k + 1))
1698        {
1699          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
1700                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
1701                    M (j - k + 2, 1);
1702          one++;
1703          two++;
1704        }
1705        else if (one.hasTerms() && one.exp() == j - k + 1)
1706        {
1707          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
1708                    bufFactors[1] [k], MOD) - M (k + 1, 1);
1709          one++;
1710        }
1711        else if (two.hasTerms() && two.exp() == j - k + 1)
1712        {
1713          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
1714                    two.coeff()), MOD) - M (k + 1, 1);
1715          two++;
1716        }
1717      }
1718      else
1719      {
1720        tmp[0] += M (k + 1, 1);
1721      }
1722    }
1723  }
1724
1725  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
1726  {
1727    if (j + 2 <= M.rows())
1728      tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
1729                         (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
1730                         - M(1,1) - M (j + 2,1);
1731  }
1732  else if (degBuf0 >= j + 1)
1733  {
1734    if (degBuf1 > 0)
1735      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
1736    else
1737      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
1738  }
1739  else if (degBuf1 >= j + 1)
1740  {
1741    if (degBuf0 > 0)
1742      tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
1743    else
1744      tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
1745  }
1746  Pi [0] += tmp[0]*xToJ*F.mvar();
1747
1748  // update Pi [l]
1749  int degPi, degBuf;
1750  for (int l= 1; l < factors.length() - 1; l++)
1751  {
1752    degPi= degree (Pi [l - 1], x);
1753    degBuf= degree (bufFactors[l + 1], x);
1754    if (degPi > 0 && degBuf > 0)
1755    {
1756      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
1757      if (j + 2 <= M.rows())
1758        M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
1759                                  MOD);
1760    }
1761
1762    if (degPi > 0 && degBuf > 0)
1763      uIZeroJ= mulMod (Pi[l -1] [0], buf[l + 1], MOD) +
1764               mulMod (uIZeroJ, bufFactors[l+1] [0], MOD);
1765    else if (degPi > 0)
1766      uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD);
1767    else if (degBuf > 0)
1768      uIZeroJ= mulMod (Pi[l - 1], buf[1], MOD);
1769    else
1770      uIZeroJ= 0;
1771
1772    Pi [l] += xToJ*uIZeroJ;
1773
1774    one= bufFactors [l + 1];
1775    two= Pi [l - 1];
1776    if (degBuf > 0 && degPi > 0)
1777    {
1778      while (one.hasTerms() && one.exp() > j) one++;
1779      while (two.hasTerms() && two.exp() > j) two++;
1780      for (k= 1; k <= (int) ceil (j/2.0); k++)
1781      {
1782        if (k != j - k + 1)
1783        {
1784          if ((one.hasTerms() && one.exp() == j - k + 1) &&
1785              (two.hasTerms() && two.exp() == j - k + 1))
1786          {
1787            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1788                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
1789                      M (j - k + 2, l + 1);
1790            one++;
1791            two++;
1792          }
1793          else if (one.hasTerms() && one.exp() == j - k + 1)
1794          {
1795            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1796                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
1797            one++;
1798          }
1799          else if (two.hasTerms() && two.exp() == j - k + 1)
1800          {
1801            tmp[l] += mulMod (bufFactors[l + 1] [k],
1802                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
1803            two++;
1804           }
1805        }
1806        else
1807          tmp[l] += M (k + 1, l + 1);
1808      }
1809    }
1810
1811    if (degPi >= j + 1 && degBuf >= j + 1)
1812    {
1813      if (j + 2 <= M.rows())
1814        tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
1815                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
1816                            , MOD) - M(1,l+1) - M (j + 2,l+1);
1817    }
1818    else if (degPi >= j + 1)
1819    {
1820      if (degBuf > 0)
1821        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
1822      else
1823        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
1824    }
1825    else if (degBuf >= j + 1)
1826    {
1827      if (degPi > 0)
1828        tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
1829      else
1830        tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
1831    }
1832
1833    Pi[l] += tmp[l]*xToJ*F.mvar();
1834  }
1835  return;
1836}
1837
1838// wrt. Variable (1)
1839CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
1840{
1841  if (degree (F, 1) <= 0)
1842   return c;
1843  else
1844  {
1845    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
1846    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
1847              - LC (result))*power (result.mvar(), degree (result));
1848    return swapvar (result, Variable (F.level() + 1), Variable (1));
1849  }
1850}
1851
1852CFList
1853nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
1854                      diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
1855                      const CFList& LCs2, bool& bad)
1856{
1857  CFList buf= factors;
1858  int k= 0;
1859  int liftBoundBivar= l[k];
1860  CFList bufbuf= factors;
1861  Variable v= Variable (2);
1862
1863  CFList MOD;
1864  MOD.append (power (Variable (2), liftBoundBivar));
1865  CFArray bufFactors= CFArray (factors.length());
1866  k= 0;
1867  CFListIterator j= eval;
1868  j++;
1869  CFListIterator iter1= LCs1;
1870  CFListIterator iter2= LCs2;
1871  iter1++;
1872  iter2++;
1873  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
1874  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
1875
1876  CFListIterator i= buf;
1877  i++;
1878  Variable y= j.getItem().mvar();
1879  if (y.level() != 3)
1880    y= Variable (3);
1881
1882  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
1883  M (1, 1)= Pi[0];
1884  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
1885    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
1886                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
1887  else if (degree (bufFactors[0], y) > 0)
1888    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
1889  else if (degree (bufFactors[1], y) > 0)
1890    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
1891
1892  CFList products;
1893  for (int i= 0; i < bufFactors.size(); i++)
1894  {
1895    if (degree (bufFactors[i], y) > 0)
1896      products.append (eval.getFirst()/bufFactors[i] [0]);
1897    else
1898      products.append (eval.getFirst()/bufFactors[i]);
1899  }
1900
1901  for (int d= 1; d < l[1]; d++)
1902  {
1903    nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
1904                        d, MOD, bad);
1905    if (bad)
1906      return CFList();
1907  }
1908  CFList result;
1909  for (k= 0; k < factors.length(); k++)
1910    result.append (bufFactors[k]);
1911  return result;
1912}
1913
1914
1915CFList
1916nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
1917                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
1918                    int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
1919                   )
1920{
1921  int k= 0;
1922  CFArray bufFactors= CFArray (factors.length());
1923  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
1924  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
1925  CFList buf= factors;
1926  Variable y= F.getLast().mvar();
1927  Variable x= F.getFirst().mvar();
1928  CanonicalForm xToLOld= power (x, lOld);
1929  Pi [0]= mod (Pi[0], xToLOld);
1930  M (1, 1)= Pi [0];
1931
1932  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
1933    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
1934                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
1935  else if (degree (bufFactors[0], y) > 0)
1936    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
1937  else if (degree (bufFactors[1], y) > 0)
1938    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
1939
1940  CFList products;
1941  CanonicalForm quot;
1942  for (int i= 0; i < bufFactors.size(); i++)
1943  {
1944    if (degree (bufFactors[i], y) > 0)
1945    {
1946      if (!fdivides (bufFactors[i] [0], F.getFirst(), quot))
1947      {
1948        bad= true;
1949        return CFList();
1950      }
1951      products.append (quot);
1952    }
1953    else
1954    {
1955      if (!fdivides (bufFactors[i], F.getFirst(), quot))
1956      {
1957        bad= true;
1958        return CFList();
1959      }
1960      products.append (quot);
1961    }
1962  }
1963
1964  for (int d= 1; d < lNew; d++)
1965  {
1966    nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
1967                        d, MOD, bad);
1968    if (bad)
1969      return CFList();
1970  }
1971
1972  CFList result;
1973  for (k= 0; k < factors.length(); k++)
1974    result.append (bufFactors[k]);
1975  return result;
1976}
1977
1978CFList
1979nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
1980                    lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
1981                    const CFArray& Pi, const CFList& diophant, bool& bad)
1982{
1983  CFList bufDiophant= diophant;
1984  CFList buf= factors;
1985  if (sort)
1986    sortList (buf, Variable (1));
1987  CFArray bufPi= Pi;
1988  CFMatrix M= CFMatrix (l[1], factors.length());
1989  CFList result= 
1990    nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
1991  if (bad)
1992    return CFList();
1993
1994  if (eval.length() == 2)
1995    return result;
1996  CFList MOD;
1997  for (int i= 0; i < 2; i++)
1998    MOD.append (power (Variable (i + 2), l[i]));
1999  CFListIterator j= eval;
2000  j++;
2001  CFList bufEval;
2002  bufEval.append (j.getItem());
2003  j++;
2004  CFListIterator jj= LCs1;
2005  CFListIterator jjj= LCs2;
2006  CFList bufLCs1, bufLCs2;
2007  jj++, jjj++;
2008  bufLCs1.append (jj.getItem());
2009  bufLCs2.append (jjj.getItem());
2010  jj++, jjj++;
2011
2012  for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
2013  {
2014    bufEval.append (j.getItem());
2015    bufLCs1.append (jj.getItem());
2016    bufLCs2.append (jjj.getItem());
2017    M= CFMatrix (l[i], factors.length());
2018    result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
2019                                 l[i - 1], l[i], bufLCs1, bufLCs2, bad);
2020    if (bad)
2021      return CFList();
2022    MOD.append (power (Variable (i + 2), l[i]));
2023    bufEval.removeFirst();
2024    bufLCs1.removeFirst();
2025    bufLCs2.removeFirst();
2026  }
2027  return result;
2028}
2029
2030CFList
2031nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
2032                      CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
2033                      int bivarLiftBound, bool& bad)
2034{
2035  CFList bufFactors2= factors;
2036
2037  Variable y= Variable (2);
2038  for (CFListIterator i= bufFactors2; i.hasItem(); i++)
2039    i.getItem()= mod (i.getItem(), y);
2040
2041  CanonicalForm bufF= F;
2042  bufF= mod (bufF, y);
2043  bufF= mod (bufF, Variable (3));
2044
2045  diophant= diophantine (bufF, bufFactors2);
2046
2047  CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
2048
2049  Pi= CFArray (bufFactors2.length() - 1);
2050
2051  CFArray bufFactors= CFArray (bufFactors2.length());
2052  CFListIterator j= LCs;
2053  int i= 0;
2054  for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
2055    bufFactors[i]= replaceLC (k.getItem(), j.getItem());
2056
2057  //initialise Pi
2058  Variable v= Variable (3);
2059  CanonicalForm yToL= power (y, bivarLiftBound);
2060  if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
2061  {
2062    M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
2063    Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
2064                       mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
2065  }
2066  else if (degree (bufFactors[0], v) > 0)
2067  {
2068    M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
2069    Pi [0]=  M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
2070  }
2071  else if (degree (bufFactors[1], v) > 0)
2072  {
2073    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
2074    Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
2075  }
2076  else
2077  {
2078    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
2079    Pi [0]= M (1,1);
2080  }
2081
2082  for (i= 1; i < Pi.size(); i++)
2083  {
2084    if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
2085    {
2086      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
2087      Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
2088                       mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
2089    }
2090    else if (degree (Pi[i-1], v) > 0)
2091    {
2092      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
2093      Pi [i]=  M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
2094    }
2095    else if (degree (bufFactors[i+1], v) > 0)
2096    {
2097      M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
2098      Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
2099    }
2100    else
2101    {
2102      M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
2103      Pi [i]= M (1,i+1);
2104    }
2105  }
2106
2107  CFList products;
2108  bufF= mod (F, Variable (3));
2109  for (CFListIterator k= factors; k.hasItem(); k++)
2110    products.append (bufF/k.getItem());
2111
2112  CFList MOD= CFList (power (v, liftBound));
2113  MOD.insert (yToL);
2114  for (int d= 1; d < liftBound; d++)
2115  {
2116    nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
2117                        MOD, bad);
2118    if (bad)
2119      return CFList();
2120  }
2121
2122  CFList result;
2123  for (i= 0; i < factors.length(); i++)
2124    result.append (bufFactors[i]);
2125  return result;
2126}
2127
2128CFList
2129nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
2130                    CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2131                    int& lNew, const CFList& MOD, bool& noOneToOne
2132                   )
2133{
2134
2135  int k= 0;
2136  CFArray bufFactors= CFArray (factors.length());
2137  CFListIterator j= LCs;
2138  for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
2139    bufFactors [k]= replaceLC (i.getItem(), j.getItem());
2140
2141  Variable y= F.getLast().mvar();
2142  Variable x= F.getFirst().mvar();
2143  CanonicalForm xToLOld= power (x, lOld);
2144
2145  Pi [0]= mod (Pi[0], xToLOld);
2146  M (1, 1)= Pi [0];
2147
2148  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2149    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2150                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2151  else if (degree (bufFactors[0], y) > 0)
2152    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2153  else if (degree (bufFactors[1], y) > 0)
2154    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2155
2156  for (int i= 1; i < Pi.size(); i++)
2157  {
2158    Pi [i]= mod (Pi [i], xToLOld);
2159    M (1, i + 1)= Pi [i];
2160
2161    if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
2162      Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
2163                 mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
2164    else if (degree (Pi[i-1], y) > 0)
2165      Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
2166    else if (degree (bufFactors[i+1], y) > 0)
2167      Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
2168  }
2169
2170  CFList products;
2171  CanonicalForm quot, bufF= F.getFirst();
2172
2173  for (int i= 0; i < bufFactors.size(); i++)
2174  {
2175    if (degree (bufFactors[i], y) > 0)
2176    {
2177      if (!fdivides (bufFactors[i] [0], bufF, quot))
2178      {
2179        noOneToOne= true;
2180        return factors;
2181      }
2182      products.append (quot);
2183    }
2184    else
2185    {
2186      if (!fdivides (bufFactors[i], bufF, quot))
2187      {
2188        noOneToOne= true;
2189        return factors;
2190      }
2191      products.append (quot);
2192    }
2193  }
2194
2195  for (int d= 1; d < lNew; d++)
2196  {
2197    nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
2198                        products, d, MOD, noOneToOne);
2199    if (noOneToOne)
2200      return CFList();
2201  }
2202
2203  CFList result;
2204  for (k= 0; k < factors.length(); k++)
2205    result.append (bufFactors[k]);
2206  return result;
2207}
2208
2209CFList
2210nonMonicHenselLift (const CFList& eval, const CFList& factors,
2211                    CFList* const& LCs, CFList& diophant, CFArray& Pi,
2212                    int* liftBound, int length, bool& noOneToOne
2213                   )
2214{
2215  CFList bufDiophant= diophant;
2216  CFList buf= factors;
2217  CFArray bufPi= Pi;
2218  CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
2219  int k= 0;
2220
2221  CFList result=
2222  nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
2223                        liftBound[1], liftBound[0], noOneToOne);
2224
2225  if (noOneToOne)
2226    return CFList();
2227
2228  if (eval.length() == 1)
2229    return result;
2230
2231  k++;
2232  CFList MOD;
2233  for (int i= 0; i < 2; i++)
2234    MOD.append (power (Variable (i + 2), liftBound[i]));
2235
2236  CFListIterator j= eval;
2237  CFList bufEval;
2238  bufEval.append (j.getItem());
2239  j++;
2240
2241  for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
2242  {
2243    bufEval.append (j.getItem());
2244    M= CFMatrix (liftBound[i], factors.length() - 1);
2245    result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
2246                                liftBound[i-1], liftBound[i], MOD, noOneToOne);
2247    if (noOneToOne)
2248      return result;
2249    MOD.append (power (Variable (i + 2), liftBound[i]));
2250    bufEval.removeFirst();
2251  }
2252
2253  return result;
2254}
2255
2256#endif
2257/* HAVE_NTL */
2258
Note: See TracBrowser for help on using the repository browser.