source: git/factory/facHensel.cc @ 9b2c8a

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