source: git/factory/facHensel.cc @ 1ddcde9

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