source: git/factory/facHensel.cc @ 4e6d2a

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