source: git/factory/facHensel.cc @ a090c88

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