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

spielwiese
Last change on this file since 4a05ed was 4a05ed, checked in by Martin Lee <martinlee84@…>, 12 years ago
added modular solving of diophantine equation for Hensel lifting over Q(a) reduce result of tryExtgcd modulo minimal polynomial git-svn-id: file:///usr/local/Singular/svn/trunk@14354 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 89.6 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file facHensel.cc
5 *
6 * This file implements functions to lift factors via Hensel lifting and
7 * functions for modular multiplication and division with remainder.
8 *
9 * ABSTRACT: Hensel lifting is described in "Efficient Multivariate
10 * Factorization over Finite Fields" by L. Bernardin & M. Monagon. Division with
11 * remainder is described in "Fast Recursive Division" by C. Burnikel and
12 * J. Ziegler. Karatsuba multiplication is described in "Modern Computer
13 * Algebra" by J. von zur Gathen and J. Gerhard.
14 *
15 * @author Martin Lee
16 *
17 * @internal @version \$Id$
18 *
19 **/
20/*****************************************************************************/
21
22#include "assert.h"
23#include "debug.h"
24#include "timing.h"
25
26#include "facHensel.h"
27#include "cf_util.h"
28#include "fac_util.h"
29#include "cf_algorithm.h"
30#include "cf_primes.h"
31
32#ifdef HAVE_NTL
33#include <NTL/lzz_pEX.h>
34#include "NTLconvert.h"
35
36static
37CFList productsNTL (const CFList& factors, const CanonicalForm& M)
38{
39  zz_p::init (getCharacteristic());
40  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
41  zz_pE::init (NTLMipo);
42  zz_pEX prod;
43  vec_zz_pEX v;
44  v.SetLength (factors.length());
45  int j= 0;
46  for (CFListIterator i= factors; i.hasItem(); i++, j++)
47  {
48    if (i.getItem().inCoeffDomain())
49      v[j]= to_zz_pEX (to_zz_pE (convertFacCF2NTLzzpX (i.getItem())));
50    else
51      v[j]= convertFacCF2NTLzz_pEX (i.getItem(), NTLMipo);
52  }
53  CFList result;
54  Variable x= Variable (1);
55  for (int j= 0; j < factors.length(); j++)
56  {
57    int k= 0;
58    set(prod);
59    for (int i= 0; i < factors.length(); i++, k++)
60    {
61      if (k == j)
62        continue;
63      prod *= v[i];
64    }
65    result.append (convertNTLzz_pEX2CF (prod, x, M.mvar()));
66  }
67  return result;
68}
69
70static
71void tryDiophantine (CFList& result, const CanonicalForm& F,
72                     const CFList& factors, const CanonicalForm& M, bool& fail)
73{
74  ASSERT (M.isUnivariate(), "expected univariate poly");
75
76  CFList bufFactors= factors;
77  bufFactors.removeFirst();
78  bufFactors.insert (factors.getFirst () (0,2));
79  CanonicalForm inv, leadingCoeff= Lc (F);
80  CFListIterator i= bufFactors;
81  if (bufFactors.getFirst().inCoeffDomain())
82  {
83    if (i.hasItem())
84      i++;
85  }
86  for (; i.hasItem(); i++)
87  {
88    tryInvert (Lc (i.getItem()), M, inv ,fail);
89    if (fail)
90      return;
91    i.getItem()= reduce (i.getItem()*inv, M);
92  }
93  bufFactors= productsNTL (bufFactors, M);
94
95  CanonicalForm buf1, buf2, buf3, S, T;
96  i= bufFactors;
97  if (i.hasItem())
98    i++;
99  buf1= bufFactors.getFirst();
100  buf2= i.getItem();
101  tryExtgcd (buf1, buf2, M, buf3, S, T, fail);
102  if (fail)
103    return;
104  result.append (S);
105  result.append (T);
106  if (i.hasItem())
107    i++;
108  for (; i.hasItem(); i++)
109  {
110    buf1= i.getItem();
111    CanonicalForm bb= buf3;
112    tryExtgcd (buf3, buf1, M, buf3, S, T, fail);
113    if (fail)
114      return;
115    CFListIterator k= factors;
116    for (CFListIterator j= result; j.hasItem(); j++, k++)
117    {
118      j.getItem() *= S;
119      CanonicalForm q;
120      j.getItem()= mod (j.getItem(), k.getItem());
121      j.getItem()= reduce (j.getItem(), M);
122    }
123    result.append (T);
124  }
125}
126
127static inline
128CFList mapinto (const CFList& L)
129{
130  CFList result;
131  for (CFListIterator i= L; i.hasItem(); i++)
132    result.append (mapinto (i.getItem()));
133  return result;
134}
135
136static inline
137int mod (const CFList& L, const CanonicalForm& p)
138{
139  for (CFListIterator i= L; i.hasItem(); i++)
140  {
141    if (mod (i.getItem(), p) == 0)
142      return 0;
143  }
144  return 1;
145}
146
147
148static inline void
149chineseRemainder (const CFList & x1, const CanonicalForm & q1,
150                  const CFList & x2, const CanonicalForm & q2,
151                  CFList & xnew, CanonicalForm & qnew)
152{
153  ASSERT (x1.length() == x2.length(), "expected lists of equal length");
154  CanonicalForm tmp1, tmp2;
155  CFListIterator j= x2;
156  for (CFListIterator i= x1; i.hasItem() && j.hasItem(); i++, j++)
157  {
158    chineseRemainder (i.getItem(), q1, j.getItem(), q2, tmp1, tmp2);
159    xnew.append (tmp1);
160  }
161  qnew= tmp2;
162}
163
164static inline
165CFList Farey (const CFList& L, const CanonicalForm& q)
166{
167  CFList result;
168  for (CFListIterator i= L; i.hasItem(); i++)
169    result.append (Farey (i.getItem(), q));
170  return result;
171}
172
173static inline
174CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
175{
176  CFList result;
177  for (CFListIterator i= L; i.hasItem(); i++)
178    result.append (replacevar (i.getItem(), a, b));
179  return result;
180}
181
182CFList
183modularDiophant (const CanonicalForm& f, const CFList& factors,
184                 const CanonicalForm& M)
185{
186  bool save_rat=!isOn (SW_RATIONAL);
187  On (SW_RATIONAL);
188  CanonicalForm F= f*bCommonDen (f);
189  CFList products= factors;
190  for (CFListIterator i= products; i.hasItem(); i++)
191    i.getItem() *= bCommonDen (i.getItem());
192  if (products.getFirst().level() == 1)
193    products.insert (Lc (F));
194  CanonicalForm bound= maxNorm (F);
195  CFList leadingCoeffs;
196  leadingCoeffs.append (lc (F));
197  CanonicalForm dummy;
198  for (CFListIterator i= products; i.hasItem(); i++)
199  {
200    leadingCoeffs.append (lc (i.getItem()));
201    dummy= maxNorm (i.getItem());
202    bound= (dummy > bound) ? dummy : bound;
203  }
204  bound *= maxNorm (Lc (F))*maxNorm (Lc(F))*bound;
205  bound *= bound*bound;
206  bound= power (bound, degree (M));
207  bound *= power (CanonicalForm (2),degree (f));
208  CanonicalForm bufBound= bound;
209  int i = cf_getNumBigPrimes() - 1;
210  int p;
211  CFList resultModP, result, newResult;
212  CanonicalForm q (0), newQ;
213  bool fail= false;
214  Variable a= M.mvar();
215  Variable b= Variable (2);
216  setReduce (M.mvar(), false);
217  CanonicalForm mipo= bCommonDen (M)*M;
218  Off (SW_RATIONAL);
219  CanonicalForm modMipo;
220  leadingCoeffs.append (lc (mipo));
221  CFList tmp1, tmp2;
222  bool equal= false;
223  int count= 0;
224  do
225  {
226    p = cf_getBigPrime( i );
227    i--;
228    while ( i >= 0 && mod( leadingCoeffs, p ) == 0)
229    {
230      p = cf_getBigPrime( i );
231      i--;
232    }
233
234    ASSERT (i >= 0, "ran out of primes"); //sic
235
236    setCharacteristic (p);
237    modMipo= mapinto (mipo);
238    modMipo /= lc (modMipo);
239    resultModP= CFList();
240    tryDiophantine (resultModP, mapinto (F), mapinto (products), modMipo, fail);
241    setCharacteristic (0);
242    if (fail)
243    {
244      fail= false;
245      continue;
246    }
247
248    if ( q.isZero() )
249    {
250      result= replacevar (mapinto(resultModP), a, b);
251      q= p;
252    }
253    else
254    {
255      result= replacevar (result, a, b);
256      newResult= CFList();
257      chineseRemainder( result, q, replacevar (mapinto (resultModP), a, b),
258                        p, newResult, newQ );
259      q= newQ;
260      result= newResult;
261      if (newQ > bound)
262      {
263        count++;
264        tmp1= replacevar (Farey (result, q), b, a);
265        if (tmp2.isEmpty())
266          tmp2= tmp1;
267        else
268        {
269          equal= true;
270          CFListIterator k= tmp1;
271          for (CFListIterator j= tmp2; j.hasItem(); j++, k++)
272          {
273            if (j.getItem() != k.getItem())
274              equal= false;
275          }
276          if (!equal)
277            tmp2= tmp1;
278        }
279        if (count > 2)
280        {
281          bound *= bufBound;
282          equal= false;
283          count= 0;
284        }
285      }
286      if (newQ > bound && equal)
287      {
288        On( SW_RATIONAL );
289        result= tmp2;
290        CFList bufResult= result;
291        setReduce (M.mvar(), true);
292        if (factors.getFirst().level() == 1)
293        {
294          result.removeFirst();
295          CFListIterator j= factors;
296          CanonicalForm denf= bCommonDen (f);
297          for (CFListIterator i= result; i.hasItem(); i++, j++)
298            i.getItem() *= Lc (j.getItem())*denf;
299        }
300        if (factors.getFirst().level() != 1 && 
301            !bCommonDen (factors.getFirst()).isOne())
302        {
303          CanonicalForm denFirst= bCommonDen (factors.getFirst());
304          for (CFListIterator i= result; i.hasItem(); i++)
305            i.getItem() *= denFirst;
306        }
307
308        CanonicalForm test= 0;
309        CFListIterator jj= factors;
310        for (CFListIterator ii= result; ii.hasItem(); ii++, jj++)
311          test += ii.getItem()*(f/jj.getItem());
312        if (!test.isOne())
313        {
314          bound *= bufBound;
315          equal= false;
316          count= 0;
317          setReduce (M.mvar(), false);
318          result= bufResult;
319          Off (SW_RATIONAL);
320        }
321        else
322          break;
323      }
324    }
325  } while (1);
326  if (save_rat) Off(SW_RATIONAL);
327  return result;
328}
329
330CanonicalForm
331mulNTL (const CanonicalForm& F, const CanonicalForm& G)
332{
333  if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
334    return F*G;
335  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
336  ASSERT (F.level() == G.level(), "expected polys of same level");
337  if (CFFactory::gettype() == GaloisFieldDomain)
338    return F*G;
339  zz_p::init (getCharacteristic());
340  Variable alpha;
341  CanonicalForm result;
342  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
343  {
344    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
345    zz_pE::init (NTLMipo);
346    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
347    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
348    mul (NTLF, NTLF, NTLG);
349    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
350  }
351  else
352  {
353    zz_pX NTLF= convertFacCF2NTLzzpX (F);
354    zz_pX NTLG= convertFacCF2NTLzzpX (G);
355    mul (NTLF, NTLF, NTLG);
356    result= convertNTLzzpX2CF(NTLF, F.mvar());
357  }
358  return result;
359}
360
361CanonicalForm
362modNTL (const CanonicalForm& F, const CanonicalForm& G)
363{
364  if (F.inCoeffDomain() && G.isUnivariate())
365    return F;
366  else if (F.inCoeffDomain() && G.inCoeffDomain())
367    return mod (F, G);
368  else if (F.isUnivariate() && G.inCoeffDomain())
369    return mod (F,G);
370
371  if (getCharacteristic() == 0)
372    return mod (F, G);
373
374  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
375  ASSERT (F.level() == G.level(), "expected polys of same level");
376  if (CFFactory::gettype() == GaloisFieldDomain)
377    return mod (F, G);
378  zz_p::init (getCharacteristic());
379  Variable alpha;
380  CanonicalForm result;
381  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
382  {
383    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
384    zz_pE::init (NTLMipo);
385    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
386    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
387    rem (NTLF, NTLF, NTLG);
388    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
389  }
390  else
391  {
392    zz_pX NTLF= convertFacCF2NTLzzpX (F);
393    zz_pX NTLG= convertFacCF2NTLzzpX (G);
394    rem (NTLF, NTLF, NTLG);
395    result= convertNTLzzpX2CF(NTLF, F.mvar());
396  }
397  return result;
398}
399
400CanonicalForm
401divNTL (const CanonicalForm& F, const CanonicalForm& G)
402{
403  if (F.inCoeffDomain() && G.isUnivariate())
404    return F;
405  else if (F.inCoeffDomain() && G.inCoeffDomain())
406    return div (F, G);
407  else if (F.isUnivariate() && G.inCoeffDomain())
408    return div (F,G);
409
410  if (getCharacteristic() == 0)
411    return div (F, G);
412
413  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
414  ASSERT (F.level() == G.level(), "expected polys of same level");
415  if (CFFactory::gettype() == GaloisFieldDomain)
416    return div (F, G);
417  zz_p::init (getCharacteristic());
418  Variable alpha;
419  CanonicalForm result;
420  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
421  {
422    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
423    zz_pE::init (NTLMipo);
424    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
425    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
426    div (NTLF, NTLF, NTLG);
427    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
428  }
429  else
430  {
431    zz_pX NTLF= convertFacCF2NTLzzpX (F);
432    zz_pX NTLG= convertFacCF2NTLzzpX (G);
433    div (NTLF, NTLF, NTLG);
434    result= convertNTLzzpX2CF(NTLF, F.mvar());
435  }
436  return result;
437}
438
439/*
440void
441divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
442           CanonicalForm& R)
443{
444  if (F.inCoeffDomain() && G.isUnivariate())
445  {
446    R= F;
447    Q= 0;
448  }
449  else if (F.inCoeffDomain() && G.inCoeffDomain())
450  {
451    divrem (F, G, Q, R);
452    return;
453  }
454  else if (F.isUnivariate() && G.inCoeffDomain())
455  {
456    divrem (F, G, Q, R);
457    return;
458  }
459
460  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
461  ASSERT (F.level() == G.level(), "expected polys of same level");
462  if (CFFactory::gettype() == GaloisFieldDomain)
463  {
464    divrem (F, G, Q, R);
465    return;
466  }
467  zz_p::init (getCharacteristic());
468  Variable alpha;
469  CanonicalForm result;
470  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
471  {
472    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
473    zz_pE::init (NTLMipo);
474    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
475    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
476    zz_pEX NTLQ;
477    zz_pEX NTLR;
478    DivRem (NTLQ, NTLR, NTLF, NTLG);
479    Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
480    R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
481    return;
482  }
483  else
484  {
485    zz_pX NTLF= convertFacCF2NTLzzpX (F);
486    zz_pX NTLG= convertFacCF2NTLzzpX (G);
487    zz_pX NTLQ;
488    zz_pX NTLR;
489    DivRem (NTLQ, NTLR, NTLF, NTLG);
490    Q= convertNTLzzpX2CF(NTLQ, F.mvar());
491    R= convertNTLzzpX2CF(NTLR, F.mvar());
492    return;
493  }
494}*/
495
496CanonicalForm mod (const CanonicalForm& F, const CFList& M)
497{
498  CanonicalForm A= F;
499  for (CFListIterator i= M; i.hasItem(); i++)
500    A= mod (A, i.getItem());
501  return A;
502}
503
504zz_pX kronSubFp (const CanonicalForm& A, int d)
505{
506  int degAy= degree (A);
507  zz_pX result;
508  result.rep.SetLength (d*(degAy + 1));
509
510  zz_p *resultp;
511  resultp= result.rep.elts();
512  zz_pX buf;
513  zz_p *bufp;
514
515  for (CFIterator i= A; i.hasTerms(); i++)
516  {
517    if (i.coeff().inCoeffDomain())
518      buf= convertFacCF2NTLzzpX (i.coeff());
519    else
520      buf= convertFacCF2NTLzzpX (i.coeff());
521
522    int k= i.exp()*d;
523    bufp= buf.rep.elts();
524    int bufRepLength= (int) buf.rep.length();
525    for (int j= 0; j < bufRepLength; j++)
526      resultp [j + k]= bufp [j];
527  }
528  result.normalize();
529
530  return result;
531}
532
533zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
534{
535  int degAy= degree (A);
536  zz_pEX result;
537  result.rep.SetLength (d*(degAy + 1));
538
539  Variable v;
540  zz_pE *resultp;
541  resultp= result.rep.elts();
542  zz_pEX buf1;
543  zz_pE *buf1p;
544  zz_pX buf2;
545  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
546
547  for (CFIterator i= A; i.hasTerms(); i++)
548  {
549    if (i.coeff().inCoeffDomain())
550    {
551      buf2= convertFacCF2NTLzzpX (i.coeff());
552      buf1= to_zz_pEX (to_zz_pE (buf2));
553    }
554    else
555      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
556
557    int k= i.exp()*d;
558    buf1p= buf1.rep.elts();
559    int buf1RepLength= (int) buf1.rep.length();
560    for (int j= 0; j < buf1RepLength; j++)
561      resultp [j + k]= buf1p [j];
562  }
563  result.normalize();
564
565  return result;
566}
567
568void
569kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
570                const Variable& alpha)
571{
572  int degAy= degree (A);
573  subA1.rep.SetLength ((long) d*(degAy + 1));
574  subA2.rep.SetLength ((long) d*(degAy + 1));
575
576  Variable v;
577  zz_pE *subA1p;
578  zz_pE *subA2p;
579  subA1p= subA1.rep.elts();
580  subA2p= subA2.rep.elts();
581  zz_pEX buf;
582  zz_pE *bufp;
583  zz_pX buf2;
584  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
585
586  for (CFIterator i= A; i.hasTerms(); i++)
587  {
588    if (i.coeff().inCoeffDomain())
589    {
590      buf2= convertFacCF2NTLzzpX (i.coeff());
591      buf= to_zz_pEX (to_zz_pE (buf2));
592    }
593    else
594      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
595
596    int k= i.exp()*d;
597    int kk= (degAy - i.exp())*d;
598    bufp= buf.rep.elts();
599    int bufRepLength= (int) buf.rep.length();
600    for (int j= 0; j < bufRepLength; j++)
601    {
602      subA1p [j + k]= bufp [j];
603      subA2p [j + kk]= bufp [j];
604    }
605  }
606  subA1.normalize();
607  subA2.normalize();
608}
609
610void
611kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
612{
613  int degAy= degree (A);
614  subA1.rep.SetLength ((long) d*(degAy + 1));
615  subA2.rep.SetLength ((long) d*(degAy + 1));
616
617  Variable v;
618  zz_p *subA1p;
619  zz_p *subA2p;
620  subA1p= subA1.rep.elts();
621  subA2p= subA2.rep.elts();
622  zz_pX buf;
623  zz_p *bufp;
624
625  for (CFIterator i= A; i.hasTerms(); i++)
626  {
627    buf= convertFacCF2NTLzzpX (i.coeff());
628
629    int k= i.exp()*d;
630    int kk= (degAy - i.exp())*d;
631    bufp= buf.rep.elts();
632    int bufRepLength= (int) buf.rep.length();
633    for (int j= 0; j < bufRepLength; j++)
634    {
635      subA1p [j + k]= bufp [j];
636      subA2p [j + kk]= bufp [j];
637    }
638  }
639  subA1.normalize();
640  subA2.normalize();
641}
642
643CanonicalForm
644reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
645              const Variable& alpha)
646{
647  Variable y= Variable (2);
648  Variable x= Variable (1);
649
650  zz_pEX f= F;
651  zz_pEX g= G;
652  int degf= deg(f);
653  int degg= deg(g);
654
655  zz_pEX buf1;
656  zz_pEX buf2;
657  zz_pEX buf3;
658
659  zz_pE *buf1p;
660  zz_pE *buf2p;
661  zz_pE *buf3p;
662  if (f.rep.length() < (long) d*(k+1)) //zero padding
663    f.rep.SetLength ((long)d*(k+1));
664
665  zz_pE *gp= g.rep.elts();
666  zz_pE *fp= f.rep.elts();
667  CanonicalForm result= 0;
668  int i= 0;
669  int lf= 0;
670  int lg= d*k;
671  int degfSubLf= degf;
672  int deggSubLg= degg-lg;
673  int repLengthBuf2;
674  int repLengthBuf1;
675  zz_pE zzpEZero= zz_pE();
676
677  while (degf >= lf || lg >= 0)
678  {
679    if (degfSubLf >= d)
680      repLengthBuf1= d;
681    else if (degfSubLf < 0)
682      repLengthBuf1= 0;
683    else
684      repLengthBuf1= degfSubLf + 1;
685    buf1.rep.SetLength((long) repLengthBuf1);
686
687    buf1p= buf1.rep.elts();
688    for (int ind= 0; ind < repLengthBuf1; ind++)
689      buf1p [ind]= fp [ind + lf];
690    buf1.normalize();
691
692    repLengthBuf1= buf1.rep.length();
693
694    if (deggSubLg >= d - 1)
695      repLengthBuf2= d - 1;
696    else if (deggSubLg < 0)
697      repLengthBuf2= 0;
698    else
699      repLengthBuf2= deggSubLg + 1;
700
701    buf2.rep.SetLength ((long) repLengthBuf2);
702    buf2p= buf2.rep.elts();
703    for (int ind= 0; ind < repLengthBuf2; ind++)
704    {
705      buf2p [ind]= gp [ind + lg];
706    }
707    buf2.normalize();
708
709    repLengthBuf2= buf2.rep.length();
710
711    buf3.rep.SetLength((long) repLengthBuf2 + d);
712    buf3p= buf3.rep.elts();
713    buf2p= buf2.rep.elts();
714    buf1p= buf1.rep.elts();
715    for (int ind= 0; ind < repLengthBuf1; ind++)
716      buf3p [ind]= buf1p [ind];
717    for (int ind= repLengthBuf1; ind < d; ind++)
718      buf3p [ind]= zzpEZero;
719    for (int ind= 0; ind < repLengthBuf2; ind++)
720      buf3p [ind + d]= buf2p [ind];
721    buf3.normalize();
722
723    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
724    i++;
725
726
727    lf= i*d;
728    degfSubLf= degf - lf;
729
730    lg= d*(k-i);
731    deggSubLg= degg - lg;
732
733    buf1p= buf1.rep.elts();
734
735    if (lg >= 0 && deggSubLg > 0)
736    {
737      if (repLengthBuf2 > degfSubLf + 1)
738        degfSubLf= repLengthBuf2 - 1;
739      int tmp= tmin (repLengthBuf1, deggSubLg);
740      for (int ind= 0; ind < tmp; ind++)
741        gp [ind + lg] -= buf1p [ind];
742    }
743
744    if (lg < 0)
745      break;
746
747    buf2p= buf2.rep.elts();
748    if (degfSubLf >= 0)
749    {
750      for (int ind= 0; ind < repLengthBuf2; ind++)
751        fp [ind + lf] -= buf2p [ind];
752    }
753  }
754
755  return result;
756}
757
758CanonicalForm
759reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
760{
761  Variable y= Variable (2);
762  Variable x= Variable (1);
763
764  zz_pX f= F;
765  zz_pX g= G;
766  int degf= deg(f);
767  int degg= deg(g);
768
769  zz_pX buf1;
770  zz_pX buf2;
771  zz_pX buf3;
772
773  zz_p *buf1p;
774  zz_p *buf2p;
775  zz_p *buf3p;
776
777  if (f.rep.length() < (long) d*(k+1)) //zero padding
778    f.rep.SetLength ((long)d*(k+1));
779
780  zz_p *gp= g.rep.elts();
781  zz_p *fp= f.rep.elts();
782  CanonicalForm result= 0;
783  int i= 0;
784  int lf= 0;
785  int lg= d*k;
786  int degfSubLf= degf;
787  int deggSubLg= degg-lg;
788  int repLengthBuf2;
789  int repLengthBuf1;
790  zz_p zzpZero= zz_p();
791  while (degf >= lf || lg >= 0)
792  {
793    if (degfSubLf >= d)
794      repLengthBuf1= d;
795    else if (degfSubLf < 0)
796      repLengthBuf1= 0;
797    else
798      repLengthBuf1= degfSubLf + 1;
799    buf1.rep.SetLength((long) repLengthBuf1);
800
801    buf1p= buf1.rep.elts();
802    for (int ind= 0; ind < repLengthBuf1; ind++)
803      buf1p [ind]= fp [ind + lf];
804    buf1.normalize();
805
806    repLengthBuf1= buf1.rep.length();
807
808    if (deggSubLg >= d - 1)
809      repLengthBuf2= d - 1;
810    else if (deggSubLg < 0)
811      repLengthBuf2= 0;
812    else
813      repLengthBuf2= deggSubLg + 1;
814
815    buf2.rep.SetLength ((long) repLengthBuf2);
816    buf2p= buf2.rep.elts();
817    for (int ind= 0; ind < repLengthBuf2; ind++)
818      buf2p [ind]= gp [ind + lg];
819
820    buf2.normalize();
821
822    repLengthBuf2= buf2.rep.length();
823
824
825    buf3.rep.SetLength((long) repLengthBuf2 + d);
826    buf3p= buf3.rep.elts();
827    buf2p= buf2.rep.elts();
828    buf1p= buf1.rep.elts();
829    for (int ind= 0; ind < repLengthBuf1; ind++)
830      buf3p [ind]= buf1p [ind];
831    for (int ind= repLengthBuf1; ind < d; ind++)
832      buf3p [ind]= zzpZero;
833    for (int ind= 0; ind < repLengthBuf2; ind++)
834      buf3p [ind + d]= buf2p [ind];
835    buf3.normalize();
836
837    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
838    i++;
839
840
841    lf= i*d;
842    degfSubLf= degf - lf;
843
844    lg= d*(k-i);
845    deggSubLg= degg - lg;
846
847    buf1p= buf1.rep.elts();
848
849    if (lg >= 0 && deggSubLg > 0)
850    {
851      if (repLengthBuf2 > degfSubLf + 1)
852        degfSubLf= repLengthBuf2 - 1;
853      int tmp= tmin (repLengthBuf1, deggSubLg);
854      for (int ind= 0; ind < tmp; ind++)
855        gp [ind + lg] -= buf1p [ind];
856    }
857    if (lg < 0)
858      break;
859
860    buf2p= buf2.rep.elts();
861    if (degfSubLf >= 0)
862    {
863      for (int ind= 0; ind < repLengthBuf2; ind++)
864        fp [ind + lf] -= buf2p [ind];
865    }
866  }
867
868  return result;
869}
870
871CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
872{
873  Variable y= Variable (2);
874  Variable x= Variable (1);
875
876  zz_pEX f= F;
877  zz_pE *fp= f.rep.elts();
878
879  zz_pEX buf;
880  zz_pE *bufp;
881  CanonicalForm result= 0;
882  int i= 0;
883  int degf= deg(f);
884  int k= 0;
885  int degfSubK;
886  int repLength;
887  while (degf >= k)
888  {
889    degfSubK= degf - k;
890    if (degfSubK >= d)
891      repLength= d;
892    else
893      repLength= degfSubK + 1;
894
895    buf.rep.SetLength ((long) repLength);
896    bufp= buf.rep.elts();
897    for (int j= 0; j < repLength; j++)
898      bufp [j]= fp [j + k];
899    buf.normalize();
900
901    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
902    i++;
903    k= d*i;
904  }
905
906  return result;
907}
908
909CanonicalForm reverseSubstFp (const zz_pX& F, int d)
910{
911  Variable y= Variable (2);
912  Variable x= Variable (1);
913
914  zz_pX f= F;
915  zz_p *fp= f.rep.elts();
916
917  zz_pX buf;
918  zz_p *bufp;
919  CanonicalForm result= 0;
920  int i= 0;
921  int degf= deg(f);
922  int k= 0;
923  int degfSubK;
924  int repLength;
925  while (degf >= k)
926  {
927    degfSubK= degf - k;
928    if (degfSubK >= d)
929      repLength= d;
930    else
931      repLength= degfSubK + 1;
932
933    buf.rep.SetLength ((long) repLength);
934    bufp= buf.rep.elts();
935    for (int j= 0; j < repLength; j++)
936      bufp [j]= fp [j + k];
937    buf.normalize();
938
939    result += convertNTLzzpX2CF (buf, x)*power (y, i);
940    i++;
941    k= d*i;
942  }
943
944  return result;
945}
946
947// assumes input to be reduced mod M and to be an element of Fq not Fp
948CanonicalForm
949mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
950                  CanonicalForm& M)
951{
952  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
953  int d2= tmax (degree (F, 2), degree (G, 2));
954
955  zz_pX F1, F2;
956  kronSubRecipro (F1, F2, F, d1);
957  zz_pX G1, G2;
958  kronSubRecipro (G1, G2, G, d1);
959
960  int k= d1*degree (M);
961  MulTrunc (F1, F1, G1, (long) k);
962
963  mul (F2, F2, G2);
964  if (deg (F2) > k)
965    F2 >>= (k - d1);
966
967  return reverseSubst (F1, F2, d1, d2);
968}
969
970//Kronecker substitution
971CanonicalForm
972mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
973              CanonicalForm& M)
974{
975  CanonicalForm A= F;
976  CanonicalForm B= G;
977
978  int degAx= degree (A, 1);
979  int degAy= degree (A, 2);
980  int degBx= degree (B, 1);
981  int degBy= degree (B, 2);
982  int d1= degAx + 1 + degBx;
983  int d2= tmax (degree (A, 2), degree (B, 2));
984
985  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
986    return mulMod2NTLFpReci (A, B, M);
987
988  zz_pX NTLA= kronSubFp (A, d1);
989  zz_pX NTLB= kronSubFp (B, d1);
990
991  int k= d1*degree (M);
992  MulTrunc (NTLA, NTLA, NTLB, (long) k);
993
994  A= reverseSubstFp (NTLA, d1);
995
996  return A;
997}
998
999// assumes input to be reduced mod M and to be an element of Fq not Fp
1000CanonicalForm
1001mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
1002                  CanonicalForm& M, const Variable& alpha)
1003{
1004  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
1005  int d2= tmax (degree (F, 2), degree (G, 2));
1006
1007  zz_pEX F1, F2;
1008  kronSubRecipro (F1, F2, F, d1, alpha);
1009  zz_pEX G1, G2;
1010  kronSubRecipro (G1, G2, G, d1, alpha);
1011
1012  int k= d1*degree (M);
1013  MulTrunc (F1, F1, G1, (long) k);
1014
1015  mul (F2, F2, G2);
1016  if (deg (F2) > k)
1017    F2 >>= (k-d1);
1018
1019  CanonicalForm result= reverseSubst (F1, F2, d1, d2, alpha);
1020
1021  return result;
1022}
1023
1024CanonicalForm
1025mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
1026              CanonicalForm& M)
1027{
1028  Variable alpha;
1029  CanonicalForm A= F;
1030  CanonicalForm B= G;
1031
1032  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
1033  {
1034    int degAx= degree (A, 1);
1035    int degAy= degree (A, 2);
1036    int degBx= degree (B, 1);
1037    int degBy= degree (B, 2);
1038    int d1= degAx + degBx + 1;
1039    int d2= tmax (degree (A, 2), degree (B, 2));
1040    zz_p::init (getCharacteristic());
1041    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1042    zz_pE::init (NTLMipo);
1043
1044    int degMipo= degree (getMipo (alpha));
1045    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
1046        (2*degAy > degree (M)))
1047      return mulMod2NTLFqReci (A, B, M, alpha);
1048
1049    zz_pEX NTLA= kronSub (A, d1, alpha);
1050    zz_pEX NTLB= kronSub (B, d1, alpha);
1051
1052    int k= d1*degree (M);
1053
1054    MulTrunc (NTLA, NTLA, NTLB, (long) k);
1055
1056    A= reverseSubst (NTLA, d1, alpha);
1057
1058    return A;
1059  }
1060  else
1061    return mulMod2NTLFp (A, B, M);
1062}
1063
1064CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
1065                       const CanonicalForm& M)
1066{
1067  if (A.isZero() || B.isZero())
1068    return 0;
1069
1070  ASSERT (M.isUnivariate(), "M must be univariate");
1071
1072  CanonicalForm F= mod (A, M);
1073  CanonicalForm G= mod (B, M);
1074  if (F.inCoeffDomain() || G.inCoeffDomain())
1075    return F*G;
1076  Variable y= M.mvar();
1077  int degF= degree (F, y);
1078  int degG= degree (G, y);
1079
1080  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
1081      (F.level() == G.level()))
1082  {
1083    CanonicalForm result= mulNTL (F, G);
1084    return mod (result, M);
1085  }
1086  else if (degF <= 1 && degG <= 1)
1087  {
1088    CanonicalForm result= F*G;
1089    return mod (result, M);
1090  }
1091
1092  int sizeF= size (F);
1093  int sizeG= size (G);
1094
1095  int fallBackToNaive= 50;
1096  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
1097    return mod (F*G, M);
1098
1099  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
1100      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
1101    return mulMod2NTLFq (F, G, M);
1102
1103  int m= (int) ceil (degree (M)/2.0);
1104  if (degF >= m || degG >= m)
1105  {
1106    CanonicalForm MLo= power (y, m);
1107    CanonicalForm MHi= power (y, degree (M) - m);
1108    CanonicalForm F0= mod (F, MLo);
1109    CanonicalForm F1= div (F, MLo);
1110    CanonicalForm G0= mod (G, MLo);
1111    CanonicalForm G1= div (G, MLo);
1112    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
1113    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
1114    CanonicalForm F0G0= mulMod2 (F0, G0, M);
1115    return F0G0 + MLo*(F0G1 + F1G0);
1116  }
1117  else
1118  {
1119    m= (int) ceil (tmax (degF, degG)/2.0);
1120    CanonicalForm yToM= power (y, m);
1121    CanonicalForm F0= mod (F, yToM);
1122    CanonicalForm F1= div (F, yToM);
1123    CanonicalForm G0= mod (G, yToM);
1124    CanonicalForm G1= div (G, yToM);
1125    CanonicalForm H00= mulMod2 (F0, G0, M);
1126    CanonicalForm H11= mulMod2 (F1, G1, M);
1127    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
1128    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
1129  }
1130  DEBOUTLN (cerr, "fatal end in mulMod2");
1131}
1132
1133CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
1134                      const CFList& MOD)
1135{
1136  if (A.isZero() || B.isZero())
1137    return 0;
1138
1139  if (MOD.length() == 1)
1140    return mulMod2 (A, B, MOD.getLast());
1141
1142  CanonicalForm M= MOD.getLast();
1143  CanonicalForm F= mod (A, M);
1144  CanonicalForm G= mod (B, M);
1145  if (F.inCoeffDomain() || G.inCoeffDomain())
1146    return F*G;
1147  Variable y= M.mvar();
1148  int degF= degree (F, y);
1149  int degG= degree (G, y);
1150
1151  if ((degF <= 1 && F.level() <= M.level()) &&
1152      (degG <= 1 && G.level() <= M.level()))
1153  {
1154    CFList buf= MOD;
1155    buf.removeLast();
1156    if (degF == 1 && degG == 1)
1157    {
1158      CanonicalForm F0= mod (F, y);
1159      CanonicalForm F1= div (F, y);
1160      CanonicalForm G0= mod (G, y);
1161      CanonicalForm G1= div (G, y);
1162      if (degree (M) > 2)
1163      {
1164        CanonicalForm H00= mulMod (F0, G0, buf);
1165        CanonicalForm H11= mulMod (F1, G1, buf);
1166        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
1167        return H11*y*y + (H01 - H00 - H11)*y + H00;
1168      }
1169      else //here degree (M) == 2
1170      {
1171        buf.append (y);
1172        CanonicalForm F0G1= mulMod (F0, G1, buf);
1173        CanonicalForm F1G0= mulMod (F1, G0, buf);
1174        CanonicalForm F0G0= mulMod (F0, G0, MOD);
1175        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
1176        return result;
1177      }
1178    }
1179    else if (degF == 1 && degG == 0)
1180      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
1181    else if (degF == 0 && degG == 1)
1182      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
1183    else
1184      return mulMod (F, G, buf);
1185  }
1186  int m= (int) ceil (degree (M)/2.0);
1187  if (degF >= m || degG >= m)
1188  {
1189    CanonicalForm MLo= power (y, m);
1190    CanonicalForm MHi= power (y, degree (M) - m);
1191    CanonicalForm F0= mod (F, MLo);
1192    CanonicalForm F1= div (F, MLo);
1193    CanonicalForm G0= mod (G, MLo);
1194    CanonicalForm G1= div (G, MLo);
1195    CFList buf= MOD;
1196    buf.removeLast();
1197    buf.append (MHi);
1198    CanonicalForm F0G1= mulMod (F0, G1, buf);
1199    CanonicalForm F1G0= mulMod (F1, G0, buf);
1200    CanonicalForm F0G0= mulMod (F0, G0, MOD);
1201    return F0G0 + MLo*(F0G1 + F1G0);
1202  }
1203  else
1204  {
1205    m= (int) ceil (tmax (degF, degG)/2.0);
1206    CanonicalForm yToM= power (y, m);
1207    CanonicalForm F0= mod (F, yToM);
1208    CanonicalForm F1= div (F, yToM);
1209    CanonicalForm G0= mod (G, yToM);
1210    CanonicalForm G1= div (G, yToM);
1211    CanonicalForm H00= mulMod (F0, G0, MOD);
1212    CanonicalForm H11= mulMod (F1, G1, MOD);
1213    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
1214    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
1215  }
1216  DEBOUTLN (cerr, "fatal end in mulMod");
1217}
1218
1219CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
1220{
1221  if (L.isEmpty())
1222    return 1;
1223  int l= L.length();
1224  if (l == 1)
1225    return mod (L.getFirst(), M);
1226  else if (l == 2) {
1227    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
1228    return result;
1229  }
1230  else
1231  {
1232    l /= 2;
1233    CFList tmp1, tmp2;
1234    CFListIterator i= L;
1235    CanonicalForm buf1, buf2;
1236    for (int j= 1; j <= l; j++, i++)
1237      tmp1.append (i.getItem());
1238    tmp2= Difference (L, tmp1);
1239    buf1= prodMod (tmp1, M);
1240    buf2= prodMod (tmp2, M);
1241    CanonicalForm result= mulMod2 (buf1, buf2, M);
1242    return result;
1243  }
1244}
1245
1246CanonicalForm prodMod (const CFList& L, const CFList& M)
1247{
1248  if (L.isEmpty())
1249    return 1;
1250  else if (L.length() == 1)
1251    return L.getFirst();
1252  else if (L.length() == 2)
1253    return mulMod (L.getFirst(), L.getLast(), M);
1254  else
1255  {
1256    int l= L.length()/2;
1257    CFListIterator i= L;
1258    CFList tmp1, tmp2;
1259    CanonicalForm buf1, buf2;
1260    for (int j= 1; j <= l; j++, i++)
1261      tmp1.append (i.getItem());
1262    tmp2= Difference (L, tmp1);
1263    buf1= prodMod (tmp1, M);
1264    buf2= prodMod (tmp2, M);
1265    return mulMod (buf1, buf2, M);
1266  }
1267}
1268
1269
1270CanonicalForm reverse (const CanonicalForm& F, int d)
1271{
1272  if (d == 0)
1273    return F;
1274  CanonicalForm A= F;
1275  Variable y= Variable (2);
1276  Variable x= Variable (1);
1277  if (degree (A, x) > 0)
1278  {
1279    A= swapvar (A, x, y);
1280    CanonicalForm result= 0;
1281    CFIterator i= A;
1282    while (d - i.exp() < 0)
1283      i++;
1284
1285    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
1286      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
1287    return result;
1288  }
1289  else
1290    return A*power (x, d);
1291}
1292
1293CanonicalForm
1294newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
1295{
1296  int l= ilog2(n);
1297
1298  CanonicalForm g= mod (F, M)[0] [0];
1299
1300  ASSERT (!g.isZero(), "expected a unit");
1301
1302  Variable alpha;
1303
1304  if (!g.isOne())
1305    g = 1/g;
1306  Variable x= Variable (1);
1307  CanonicalForm result;
1308  int exp= 0;
1309  if (n & 1)
1310  {
1311    result= g;
1312    exp= 1;
1313  }
1314  CanonicalForm h;
1315
1316  for (int i= 1; i <= l; i++)
1317  {
1318    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
1319    h= mod (h, power (x, (1 << i)) - 1);
1320    h= div (h, power (x, (1 << (i - 1))));
1321    h= mod (h, M);
1322    g -= power (x, (1 << (i - 1)))*
1323         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
1324
1325    if (n & (1 << i))
1326    {
1327      if (exp)
1328      {
1329        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
1330        h= mod (h, power (x, exp + (1 << i)) - 1);
1331        h= div (h, power (x, exp));
1332        h= mod (h, M);
1333        result -= power(x, exp)*mod (mulMod2 (g, h, M),
1334                                       power (x, (1 << i)));
1335        exp += (1 << i);
1336      }
1337      else
1338      {
1339        exp= (1 << i);
1340        result= g;
1341      }
1342    }
1343  }
1344
1345  return result;
1346}
1347
1348CanonicalForm
1349newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
1350           M)
1351{
1352  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
1353  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
1354
1355  CanonicalForm A= mod (F, M);
1356  CanonicalForm B= mod (G, M);
1357
1358  Variable x= Variable (1);
1359  int degA= degree (A, x);
1360  int degB= degree (B, x);
1361  int m= degA - degB;
1362  if (m < 0)
1363    return 0;
1364
1365  CanonicalForm Q;
1366  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
1367  {
1368    CanonicalForm R;
1369    divrem2 (A, B, Q, R, M);
1370  }
1371  else
1372  {
1373    CanonicalForm R= reverse (A, degA);
1374    CanonicalForm revB= reverse (B, degB);
1375    revB= newtonInverse (revB, m + 1, M);
1376    Q= mulMod2 (R, revB, M);
1377    Q= mod (Q, power (x, m + 1));
1378    Q= reverse (Q, m);
1379  }
1380
1381  return Q;
1382}
1383
1384void
1385newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1386              CanonicalForm& R, const CanonicalForm& M)
1387{
1388  CanonicalForm A= mod (F, M);
1389  CanonicalForm B= mod (G, M);
1390  Variable x= Variable (1);
1391  int degA= degree (A, x);
1392  int degB= degree (B, x);
1393  int m= degA - degB;
1394
1395  if (m < 0)
1396  {
1397    R= A;
1398    Q= 0;
1399    return;
1400  }
1401
1402  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
1403  {
1404     divrem2 (A, B, Q, R, M);
1405  }
1406  else
1407  {
1408    R= reverse (A, degA);
1409
1410    CanonicalForm revB= reverse (B, degB);
1411    revB= newtonInverse (revB, m + 1, M);
1412    Q= mulMod2 (R, revB, M);
1413
1414    Q= mod (Q, power (x, m + 1));
1415    Q= reverse (Q, m);
1416
1417    R= A - mulMod2 (Q, B, M);
1418  }
1419}
1420
1421static inline
1422CFList split (const CanonicalForm& F, const int m, const Variable& x)
1423{
1424  CanonicalForm A= F;
1425  CanonicalForm buf= 0;
1426  bool swap= false;
1427  if (degree (A, x) <= 0)
1428    return CFList(A);
1429  else if (x.level() != A.level())
1430  {
1431    swap= true;
1432    A= swapvar (A, x, A.mvar());
1433  }
1434
1435  int j= (int) floor ((double) degree (A)/ m);
1436  CFList result;
1437  CFIterator i= A;
1438  for (; j >= 0; j--)
1439  {
1440    while (i.hasTerms() && i.exp() - j*m >= 0)
1441    {
1442      if (swap)
1443        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
1444      else
1445        buf += i.coeff()*power (x, i.exp() - j*m);
1446      i++;
1447    }
1448    if (swap)
1449      result.append (swapvar (buf, x, F.mvar()));
1450    else
1451      result.append (buf);
1452    buf= 0;
1453  }
1454  return result;
1455}
1456
1457static inline
1458void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1459               CanonicalForm& R, const CFList& M);
1460
1461static inline
1462void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1463               CanonicalForm& R, const CFList& M)
1464{
1465  CanonicalForm A= mod (F, M);
1466  CanonicalForm B= mod (G, M);
1467  Variable x= Variable (1);
1468  int degB= degree (B, x);
1469  int degA= degree (A, x);
1470  if (degA < degB)
1471  {
1472    Q= 0;
1473    R= A;
1474    return;
1475  }
1476  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
1477  if (degB < 1)
1478  {
1479    divrem (A, B, Q, R);
1480    Q= mod (Q, M);
1481    R= mod (R, M);
1482    return;
1483  }
1484
1485  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
1486  CFList splitA= split (A, m, x);
1487  if (splitA.length() == 3)
1488    splitA.insert (0);
1489  if (splitA.length() == 2)
1490  {
1491    splitA.insert (0);
1492    splitA.insert (0);
1493  }
1494  if (splitA.length() == 1)
1495  {
1496    splitA.insert (0);
1497    splitA.insert (0);
1498    splitA.insert (0);
1499  }
1500
1501  CanonicalForm xToM= power (x, m);
1502
1503  CFListIterator i= splitA;
1504  CanonicalForm H= i.getItem();
1505  i++;
1506  H *= xToM;
1507  H += i.getItem();
1508  i++;
1509  H *= xToM;
1510  H += i.getItem();
1511  i++;
1512
1513  divrem32 (H, B, Q, R, M);
1514
1515  CFList splitR= split (R, m, x);
1516  if (splitR.length() == 1)
1517    splitR.insert (0);
1518
1519  H= splitR.getFirst();
1520  H *= xToM;
1521  H += splitR.getLast();
1522  H *= xToM;
1523  H += i.getItem();
1524
1525  CanonicalForm bufQ;
1526  divrem32 (H, B, bufQ, R, M);
1527
1528  Q *= xToM;
1529  Q += bufQ;
1530  return;
1531}
1532
1533static inline
1534void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1535               CanonicalForm& R, const CFList& M)
1536{
1537  CanonicalForm A= mod (F, M);
1538  CanonicalForm B= mod (G, M);
1539  Variable x= Variable (1);
1540  int degB= degree (B, x);
1541  int degA= degree (A, x);
1542  if (degA < degB)
1543  {
1544    Q= 0;
1545    R= A;
1546    return;
1547  }
1548  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
1549  if (degB < 1)
1550  {
1551    divrem (A, B, Q, R);
1552    Q= mod (Q, M);
1553    R= mod (R, M);
1554    return;
1555  }
1556  int m= (int) ceil ((double) (degB + 1)/ 2.0);
1557
1558  CFList splitA= split (A, m, x);
1559  CFList splitB= split (B, m, x);
1560
1561  if (splitA.length() == 2)
1562  {
1563    splitA.insert (0);
1564  }
1565  if (splitA.length() == 1)
1566  {
1567    splitA.insert (0);
1568    splitA.insert (0);
1569  }
1570  CanonicalForm xToM= power (x, m);
1571
1572  CanonicalForm H;
1573  CFListIterator i= splitA;
1574  i++;
1575
1576  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
1577  {
1578    H= splitA.getFirst()*xToM + i.getItem();
1579    divrem21 (H, splitB.getFirst(), Q, R, M);
1580  }
1581  else
1582  {
1583    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
1584       splitB.getFirst()*xToM;
1585    Q= xToM - 1;
1586  }
1587
1588  H= mulMod (Q, splitB.getLast(), M);
1589
1590  R= R*xToM + splitA.getLast() - H;
1591
1592  while (degree (R, x) >= degB)
1593  {
1594    xToM= power (x, degree (R, x) - degB);
1595    Q += LC (R, x)*xToM;
1596    R -= mulMod (LC (R, x), B, M)*xToM;
1597    Q= mod (Q, M);
1598    R= mod (R, M);
1599  }
1600
1601  return;
1602}
1603
1604void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1605              CanonicalForm& R, const CanonicalForm& M)
1606{
1607  CanonicalForm A= mod (F, M);
1608  CanonicalForm B= mod (G, M);
1609
1610  if (B.inCoeffDomain())
1611  {
1612    divrem (A, B, Q, R);
1613    return;
1614  }
1615  if (A.inCoeffDomain() && !B.inCoeffDomain())
1616  {
1617    Q= 0;
1618    R= A;
1619    return;
1620  }
1621
1622  if (B.level() < A.level())
1623  {
1624    divrem (A, B, Q, R);
1625    return;
1626  }
1627  if (A.level() > B.level())
1628  {
1629    R= A;
1630    Q= 0;
1631    return;
1632  }
1633  if (B.level() == 1 && B.isUnivariate())
1634  {
1635    divrem (A, B, Q, R);
1636    return;
1637  }
1638  if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
1639  {
1640    Q= 0;
1641    R= A;
1642    return;
1643  }
1644
1645  Variable x= Variable (1);
1646  int degB= degree (B, x);
1647  if (degB > degree (A, x))
1648  {
1649    Q= 0;
1650    R= A;
1651    return;
1652  }
1653
1654  CFList splitA= split (A, degB, x);
1655
1656  CanonicalForm xToDegB= power (x, degB);
1657  CanonicalForm H, bufQ;
1658  Q= 0;
1659  CFListIterator i= splitA;
1660  H= i.getItem()*xToDegB;
1661  i++;
1662  H += i.getItem();
1663  CFList buf;
1664  while (i.hasItem())
1665  {
1666    buf= CFList (M);
1667    divrem21 (H, B, bufQ, R, buf);
1668    i++;
1669    if (i.hasItem())
1670      H= R*xToDegB + i.getItem();
1671    Q *= xToDegB;
1672    Q += bufQ;
1673  }
1674  return;
1675}
1676
1677void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
1678             CanonicalForm& R, const CFList& MOD)
1679{
1680  CanonicalForm A= mod (F, MOD);
1681  CanonicalForm B= mod (G, MOD);
1682  Variable x= Variable (1);
1683  int degB= degree (B, x);
1684  if (degB > degree (A, x))
1685  {
1686    Q= 0;
1687    R= A;
1688    return;
1689  }
1690
1691  if (degB <= 0)
1692  {
1693    divrem (A, B, Q, R);
1694    Q= mod (Q, MOD);
1695    R= mod (R, MOD);
1696    return;
1697  }
1698  CFList splitA= split (A, degB, x);
1699
1700  CanonicalForm xToDegB= power (x, degB);
1701  CanonicalForm H, bufQ;
1702  Q= 0;
1703  CFListIterator i= splitA;
1704  H= i.getItem()*xToDegB;
1705  i++;
1706  H += i.getItem();
1707  while (i.hasItem())
1708  {
1709    divrem21 (H, B, bufQ, R, MOD);
1710    i++;
1711    if (i.hasItem())
1712      H= R*xToDegB + i.getItem();
1713    Q *= xToDegB;
1714    Q += bufQ;
1715  }
1716  return;
1717}
1718
1719void sortList (CFList& list, const Variable& x)
1720{
1721  int l= 1;
1722  int k= 1;
1723  CanonicalForm buf;
1724  CFListIterator m;
1725  for (CFListIterator i= list; l <= list.length(); i++, l++)
1726  {
1727    for (CFListIterator j= list; k <= list.length() - l; k++)
1728    {
1729      m= j;
1730      m++;
1731      if (degree (j.getItem(), x) > degree (m.getItem(), x))
1732      {
1733        buf= m.getItem();
1734        m.getItem()= j.getItem();
1735        j.getItem()= buf;
1736        j++;
1737        j.getItem()= m.getItem();
1738      }
1739      else
1740        j++;
1741    }
1742    k= 1;
1743  }
1744}
1745
1746static inline
1747CFList diophantine (const CanonicalForm& F, const CFList& factors)
1748{
1749  if (getCharacteristic() == 0)
1750  {
1751    Variable v;
1752    bool hasAlgVar= hasFirstAlgVar (F, v);
1753    for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++)
1754      hasAlgVar= hasFirstAlgVar (i.getItem(), v);
1755    if (hasAlgVar)
1756    {
1757      CFList result= modularDiophant (F, factors, getMipo (v));
1758      return result;
1759    }
1760  }
1761
1762  CanonicalForm buf1, buf2, buf3, S, T;
1763  CFListIterator i= factors;
1764  CFList result;
1765  if (i.hasItem())
1766    i++;
1767  buf1= F/factors.getFirst();
1768  buf2= divNTL (F, i.getItem());
1769  buf3= extgcd (buf1, buf2, S, T);
1770  result.append (S);
1771  result.append (T);
1772  if (i.hasItem())
1773    i++;
1774  for (; i.hasItem(); i++)
1775  {
1776    buf1= divNTL (F, i.getItem());
1777    buf3= extgcd (buf3, buf1, S, T);
1778    CFListIterator k= factors;
1779    for (CFListIterator j= result; j.hasItem(); j++, k++)
1780    {
1781      j.getItem()= mulNTL (j.getItem(), S);
1782      j.getItem()= modNTL (j.getItem(), k.getItem());
1783    }
1784    result.append (T);
1785  }
1786  return result;
1787}
1788
1789void
1790henselStep12 (const CanonicalForm& F, const CFList& factors,
1791              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1792              CFArray& Pi, int j)
1793{
1794  CanonicalForm E;
1795  CanonicalForm xToJ= power (F.mvar(), j);
1796  Variable x= F.mvar();
1797  // compute the error
1798  if (j == 1)
1799    E= F[j];
1800  else
1801  {
1802    if (degree (Pi [factors.length() - 2], x) > 0)
1803      E= F[j] - Pi [factors.length() - 2] [j];
1804    else
1805      E= F[j];
1806  }
1807
1808  CFArray buf= CFArray (diophant.length());
1809  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1810  int k= 0;
1811  CanonicalForm remainder;
1812  // actual lifting
1813  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1814  {
1815    if (degree (bufFactors[k], x) > 0)
1816    {
1817      if (k > 0)
1818        remainder= modNTL (E, bufFactors[k] [0]);
1819      else
1820        remainder= E;
1821    }
1822    else
1823      remainder= modNTL (E, bufFactors[k]);
1824
1825    buf[k]= mulNTL (i.getItem(), remainder);
1826    if (degree (bufFactors[k], x) > 0)
1827      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
1828    else
1829      buf[k]= modNTL (buf[k], bufFactors[k]);
1830  }
1831  for (k= 1; k < factors.length(); k++)
1832    bufFactors[k] += xToJ*buf[k];
1833
1834  // update Pi [0]
1835  int degBuf0= degree (bufFactors[0], x);
1836  int degBuf1= degree (bufFactors[1], x);
1837  if (degBuf0 > 0 && degBuf1 > 0)
1838    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
1839  CanonicalForm uIZeroJ;
1840  if (j == 1)
1841  {
1842    if (degBuf0 > 0 && degBuf1 > 0)
1843      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
1844                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
1845    else if (degBuf0 > 0)
1846      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
1847    else if (degBuf1 > 0)
1848      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
1849    else
1850      uIZeroJ= 0;
1851    Pi [0] += xToJ*uIZeroJ;
1852  }
1853  else
1854  {
1855    if (degBuf0 > 0 && degBuf1 > 0)
1856      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
1857                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
1858    else if (degBuf0 > 0)
1859      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
1860    else if (degBuf1 > 0)
1861      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
1862    else
1863      uIZeroJ= 0;
1864    Pi [0] += xToJ*uIZeroJ;
1865  }
1866  CFArray tmp= CFArray (factors.length() - 1);
1867  for (k= 0; k < factors.length() - 1; k++)
1868    tmp[k]= 0;
1869  CFIterator one, two;
1870  one= bufFactors [0];
1871  two= bufFactors [1];
1872  if (degBuf0 > 0 && degBuf1 > 0)
1873  {
1874    for (k= 1; k <= (int) ceil (j/2.0); k++)
1875    {
1876      if (k != j - k + 1)
1877      {
1878        if ((one.hasTerms() && one.exp() == j - k + 1) && (two.hasTerms() && two.exp() == j - k + 1))
1879        {
1880          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
1881                     two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
1882          one++;
1883          two++;
1884        }
1885        else if (one.hasTerms() && one.exp() == j - k + 1)
1886        {
1887          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) -
1888                     M (k + 1, 1);
1889          one++;
1890        }
1891        else if (two.hasTerms() && two.exp() == j - k + 1)
1892        {
1893          tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) -
1894                    M (k + 1, 1);
1895          two++;
1896        }
1897      }
1898      else
1899      {
1900        tmp[0] += M (k + 1, 1);
1901      }
1902    }
1903  }
1904  Pi [0] += tmp[0]*xToJ*F.mvar();
1905
1906  // update Pi [l]
1907  int degPi, degBuf;
1908  for (int l= 1; l < factors.length() - 1; l++)
1909  {
1910    degPi= degree (Pi [l - 1], x);
1911    degBuf= degree (bufFactors[l + 1], x);
1912    if (degPi > 0 && degBuf > 0)
1913      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
1914    if (j == 1)
1915    {
1916      if (degPi > 0 && degBuf > 0)
1917        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
1918                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
1919                  M (1, l + 1));
1920      else if (degPi > 0)
1921        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
1922      else if (degBuf > 0)
1923        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
1924    }
1925    else
1926    {
1927      if (degPi > 0 && degBuf > 0)
1928      {
1929        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
1930        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
1931      }
1932      else if (degPi > 0)
1933        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
1934      else if (degBuf > 0)
1935      {
1936        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
1937        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
1938      }
1939      Pi[l] += xToJ*uIZeroJ;
1940    }
1941    one= bufFactors [l + 1];
1942    two= Pi [l - 1];
1943    if (two.hasTerms() && two.exp() == j + 1)
1944    {
1945      if (degBuf > 0 && degPi > 0)
1946      {
1947          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
1948          two++;
1949      }
1950      else if (degPi > 0)
1951      {
1952          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
1953          two++;
1954      }
1955    }
1956    if (degBuf > 0 && degPi > 0)
1957    {
1958      for (k= 1; k <= (int) ceil (j/2.0); k++)
1959      {
1960        if (k != j - k + 1)
1961        {
1962          if ((one.hasTerms() && one.exp() == j - k + 1) &&
1963              (two.hasTerms() && two.exp() == j - k + 1))
1964          {
1965            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
1966                       two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
1967            one++;
1968            two++;
1969          }
1970          else if (one.hasTerms() && one.exp() == j - k + 1)
1971          {
1972            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) -
1973                       M (k + 1, l + 1);
1974            one++;
1975          }
1976          else if (two.hasTerms() && two.exp() == j - k + 1)
1977          {
1978            tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) -
1979                      M (k + 1, l + 1);
1980            two++;
1981          }
1982        }
1983        else
1984          tmp[l] += M (k + 1, l + 1);
1985      }
1986    }
1987    Pi[l] += tmp[l]*xToJ*F.mvar();
1988  }
1989  return;
1990}
1991
1992void
1993henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1994              CFList& diophant, CFMatrix& M, bool sort)
1995{
1996  if (sort)
1997    sortList (factors, Variable (1));
1998  Pi= CFArray (factors.length() - 1);
1999  CFListIterator j= factors;
2000  diophant= diophantine (F[0], factors);
2001  DEBOUTLN (cerr, "diophant= " << diophant);
2002  j++;
2003  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()));
2004  M (1, 1)= Pi [0];
2005  int i= 1;
2006  if (j.hasItem())
2007    j++;
2008  for (; j.hasItem(); j++, i++)
2009  {
2010    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
2011    M (1, i + 1)= Pi [i];
2012  }
2013  CFArray bufFactors= CFArray (factors.length());
2014  i= 0;
2015  for (CFListIterator k= factors; k.hasItem(); i++, k++)
2016  {
2017    if (i == 0)
2018      bufFactors[i]= mod (k.getItem(), F.mvar());
2019    else
2020      bufFactors[i]= k.getItem();
2021  }
2022  for (i= 1; i < l; i++)
2023    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
2024
2025  CFListIterator k= factors;
2026  for (i= 0; i < factors.length (); i++, k++)
2027    k.getItem()= bufFactors[i];
2028  factors.removeFirst();
2029  return;
2030}
2031
2032void
2033henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
2034                    end, CFArray& Pi, const CFList& diophant, CFMatrix& M)
2035{
2036  CFArray bufFactors= CFArray (factors.length());
2037  int i= 0;
2038  CanonicalForm xToStart= power (F.mvar(), start);
2039  for (CFListIterator k= factors; k.hasItem(); k++, i++)
2040  {
2041    if (i == 0)
2042      bufFactors[i]= mod (k.getItem(), xToStart);
2043    else
2044      bufFactors[i]= k.getItem();
2045  }
2046  for (i= start; i < end; i++)
2047    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
2048
2049  CFListIterator k= factors;
2050  for (i= 0; i < factors.length(); k++, i++)
2051    k.getItem()= bufFactors [i];
2052  factors.removeFirst();
2053  return;
2054}
2055
2056static inline
2057CFList
2058biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
2059{
2060  Variable y= F.mvar();
2061  CFList result;
2062  if (y.level() == 1)
2063  {
2064    result= diophantine (F, factors);
2065    return result;
2066  }
2067  else
2068  {
2069    CFList buf= factors;
2070    for (CFListIterator i= buf; i.hasItem(); i++)
2071      i.getItem()= mod (i.getItem(), y);
2072    CanonicalForm A= mod (F, y);
2073    int bufD= 1;
2074    CFList recResult= biDiophantine (A, buf, bufD);
2075    CanonicalForm e= 1;
2076    CFList p;
2077    CFArray bufFactors= CFArray (factors.length());
2078    CanonicalForm yToD= power (y, d);
2079    int k= 0;
2080    for (CFListIterator i= factors; i.hasItem(); i++, k++)
2081    {
2082      bufFactors [k]= i.getItem();
2083    }
2084    CanonicalForm b, quot;
2085    for (k= 0; k < factors.length(); k++) //TODO compute b's faster
2086    {
2087      b= 1;
2088      if (fdivides (bufFactors[k], F, quot))
2089        b= quot;
2090      else
2091      {
2092        for (int l= 0; l < factors.length(); l++)
2093        {
2094          if (l == k)
2095            continue;
2096          else
2097          {
2098            b= mulMod2 (b, bufFactors[l], yToD);
2099          }
2100        }
2101      }
2102      p.append (b);
2103    }
2104
2105    CFListIterator j= p;
2106    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2107      e -= i.getItem()*j.getItem();
2108
2109    if (e.isZero())
2110      return recResult;
2111    CanonicalForm coeffE;
2112    CFList s;
2113    result= recResult;
2114    CanonicalForm g;
2115    for (int i= 1; i < d; i++)
2116    {
2117      if (degree (e, y) > 0)
2118        coeffE= e[i];
2119      else
2120        coeffE= 0;
2121      if (!coeffE.isZero())
2122      {
2123        CFListIterator k= result;
2124        CFListIterator l= p;
2125        int ii= 0;
2126        j= recResult;
2127        for (; j.hasItem(); j++, k++, l++, ii++)
2128        {
2129          g= coeffE*j.getItem();
2130          if (degree (bufFactors[ii], y) <= 0)
2131            g= mod (g, bufFactors[ii]);
2132          else
2133            g= mod (g, bufFactors[ii][0]);
2134          k.getItem() += g*power (y, i);
2135          e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
2136          DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
2137                    mod (e, power (y, i + 1)));
2138        }
2139      }
2140      if (e.isZero())
2141        break;
2142    }
2143
2144    DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
2145
2146#ifdef DEBUGOUTPUT
2147    CanonicalForm test= 0;
2148    j= p;
2149    for (CFListIterator i= result; i.hasItem(); i++, j++)
2150      test += mod (i.getItem()*j.getItem(), power (y, d));
2151    DEBOUTLN (cerr, "test= " << test);
2152#endif
2153    return result;
2154  }
2155}
2156
2157static inline
2158CFList
2159multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
2160                     const CFList& recResult, const CFList& M, const int d)
2161{
2162  Variable y= F.mvar();
2163  CFList result;
2164  CFListIterator i;
2165  CanonicalForm e= 1;
2166  CFListIterator j= factors;
2167  CFList p;
2168  CFArray bufFactors= CFArray (factors.length());
2169  CanonicalForm yToD= power (y, d);
2170  int k= 0;
2171  for (CFListIterator i= factors; i.hasItem(); i++, k++)
2172    bufFactors [k]= i.getItem();
2173  CanonicalForm b, quot;
2174  CFList buf= M;
2175  buf.removeLast();
2176  buf.append (yToD);
2177  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
2178  {
2179    b= 1;
2180    if (fdivides (bufFactors[k], F, quot))
2181      b= quot;
2182    else
2183    {
2184      for (int l= 0; l < factors.length(); l++)
2185      {
2186        if (l == k)
2187          continue;
2188        else
2189        {
2190          b= mulMod (b, bufFactors[l], buf);
2191        }
2192      }
2193    }
2194    p.append (b);
2195  }
2196  j= p;
2197  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2198    e -= mulMod (i.getItem(), j.getItem(), M);
2199
2200  if (e.isZero())
2201    return recResult;
2202  CanonicalForm coeffE;
2203  CFList s;
2204  result= recResult;
2205  CanonicalForm g;
2206  for (int i= 1; i < d; i++)
2207  {
2208    if (degree (e, y) > 0)
2209      coeffE= e[i];
2210    else
2211      coeffE= 0;
2212    if (!coeffE.isZero())
2213    {
2214      CFListIterator k= result;
2215      CFListIterator l= p;
2216      j= recResult;
2217      int ii= 0;
2218      CanonicalForm dummy;
2219      for (; j.hasItem(); j++, k++, l++, ii++)
2220      {
2221        g= mulMod (coeffE, j.getItem(), M);
2222        if (degree (bufFactors[ii], y) <= 0)
2223          divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
2224                  g, M);
2225        else
2226          divrem (g, bufFactors[ii][0], dummy, g, M);
2227        k.getItem() += g*power (y, i);
2228        e -= mulMod (g*power (y, i), l.getItem(), M);
2229      }
2230    }
2231
2232    if (e.isZero())
2233      break;
2234  }
2235
2236#ifdef DEBUGOUTPUT
2237    CanonicalForm test= 0;
2238    j= p;
2239    for (CFListIterator i= result; i.hasItem(); i++, j++)
2240      test += mod (i.getItem()*j.getItem(), power (y, d));
2241    DEBOUTLN (cerr, "test= " << test);
2242#endif
2243  return result;
2244}
2245
2246static inline
2247void
2248henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
2249            const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
2250            const CFList& MOD)
2251{
2252  CanonicalForm E;
2253  CanonicalForm xToJ= power (F.mvar(), j);
2254  Variable x= F.mvar();
2255  // compute the error
2256  if (j == 1)
2257  {
2258    E= F[j];
2259#ifdef DEBUGOUTPUT
2260    CanonicalForm test= 1;
2261    for (int i= 0; i < factors.length(); i++)
2262    {
2263      if (i == 0)
2264        test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
2265      else
2266        test= mulMod (test, bufFactors[i], MOD);
2267    }
2268    CanonicalForm test2= mod (F-test, xToJ);
2269
2270    test2= mod (test2, MOD);
2271    DEBOUTLN (cerr, "test= " << test2);
2272#endif
2273  }
2274  else
2275  {
2276#ifdef DEBUGOUTPUT
2277    CanonicalForm test= 1;
2278    for (int i= 0; i < factors.length(); i++)
2279    {
2280      if (i == 0)
2281        test *= mod (bufFactors [i], power (x, j));
2282      else
2283        test *= bufFactors[i];
2284    }
2285    test= mod (test, power (x, j));
2286    test= mod (test, MOD);
2287    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
2288    DEBOUTLN (cerr, "test= " << test2);
2289#endif
2290
2291    if (degree (Pi [factors.length() - 2], x) > 0)
2292      E= F[j] - Pi [factors.length() - 2] [j];
2293    else
2294      E= F[j];
2295  }
2296
2297  CFArray buf= CFArray (diophant.length());
2298  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
2299  int k= 0;
2300  // actual lifting
2301  CanonicalForm dummy, rest1;
2302  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
2303  {
2304    if (degree (bufFactors[k], x) > 0)
2305    {
2306      if (k > 0)
2307        divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
2308      else
2309        rest1= E;
2310    }
2311    else
2312      divrem (E, bufFactors[k], dummy, rest1, MOD);
2313
2314    buf[k]= mulMod (i.getItem(), rest1, MOD);
2315
2316    if (degree (bufFactors[k], x) > 0)
2317      divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
2318    else
2319      divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
2320  }
2321  for (k= 1; k < factors.length(); k++)
2322    bufFactors[k] += xToJ*buf[k];
2323
2324  // update Pi [0]
2325  int degBuf0= degree (bufFactors[0], x);
2326  int degBuf1= degree (bufFactors[1], x);
2327  if (degBuf0 > 0 && degBuf1 > 0)
2328    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
2329  CanonicalForm uIZeroJ;
2330  if (j == 1)
2331  {
2332    if (degBuf0 > 0 && degBuf1 > 0)
2333      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
2334                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
2335    else if (degBuf0 > 0)
2336      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
2337    else if (degBuf1 > 0)
2338      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
2339    else
2340      uIZeroJ= 0;
2341    Pi [0] += xToJ*uIZeroJ;
2342  }
2343  else
2344  {
2345    if (degBuf0 > 0 && degBuf1 > 0)
2346      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
2347                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
2348    else if (degBuf0 > 0)
2349      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
2350    else if (degBuf1 > 0)
2351      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
2352    else
2353      uIZeroJ= 0;
2354    Pi [0] += xToJ*uIZeroJ;
2355  }
2356
2357  CFArray tmp= CFArray (factors.length() - 1);
2358  for (k= 0; k < factors.length() - 1; k++)
2359    tmp[k]= 0;
2360  CFIterator one, two;
2361  one= bufFactors [0];
2362  two= bufFactors [1];
2363  if (degBuf0 > 0 && degBuf1 > 0)
2364  {
2365    for (k= 1; k <= (int) ceil (j/2.0); k++)
2366    {
2367      if (k != j - k + 1)
2368      {
2369        if ((one.hasTerms() && one.exp() == j - k + 1) &&
2370            (two.hasTerms() && two.exp() == j - k + 1))
2371        {
2372          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2373                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
2374                    M (j - k + 2, 1);
2375          one++;
2376          two++;
2377        }
2378        else if (one.hasTerms() && one.exp() == j - k + 1)
2379        {
2380          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2381                    bufFactors[1] [k], MOD) - M (k + 1, 1);
2382          one++;
2383        }
2384        else if (two.hasTerms() && two.exp() == j - k + 1)
2385        {
2386          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
2387                    two.coeff()), MOD) - M (k + 1, 1);
2388          two++;
2389        }
2390      }
2391      else
2392      {
2393        tmp[0] += M (k + 1, 1);
2394      }
2395    }
2396  }
2397  Pi [0] += tmp[0]*xToJ*F.mvar();
2398
2399  // update Pi [l]
2400  int degPi, degBuf;
2401  for (int l= 1; l < factors.length() - 1; l++)
2402  {
2403    degPi= degree (Pi [l - 1], x);
2404    degBuf= degree (bufFactors[l + 1], x);
2405    if (degPi > 0 && degBuf > 0)
2406      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
2407    if (j == 1)
2408    {
2409      if (degPi > 0 && degBuf > 0)
2410        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
2411                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
2412                  M (1, l + 1));
2413      else if (degPi > 0)
2414        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
2415      else if (degBuf > 0)
2416        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
2417    }
2418    else
2419    {
2420      if (degPi > 0 && degBuf > 0)
2421      {
2422        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
2423        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
2424      }
2425      else if (degPi > 0)
2426        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
2427      else if (degBuf > 0)
2428      {
2429        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
2430        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
2431      }
2432      Pi[l] += xToJ*uIZeroJ;
2433    }
2434    one= bufFactors [l + 1];
2435    two= Pi [l - 1];
2436    if (two.hasTerms() && two.exp() == j + 1)
2437    {
2438      if (degBuf > 0 && degPi > 0)
2439      {
2440          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
2441          two++;
2442      }
2443      else if (degPi > 0)
2444      {
2445          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
2446          two++;
2447      }
2448    }
2449    if (degBuf > 0 && degPi > 0)
2450    {
2451      for (k= 1; k <= (int) ceil (j/2.0); k++)
2452      {
2453        if (k != j - k + 1)
2454        {
2455          if ((one.hasTerms() && one.exp() == j - k + 1) &&
2456              (two.hasTerms() && two.exp() == j - k + 1))
2457          {
2458            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2459                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
2460                      M (j - k + 2, l + 1);
2461            one++;
2462            two++;
2463          }
2464          else if (one.hasTerms() && one.exp() == j - k + 1)
2465          {
2466            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2467                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
2468            one++;
2469          }
2470          else if (two.hasTerms() && two.exp() == j - k + 1)
2471          {
2472            tmp[l] += mulMod (bufFactors[l + 1] [k],
2473                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
2474            two++;
2475          }
2476        }
2477        else
2478          tmp[l] += M (k + 1, l + 1);
2479      }
2480    }
2481    Pi[l] += tmp[l]*xToJ*F.mvar();
2482  }
2483
2484  return;
2485}
2486
2487CFList
2488henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
2489              diophant, CFArray& Pi, CFMatrix& M)
2490{
2491  CFList buf= factors;
2492  int k= 0;
2493  int liftBoundBivar= l[k];
2494  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
2495  CFList MOD;
2496  MOD.append (power (Variable (2), liftBoundBivar));
2497  CFArray bufFactors= CFArray (factors.length());
2498  k= 0;
2499  CFListIterator j= eval;
2500  j++;
2501  buf.removeFirst();
2502  buf.insert (LC (j.getItem(), 1));
2503  for (CFListIterator i= buf; i.hasItem(); i++, k++)
2504    bufFactors[k]= i.getItem();
2505  Pi= CFArray (factors.length() - 1);
2506  CFListIterator i= buf;
2507  i++;
2508  Variable y= j.getItem().mvar();
2509  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
2510  M (1, 1)= Pi [0];
2511  k= 1;
2512  if (i.hasItem())
2513    i++;
2514  for (; i.hasItem(); i++, k++)
2515  {
2516    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
2517    M (1, k + 1)= Pi [k];
2518  }
2519
2520  for (int d= 1; d < l[1]; d++)
2521    henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
2522  CFList result;
2523  for (k= 1; k < factors.length(); k++)
2524    result.append (bufFactors[k]);
2525  return result;
2526}
2527
2528void
2529henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
2530                  CFArray& Pi, const CFList& diophant, CFMatrix& M,
2531                  const CFList& MOD)
2532{
2533  CFArray bufFactors= CFArray (factors.length());
2534  int i= 0;
2535  CanonicalForm xToStart= power (F.mvar(), start);
2536  for (CFListIterator k= factors; k.hasItem(); k++, i++)
2537  {
2538    if (i == 0)
2539      bufFactors[i]= mod (k.getItem(), xToStart);
2540    else
2541      bufFactors[i]= k.getItem();
2542  }
2543  for (i= start; i < end; i++)
2544    henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD);
2545
2546  CFListIterator k= factors;
2547  for (i= 0; i < factors.length(); k++, i++)
2548    k.getItem()= bufFactors [i];
2549  factors.removeFirst();
2550  return;
2551}
2552
2553CFList
2554henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
2555            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
2556            lNew)
2557{
2558  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
2559  int k= 0;
2560  CFArray bufFactors= CFArray (factors.length());
2561  for (CFListIterator i= factors; i.hasItem(); i++, k++)
2562  {
2563    if (k == 0)
2564      bufFactors[k]= LC (F.getLast(), 1);
2565    else
2566      bufFactors[k]= i.getItem();
2567  }
2568  CFList buf= factors;
2569  buf.removeFirst();
2570  buf.insert (LC (F.getLast(), 1));
2571  CFListIterator i= buf;
2572  i++;
2573  Variable y= F.getLast().mvar();
2574  Variable x= F.getFirst().mvar();
2575  CanonicalForm xToLOld= power (x, lOld);
2576  Pi [0]= mod (Pi[0], xToLOld);
2577  M (1, 1)= Pi [0];
2578  k= 1;
2579  if (i.hasItem())
2580    i++;
2581  for (; i.hasItem(); i++, k++)
2582  {
2583    Pi [k]= mod (Pi [k], xToLOld);
2584    M (1, k + 1)= Pi [k];
2585  }
2586
2587  for (int d= 1; d < lNew; d++)
2588    henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
2589  CFList result;
2590  for (k= 1; k < factors.length(); k++)
2591    result.append (bufFactors[k]);
2592  return result;
2593}
2594
2595CFList
2596henselLift (const CFList& eval, const CFList& factors, const int* l, const int
2597            lLength, bool sort)
2598{
2599  CFList diophant;
2600  CFList buf= factors;
2601  buf.insert (LC (eval.getFirst(), 1));
2602  if (sort)
2603    sortList (buf, Variable (1));
2604  CFArray Pi;
2605  CFMatrix M= CFMatrix (l[1], factors.length());
2606  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
2607  if (eval.length() == 2)
2608    return result;
2609  CFList MOD;
2610  for (int i= 0; i < 2; i++)
2611    MOD.append (power (Variable (i + 2), l[i]));
2612  CFListIterator j= eval;
2613  j++;
2614  CFList bufEval;
2615  bufEval.append (j.getItem());
2616  j++;
2617
2618  for (int i= 2; i < lLength && j.hasItem(); i++, j++)
2619  {
2620    result.insert (LC (bufEval.getFirst(), 1));
2621    bufEval.append (j.getItem());
2622    M= CFMatrix (l[i], factors.length());
2623    result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
2624    MOD.append (power (Variable (i + 2), l[i]));
2625    bufEval.removeFirst();
2626  }
2627  return result;
2628}
2629
2630void
2631henselStep122 (const CanonicalForm& F, const CFList& factors,
2632              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
2633              CFArray& Pi, int j, const CFArray& LCs)
2634{
2635  Variable x= F.mvar();
2636  CanonicalForm xToJ= power (x, j);
2637
2638  CanonicalForm E;
2639  // compute the error
2640  if (degree (Pi [factors.length() - 2], x) > 0)
2641    E= F[j] - Pi [factors.length() - 2] [j];
2642  else
2643    E= F[j];
2644
2645  CFArray buf= CFArray (diophant.length());
2646
2647  int k= 0;
2648  CanonicalForm remainder;
2649  // actual lifting
2650  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
2651  {
2652    if (degree (bufFactors[k], x) > 0)
2653      remainder= modNTL (E, bufFactors[k] [0]);
2654    else
2655      remainder= modNTL (E, bufFactors[k]);
2656    buf[k]= mulNTL (i.getItem(), remainder);
2657    if (degree (bufFactors[k], x) > 0)
2658      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
2659    else
2660      buf[k]= modNTL (buf[k], bufFactors[k]);
2661  }
2662
2663  for (k= 0; k < factors.length(); k++)
2664    bufFactors[k] += xToJ*buf[k];
2665
2666  // update Pi [0]
2667  int degBuf0= degree (bufFactors[0], x);
2668  int degBuf1= degree (bufFactors[1], x);
2669  if (degBuf0 > 0 && degBuf1 > 0)
2670  {
2671    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
2672    if (j + 2 <= M.rows())
2673      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
2674  }
2675  CanonicalForm uIZeroJ;
2676  if (degBuf0 > 0 && degBuf1 > 0)
2677    uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]);
2678  else if (degBuf0 > 0)
2679    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
2680  else if (degBuf1 > 0)
2681    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
2682  else
2683    uIZeroJ= 0;
2684  Pi [0] += xToJ*uIZeroJ;
2685
2686  CFArray tmp= CFArray (factors.length() - 1);
2687  for (k= 0; k < factors.length() - 1; k++)
2688    tmp[k]= 0;
2689  CFIterator one, two;
2690  one= bufFactors [0];
2691  two= bufFactors [1];
2692  if (degBuf0 > 0 && degBuf1 > 0)
2693  {
2694    while (one.hasTerms() && one.exp() > j) one++;
2695    while (two.hasTerms() && two.exp() > j) two++;
2696    for (k= 1; k <= (int) ceil (j/2.0); k++)
2697    {
2698      if (one.hasTerms() && two.hasTerms())
2699      {
2700        if (k != j - k + 1)
2701        {
2702          if ((one.hasTerms() && one.exp() == j - k + 1) && +
2703              (two.hasTerms() && two.exp() == j - k + 1))
2704        {
2705          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
2706                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
2707          one++;
2708          two++;
2709        }
2710        else if (one.hasTerms() && one.exp() == j - k + 1)
2711        {
2712          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
2713                      M (k + 1, 1);
2714          one++;
2715        }
2716        else if (two.hasTerms() && two.exp() == j - k + 1)
2717        {
2718          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
2719                    M (k + 1, 1);
2720          two++;
2721        }
2722      }
2723      else
2724        tmp[0] += M (k + 1, 1);
2725      }
2726    }
2727  }
2728
2729  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
2730  {
2731    if (j + 2 <= M.rows())
2732      tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
2733                         (bufFactors [1] [j + 1] + bufFactors [1] [0]))
2734                         - M(1,1) - M (j + 2,1);
2735  }
2736  else if (degBuf0 >= j + 1)
2737  {
2738    if (degBuf1 > 0)
2739      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
2740    else
2741      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]);
2742  }
2743  else if (degBuf1 >= j + 1)
2744  {
2745    if (degBuf0 > 0)
2746      tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
2747    else
2748      tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]);
2749  }
2750
2751  Pi [0] += tmp[0]*xToJ*F.mvar();
2752  return;
2753}
2754
2755void
2756henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
2757              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
2758{
2759  if (sort)
2760    sortList (factors, Variable (1));
2761  Pi= CFArray (factors.length() - 2);
2762  CFList bufFactors2= factors;
2763  bufFactors2.removeFirst();
2764  CanonicalForm s,t;
2765  extgcd (bufFactors2.getFirst(), bufFactors2.getLast(), s, t);
2766  diophant= CFList();
2767  diophant.append (t);
2768  diophant.append (s);
2769  DEBOUTLN (cerr, "diophant= " << diophant);
2770
2771  CFArray bufFactors= CFArray (bufFactors2.length());
2772  int i= 0;
2773  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
2774    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
2775
2776  Variable x= F.mvar();
2777  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
2778  {
2779    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
2780    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
2781                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
2782  }
2783  else if (degree (bufFactors[0], x) > 0)
2784  {
2785    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
2786    Pi [0]= M (1, 1) +
2787            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
2788  }
2789  else if (degree (bufFactors[1], x) > 0)
2790  {
2791    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
2792    Pi [0]= M (1, 1) +
2793            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
2794  }
2795  else
2796  {
2797    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
2798    Pi [0]= M (1, 1);
2799  }
2800
2801  for (i= 1; i < l; i++)
2802    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
2803
2804  factors= CFList();
2805  for (i= 0; i < bufFactors.size(); i++)
2806    factors.append (bufFactors[i]);
2807  return;
2808}
2809
2810/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
2811/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
2812static inline
2813CFList
2814diophantine (const CFList& recResult, const CFList& factors,
2815             const CFList& products, const CFList& M, const CanonicalForm& E,
2816             bool& bad)
2817{
2818  if (M.isEmpty())
2819  {
2820    CFList result;
2821    CFListIterator j= factors;
2822    CanonicalForm buf;
2823    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2824    {
2825      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
2826              "constant or univariate poly expected");
2827      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
2828              "constant or univariate poly expected");
2829      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
2830              "constant or univariate poly expected");
2831      buf= mulNTL (E, i.getItem());
2832      result.append (modNTL (buf, j.getItem()));
2833    }
2834    return result;
2835  }
2836  Variable y= M.getLast().mvar();
2837  CFList bufFactors= factors;
2838  for (CFListIterator i= bufFactors; i.hasItem(); i++)
2839    i.getItem()= mod (i.getItem(), y);
2840  CFList bufProducts= products;
2841  for (CFListIterator i= bufProducts; i.hasItem(); i++)
2842    i.getItem()= mod (i.getItem(), y);
2843  CFList buf= M;
2844  buf.removeLast();
2845  CanonicalForm bufE= mod (E, y);
2846  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2847                                      bufE, bad);
2848
2849  if (bad)
2850    return CFList();
2851
2852  CanonicalForm e= E;
2853  CFListIterator j= products;
2854  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
2855    e -= i.getItem()*j.getItem();
2856
2857  CFList result= recDiophantine;
2858  int d= degree (M.getLast());
2859  CanonicalForm coeffE;
2860  for (int i= 1; i < d; i++)
2861  {
2862    if (degree (e, y) > 0)
2863      coeffE= e[i];
2864    else
2865      coeffE= 0;
2866    if (!coeffE.isZero())
2867    {
2868      CFListIterator k= result;
2869      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2870                                   coeffE, bad);
2871      if (bad)
2872        return CFList();
2873      CFListIterator l= products;
2874      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
2875      {
2876        k.getItem() += j.getItem()*power (y, i);
2877        e -= j.getItem()*power (y, i)*l.getItem();
2878      }
2879    }
2880    if (e.isZero())
2881      break;
2882  }
2883  if (!e.isZero())
2884  {
2885    bad= true;
2886    return CFList();
2887  }
2888  return result;
2889}
2890
2891static inline
2892void
2893henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
2894            const CFList& diophant, CFMatrix& M, CFArray& Pi,
2895            const CFList& products, int j, const CFList& MOD, bool& noOneToOne)
2896{
2897  CanonicalForm E;
2898  CanonicalForm xToJ= power (F.mvar(), j);
2899  Variable x= F.mvar();
2900
2901  // compute the error
2902#ifdef DEBUGOUTPUT
2903    CanonicalForm test= 1;
2904    for (int i= 0; i < factors.length(); i++)
2905    {
2906      if (i == 0)
2907        test *= mod (bufFactors [i], power (x, j));
2908      else
2909        test *= bufFactors[i];
2910    }
2911    test= mod (test, power (x, j));
2912    test= mod (test, MOD);
2913    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
2914    DEBOUTLN (cerr, "test= " << test2);
2915#endif
2916
2917  if (degree (Pi [factors.length() - 2], x) > 0)
2918    E= F[j] - Pi [factors.length() - 2] [j];
2919  else
2920    E= F[j];
2921
2922  CFArray buf= CFArray (diophant.length());
2923
2924  // actual lifting
2925  CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
2926                                    noOneToOne);
2927
2928  if (noOneToOne)
2929    return;
2930
2931  int k= 0;
2932  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
2933  {
2934    buf[k]= i.getItem();
2935    bufFactors[k] += xToJ*i.getItem();
2936  }
2937
2938  // update Pi [0]
2939  int degBuf0= degree (bufFactors[0], x);
2940  int degBuf1= degree (bufFactors[1], x);
2941  if (degBuf0 > 0 && degBuf1 > 0)
2942  {
2943    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
2944    if (j + 2 <= M.rows())
2945      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
2946  }
2947  CanonicalForm uIZeroJ;
2948
2949  if (degBuf0 > 0 && degBuf1 > 0)
2950    uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
2951             mulMod (bufFactors[1] [0], buf[0], MOD);
2952  else if (degBuf0 > 0)
2953    uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
2954  else if (degBuf1 > 0)
2955    uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
2956  else
2957    uIZeroJ= 0;
2958  Pi [0] += xToJ*uIZeroJ;
2959
2960  CFArray tmp= CFArray (factors.length() - 1);
2961  for (k= 0; k < factors.length() - 1; k++)
2962    tmp[k]= 0;
2963  CFIterator one, two;
2964  one= bufFactors [0];
2965  two= bufFactors [1];
2966  if (degBuf0 > 0 && degBuf1 > 0)
2967  {
2968    while (one.hasTerms() && one.exp() > j) one++;
2969    while (two.hasTerms() && two.exp() > j) two++;
2970    for (k= 1; k <= (int) ceil (j/2.0); k++)
2971    {
2972      if (k != j - k + 1)
2973      {
2974        if ((one.hasTerms() && one.exp() == j - k + 1) &&
2975            (two.hasTerms() && two.exp() == j - k + 1))
2976        {
2977          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2978                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
2979                    M (j - k + 2, 1);
2980          one++;
2981          two++;
2982        }
2983        else if (one.hasTerms() && one.exp() == j - k + 1)
2984        {
2985          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2986                    bufFactors[1] [k], MOD) - M (k + 1, 1);
2987          one++;
2988        }
2989        else if (two.hasTerms() && two.exp() == j - k + 1)
2990        {
2991          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
2992                    two.coeff()), MOD) - M (k + 1, 1);
2993          two++;
2994        }
2995      }
2996      else
2997      {
2998        tmp[0] += M (k + 1, 1);
2999      }
3000    }
3001  }
3002
3003  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
3004  {
3005    if (j + 2 <= M.rows())
3006      tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
3007                         (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
3008                         - M(1,1) - M (j + 2,1);
3009  }
3010  else if (degBuf0 >= j + 1)
3011  {
3012    if (degBuf1 > 0)
3013      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
3014    else
3015      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
3016  }
3017  else if (degBuf1 >= j + 1)
3018  {
3019    if (degBuf0 > 0)
3020      tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
3021    else
3022      tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
3023  }
3024  Pi [0] += tmp[0]*xToJ*F.mvar();
3025
3026  // update Pi [l]
3027  int degPi, degBuf;
3028  for (int l= 1; l < factors.length() - 1; l++)
3029  {
3030    degPi= degree (Pi [l - 1], x);
3031    degBuf= degree (bufFactors[l + 1], x);
3032    if (degPi > 0 && degBuf > 0)
3033    {
3034      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
3035      if (j + 2 <= M.rows())
3036        M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
3037                                  MOD);
3038    }
3039
3040    if (degPi > 0 && degBuf > 0)
3041      uIZeroJ= mulMod (Pi[l -1] [0], buf[l + 1], MOD) +
3042               mulMod (uIZeroJ, bufFactors[l+1] [0], MOD);
3043    else if (degPi > 0)
3044      uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD);
3045    else if (degBuf > 0)
3046      uIZeroJ= mulMod (Pi[l - 1], buf[1], MOD);
3047    else
3048      uIZeroJ= 0;
3049
3050    Pi [l] += xToJ*uIZeroJ;
3051
3052    one= bufFactors [l + 1];
3053    two= Pi [l - 1];
3054    if (degBuf > 0 && degPi > 0)
3055    {
3056      while (one.hasTerms() && one.exp() > j) one++;
3057      while (two.hasTerms() && two.exp() > j) two++;
3058      for (k= 1; k <= (int) ceil (j/2.0); k++)
3059      {
3060        if (k != j - k + 1)
3061        {
3062          if ((one.hasTerms() && one.exp() == j - k + 1) &&
3063              (two.hasTerms() && two.exp() == j - k + 1))
3064          {
3065            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
3066                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
3067                      M (j - k + 2, l + 1);
3068            one++;
3069            two++;
3070          }
3071          else if (one.hasTerms() && one.exp() == j - k + 1)
3072          {
3073            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
3074                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
3075            one++;
3076          }
3077          else if (two.hasTerms() && two.exp() == j - k + 1)
3078          {
3079            tmp[l] += mulMod (bufFactors[l + 1] [k],
3080                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
3081            two++;
3082           }
3083        }
3084        else
3085          tmp[l] += M (k + 1, l + 1);
3086      }
3087    }
3088
3089    if (degPi >= j + 1 && degBuf >= j + 1)
3090    {
3091      if (j + 2 <= M.rows())
3092        tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
3093                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
3094                            , MOD) - M(1,l+1) - M (j + 2,l+1);
3095    }
3096    else if (degPi >= j + 1)
3097    {
3098      if (degBuf > 0)
3099        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
3100      else
3101        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
3102    }
3103    else if (degBuf >= j + 1)
3104    {
3105      if (degPi > 0)
3106        tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
3107      else
3108        tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
3109    }
3110
3111    Pi[l] += tmp[l]*xToJ*F.mvar();
3112  }
3113  return;
3114}
3115
3116// wrt. Variable (1)
3117CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
3118{
3119  if (degree (F, 1) <= 0)
3120   return c;
3121  else
3122  {
3123    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
3124    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
3125              - LC (result))*power (result.mvar(), degree (result));
3126    return swapvar (result, Variable (F.level() + 1), Variable (1));
3127  }
3128}
3129
3130CFList
3131henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
3132              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
3133              const CFList& LCs2, bool& bad)
3134{
3135  CFList buf= factors;
3136  int k= 0;
3137  int liftBoundBivar= l[k];
3138  CFList bufbuf= factors;
3139  Variable v= Variable (2);
3140
3141  CFList MOD;
3142  MOD.append (power (Variable (2), liftBoundBivar));
3143  CFArray bufFactors= CFArray (factors.length());
3144  k= 0;
3145  CFListIterator j= eval;
3146  j++;
3147  CFListIterator iter1= LCs1;
3148  CFListIterator iter2= LCs2;
3149  iter1++;
3150  iter2++;
3151  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
3152  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
3153
3154  CFListIterator i= buf;
3155  i++;
3156  Variable y= j.getItem().mvar();
3157  if (y.level() != 3)
3158    y= Variable (3);
3159
3160  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
3161  M (1, 1)= Pi[0];
3162  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
3163    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
3164                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
3165  else if (degree (bufFactors[0], y) > 0)
3166    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
3167  else if (degree (bufFactors[1], y) > 0)
3168    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
3169
3170  CFList products;
3171  for (int i= 0; i < bufFactors.size(); i++)
3172  {
3173    if (degree (bufFactors[i], y) > 0)
3174      products.append (eval.getFirst()/bufFactors[i] [0]);
3175    else
3176      products.append (eval.getFirst()/bufFactors[i]);
3177  }
3178
3179  for (int d= 1; d < l[1]; d++)
3180  {
3181    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
3182    if (bad)
3183      return CFList();
3184  }
3185  CFList result;
3186  for (k= 0; k < factors.length(); k++)
3187    result.append (bufFactors[k]);
3188  return result;
3189}
3190
3191
3192CFList
3193henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
3194            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
3195            lNew, const CFList& LCs1, const CFList& LCs2, bool& bad)
3196{
3197  int k= 0;
3198  CFArray bufFactors= CFArray (factors.length());
3199  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
3200  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
3201  CFList buf= factors;
3202  Variable y= F.getLast().mvar();
3203  Variable x= F.getFirst().mvar();
3204  CanonicalForm xToLOld= power (x, lOld);
3205  Pi [0]= mod (Pi[0], xToLOld);
3206  M (1, 1)= Pi [0];
3207
3208  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
3209    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
3210                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
3211  else if (degree (bufFactors[0], y) > 0)
3212    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
3213  else if (degree (bufFactors[1], y) > 0)
3214    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
3215
3216  CFList products;
3217  CanonicalForm quot;
3218  for (int i= 0; i < bufFactors.size(); i++)
3219  {
3220    if (degree (bufFactors[i], y) > 0)
3221    {
3222      if (!fdivides (bufFactors[i] [0], F.getFirst(), quot))
3223      {
3224        bad= true;
3225        return CFList();
3226      }
3227      products.append (quot);
3228    }
3229    else
3230    {
3231      if (!fdivides (bufFactors[i], F.getFirst(), quot))
3232      {
3233        bad= true;
3234        return CFList();
3235      }
3236      products.append (quot);
3237    }
3238  }
3239
3240  for (int d= 1; d < lNew; d++)
3241  {
3242    henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
3243    if (bad)
3244      return CFList();
3245  }
3246
3247  CFList result;
3248  for (k= 0; k < factors.length(); k++)
3249    result.append (bufFactors[k]);
3250  return result;
3251}
3252
3253CFList
3254henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
3255             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
3256             const CFArray& Pi, const CFList& diophant, bool& bad)
3257{
3258  CFList bufDiophant= diophant;
3259  CFList buf= factors;
3260  if (sort)
3261    sortList (buf, Variable (1));
3262  CFArray bufPi= Pi;
3263  CFMatrix M= CFMatrix (l[1], factors.length());
3264  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2,
3265                               bad);
3266  if (bad)
3267    return CFList();
3268
3269  if (eval.length() == 2)
3270    return result;
3271  CFList MOD;
3272  for (int i= 0; i < 2; i++)
3273    MOD.append (power (Variable (i + 2), l[i]));
3274  CFListIterator j= eval;
3275  j++;
3276  CFList bufEval;
3277  bufEval.append (j.getItem());
3278  j++;
3279  CFListIterator jj= LCs1;
3280  CFListIterator jjj= LCs2;
3281  CFList bufLCs1, bufLCs2;
3282  jj++, jjj++;
3283  bufLCs1.append (jj.getItem());
3284  bufLCs2.append (jjj.getItem());
3285  jj++, jjj++;
3286
3287  for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
3288  {
3289    bufEval.append (j.getItem());
3290    bufLCs1.append (jj.getItem());
3291    bufLCs2.append (jjj.getItem());
3292    M= CFMatrix (l[i], factors.length());
3293    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
3294                         l[i], bufLCs1, bufLCs2, bad);
3295    if (bad)
3296      return CFList();
3297    MOD.append (power (Variable (i + 2), l[i]));
3298    bufEval.removeFirst();
3299    bufLCs1.removeFirst();
3300    bufLCs2.removeFirst();
3301  }
3302  return result;
3303}
3304
3305CFList
3306nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
3307                      CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
3308                      int bivarLiftBound, bool& bad)
3309{
3310  CFList bufFactors2= factors;
3311
3312  Variable y= Variable (2);
3313  for (CFListIterator i= bufFactors2; i.hasItem(); i++)
3314    i.getItem()= mod (i.getItem(), y);
3315
3316  CanonicalForm bufF= F;
3317  bufF= mod (bufF, y);
3318  bufF= mod (bufF, Variable (3));
3319
3320  diophant= diophantine (bufF, bufFactors2);
3321
3322  CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
3323
3324  Pi= CFArray (bufFactors2.length() - 1);
3325
3326  CFArray bufFactors= CFArray (bufFactors2.length());
3327  CFListIterator j= LCs;
3328  int i= 0;
3329  for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
3330    bufFactors[i]= replaceLC (k.getItem(), j.getItem());
3331
3332  //initialise Pi
3333  Variable v= Variable (3);
3334  CanonicalForm yToL= power (y, bivarLiftBound);
3335  if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
3336  {
3337    M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
3338    Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
3339                       mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
3340  }
3341  else if (degree (bufFactors[0], v) > 0)
3342  {
3343    M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
3344    Pi [0]=  M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
3345  }
3346  else if (degree (bufFactors[1], v) > 0)
3347  {
3348    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
3349    Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
3350  }
3351  else
3352  {
3353    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
3354    Pi [0]= M (1,1);
3355  }
3356
3357  for (i= 1; i < Pi.size(); i++)
3358  {
3359    if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
3360    {
3361      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
3362      Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
3363                       mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
3364    }
3365    else if (degree (Pi[i-1], v) > 0)
3366    {
3367      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
3368      Pi [i]=  M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
3369    }
3370    else if (degree (bufFactors[i+1], v) > 0)
3371    {
3372      M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
3373      Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
3374    }
3375    else
3376    {
3377      M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
3378      Pi [i]= M (1,i+1);
3379    }
3380  }
3381
3382  CFList products;
3383  bufF= mod (F, Variable (3));
3384  for (CFListIterator k= factors; k.hasItem(); k++)
3385    products.append (bufF/k.getItem());
3386
3387  CFList MOD= CFList (power (v, liftBound));
3388  MOD.insert (yToL);
3389  for (int d= 1; d < liftBound; d++)
3390  {
3391    henselStep2 (F, factors, bufFactors, diophant, M, Pi, products, d, MOD, bad);
3392    if (bad)
3393      return CFList();
3394  }
3395
3396  CFList result;
3397  for (i= 0; i < factors.length(); i++)
3398    result.append (bufFactors[i]);
3399  return result;
3400}
3401
3402CFList
3403nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
3404                    CFList& diophant, CFArray& Pi, CFMatrix& M, const int lOld,
3405                    int& lNew, const CFList& MOD, bool& noOneToOne
3406                   )
3407{
3408
3409  int k= 0;
3410  CFArray bufFactors= CFArray (factors.length());
3411  CFListIterator j= LCs;
3412  for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
3413    bufFactors [k]= replaceLC (i.getItem(), j.getItem());
3414
3415  Variable y= F.getLast().mvar();
3416  Variable x= F.getFirst().mvar();
3417  CanonicalForm xToLOld= power (x, lOld);
3418
3419  Pi [0]= mod (Pi[0], xToLOld);
3420  M (1, 1)= Pi [0];
3421
3422  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
3423    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
3424                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
3425  else if (degree (bufFactors[0], y) > 0)
3426    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
3427  else if (degree (bufFactors[1], y) > 0)
3428    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
3429
3430  for (int i= 1; i < Pi.size(); i++)
3431  {
3432    Pi [i]= mod (Pi [i], xToLOld);
3433    M (1, i + 1)= Pi [i];
3434
3435    if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
3436      Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
3437                 mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
3438    else if (degree (Pi[i-1], y) > 0)
3439      Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
3440    else if (degree (bufFactors[i+1], y) > 0)
3441      Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
3442  }
3443
3444  CFList products;
3445  CanonicalForm quot, bufF= F.getFirst();
3446
3447  for (int i= 0; i < bufFactors.size(); i++)
3448  {
3449    if (degree (bufFactors[i], y) > 0)
3450    {
3451      if (!fdivides (bufFactors[i] [0], bufF, quot))
3452      {
3453        noOneToOne= true;
3454        return factors;
3455      }
3456      products.append (quot);
3457    }
3458    else
3459    {
3460      if (!fdivides (bufFactors[i], bufF, quot))
3461      {
3462        noOneToOne= true;
3463        return factors;
3464      }
3465      products.append (quot);
3466    }
3467  }
3468
3469  for (int d= 1; d < lNew; d++)
3470  {
3471    henselStep2 (F.getLast(), factors, bufFactors, diophant, M, Pi, products, d,
3472                 MOD, noOneToOne);
3473    if (noOneToOne)
3474      return CFList();
3475  }
3476
3477  CFList result;
3478  for (k= 0; k < factors.length(); k++)
3479    result.append (bufFactors[k]);
3480  return result;
3481}
3482
3483CFList
3484nonMonicHenselLift (const CFList& eval, const CFList& factors,
3485                    CFList* const& LCs, CFList& diophant, CFArray& Pi,
3486                    int* liftBound, int length, bool& noOneToOne
3487                   )
3488{
3489  CFList bufDiophant= diophant;
3490  CFList buf= factors;
3491  CFArray bufPi= Pi;
3492  CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
3493  int k= 0;
3494
3495  CFList result=
3496  nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
3497                        liftBound[1], liftBound[0], noOneToOne);
3498
3499  if (noOneToOne)
3500    return CFList();
3501
3502  if (eval.length() == 1)
3503    return result;
3504
3505  k++;
3506  CFList MOD;
3507  for (int i= 0; i < 2; i++)
3508    MOD.append (power (Variable (i + 2), liftBound[i]));
3509
3510  CFListIterator j= eval;
3511  CFList bufEval;
3512  bufEval.append (j.getItem());
3513  j++;
3514
3515  for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
3516  {
3517    bufEval.append (j.getItem());
3518    M= CFMatrix (liftBound[i], factors.length() - 1);
3519    result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
3520                                liftBound[i-1], liftBound[i], MOD, noOneToOne);
3521    if (noOneToOne)
3522      return result;
3523    MOD.append (power (Variable (i + 2), liftBound[i]));
3524    bufEval.removeFirst();
3525  }
3526
3527  return result;
3528}
3529
3530#endif
3531/* HAVE_NTL */
3532
Note: See TracBrowser for help on using the repository browser.