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

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