source: git/factory/facHensel.cc @ 96f9fdf

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