source: git/factory/facHensel.cc @ 935e83

spielwiese
Last change on this file since 935e83 was 935e83, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: deleted some commented out debug stuff and formatted code
  • Property mode set to 100644
File size: 117.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 "config.h"
23
24#include "cf_assert.h"
25#include "debug.h"
26#include "timing.h"
27
28#include "facHensel.h"
29#include "cf_util.h"
30#include "fac_util.h"
31#include "cf_algorithm.h"
32#include "cf_primes.h"
33
34#ifdef HAVE_NTL
35#include <NTL/lzz_pEX.h>
36#include "NTLconvert.h"
37
38#ifdef HAVE_FLINT
39#include "FLINTconvert.h"
40#endif
41
42#ifdef HAVE_FLINT
43void kronSub (fmpz_poly_t result, const CanonicalForm& A, int d)
44{
45  int degAy= degree (A);
46  fmpz_poly_init2 (result, d*(degAy + 1));
47  _fmpz_poly_set_length (result, d*(degAy + 1));
48  CFIterator j;
49  for (CFIterator i= A; i.hasTerms(); i++)
50  {
51    if (i.coeff().inBaseDomain())
52      convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d), i.coeff());
53    else
54      for (j= i.coeff(); j.hasTerms(); j++)
55        convertCF2Fmpz (fmpz_poly_get_coeff_ptr (result, i.exp()*d+j.exp()),
56                        j.coeff());
57  }
58  _fmpz_poly_normalise(result);
59}
60
61
62CanonicalForm
63reverseSubst (const fmpz_poly_t F, int d, const Variable& alpha,
64              const CanonicalForm& den)
65{
66  Variable x= Variable (1);
67
68  CanonicalForm result= 0;
69  int i= 0;
70  int degf= fmpz_poly_degree (F);
71  int k= 0;
72  int degfSubK;
73  int repLength, j;
74  CanonicalForm coeff;
75  fmpz* tmp;
76  while (degf >= k)
77  {
78    coeff= 0;
79    degfSubK= degf - k;
80    if (degfSubK >= d)
81      repLength= d;
82    else
83      repLength= degfSubK + 1;
84
85    for (j= 0; j < repLength; j++)
86    {
87      tmp= fmpz_poly_get_coeff_ptr (F, j+k);
88      if (!fmpz_is_zero (tmp))
89      {
90        CanonicalForm ff= convertFmpz2CF (tmp)/den;
91        coeff += ff*power (alpha, j);
92      }
93    }
94    result += coeff*power (x, i);
95    i++;
96    k= d*i;
97  }
98  return result;
99}
100
101CanonicalForm
102mulFLINTQa (const CanonicalForm& F, const CanonicalForm& G,
103            const Variable& alpha)
104{
105  CanonicalForm A= F;
106  CanonicalForm B= G;
107
108  CanonicalForm denA= bCommonDen (A);
109  CanonicalForm denB= bCommonDen (B);
110
111  A *= denA;
112  B *= denB;
113  int degAa= degree (A, alpha);
114  int degBa= degree (B, alpha);
115  int d= degAa + 1 + degBa;
116
117  fmpz_poly_t FLINTA,FLINTB;
118  fmpz_poly_init (FLINTA);
119  fmpz_poly_init (FLINTB);
120  kronSub (FLINTA, A, d);
121  kronSub (FLINTB, B, d);
122
123  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
124
125  denA *= denB;
126  A= reverseSubst (FLINTA, d, alpha, denA);
127
128  fmpz_poly_clear (FLINTA);
129  fmpz_poly_clear (FLINTB);
130  return A;
131}
132
133CanonicalForm
134mulFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
135{
136  CanonicalForm A= F;
137  CanonicalForm B= G;
138
139  CanonicalForm denA= bCommonDen (A);
140  CanonicalForm denB= bCommonDen (B);
141
142  A *= denA;
143  B *= denB;
144  fmpz_poly_t FLINTA,FLINTB;
145  convertFacCF2Fmpz_poly_t (FLINTA, A);
146  convertFacCF2Fmpz_poly_t (FLINTB, B);
147  fmpz_poly_mul (FLINTA, FLINTA, FLINTB);
148  denA *= denB;
149  A= convertFmpz_poly_t2FacCF (FLINTA, F.mvar());
150  A /= denA;
151  fmpz_poly_clear (FLINTA);
152  fmpz_poly_clear (FLINTB);
153
154  return A;
155}
156
157/*CanonicalForm
158mulFLINTQ2 (const CanonicalForm& F, const CanonicalForm& G)
159{
160  CanonicalForm A= F;
161  CanonicalForm B= G;
162
163  fmpq_poly_t FLINTA,FLINTB;
164  convertFacCF2Fmpq_poly_t (FLINTA, A);
165  convertFacCF2Fmpq_poly_t (FLINTB, B);
166
167  fmpq_poly_mul (FLINTA, FLINTA, FLINTB);
168  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
169  fmpq_poly_clear (FLINTA);
170  fmpq_poly_clear (FLINTB);
171  return A;
172}*/
173
174CanonicalForm
175divFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
176{
177  CanonicalForm A= F;
178  CanonicalForm B= G;
179
180  fmpq_poly_t FLINTA,FLINTB;
181  convertFacCF2Fmpq_poly_t (FLINTA, A);
182  convertFacCF2Fmpq_poly_t (FLINTB, B);
183
184  fmpq_poly_div (FLINTA, FLINTA, FLINTB);
185  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
186
187  fmpq_poly_clear (FLINTA);
188  fmpq_poly_clear (FLINTB);
189  return A;
190}
191
192CanonicalForm
193modFLINTQ (const CanonicalForm& F, const CanonicalForm& G)
194{
195  CanonicalForm A= F;
196  CanonicalForm B= G;
197
198  fmpq_poly_t FLINTA,FLINTB;
199  convertFacCF2Fmpq_poly_t (FLINTA, A);
200  convertFacCF2Fmpq_poly_t (FLINTB, B);
201
202  fmpq_poly_rem (FLINTA, FLINTA, FLINTB);
203  A= convertFmpq_poly_t2FacCF (FLINTA, F.mvar());
204
205  fmpq_poly_clear (FLINTA);
206  fmpq_poly_clear (FLINTB);
207  return A;
208}
209
210CanonicalForm
211mulFLINTQaTrunc (const CanonicalForm& F, const CanonicalForm& G,
212                 const Variable& alpha, int m)
213{
214  CanonicalForm A= F;
215  CanonicalForm B= G;
216
217  CanonicalForm denA= bCommonDen (A);
218  CanonicalForm denB= bCommonDen (B);
219
220  A *= denA;
221  B *= denB;
222
223  int degAa= degree (A, alpha);
224  int degBa= degree (B, alpha);
225  int d= degAa + 1 + degBa;
226
227  fmpz_poly_t FLINTA,FLINTB;
228  fmpz_poly_init (FLINTA);
229  fmpz_poly_init (FLINTB);
230  kronSub (FLINTA, A, d);
231  kronSub (FLINTB, B, d);
232
233  int k= d*m;
234  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, k);
235
236  denA *= denB;
237  A= reverseSubst (FLINTA, d, alpha, denA);
238  fmpz_poly_clear (FLINTA);
239  fmpz_poly_clear (FLINTB);
240  return A;
241}
242
243#endif
244
245static
246CFList productsNTL (const CFList& factors, const CanonicalForm& M)
247{
248  zz_p::init (getCharacteristic());
249  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
250  zz_pE::init (NTLMipo);
251  zz_pEX prod;
252  vec_zz_pEX v;
253  v.SetLength (factors.length());
254  int j= 0;
255  for (CFListIterator i= factors; i.hasItem(); i++, j++)
256  {
257    if (i.getItem().inCoeffDomain())
258      v[j]= to_zz_pEX (to_zz_pE (convertFacCF2NTLzzpX (i.getItem())));
259    else
260      v[j]= convertFacCF2NTLzz_pEX (i.getItem(), NTLMipo);
261  }
262  CFList result;
263  Variable x= Variable (1);
264  for (int j= 0; j < factors.length(); j++)
265  {
266    int k= 0;
267    set(prod);
268    for (int i= 0; i < factors.length(); i++, k++)
269    {
270      if (k == j)
271        continue;
272      prod *= v[i];
273    }
274    result.append (convertNTLzz_pEX2CF (prod, x, M.mvar()));
275  }
276  return result;
277}
278
279static
280void tryDiophantine (CFList& result, const CanonicalForm& F,
281                     const CFList& factors, const CanonicalForm& M, bool& fail)
282{
283  ASSERT (M.isUnivariate(), "expected univariate poly");
284
285  CFList bufFactors= factors;
286  bufFactors.removeFirst();
287  bufFactors.insert (factors.getFirst () (0,2));
288  CanonicalForm inv, leadingCoeff= Lc (F);
289  CFListIterator i= bufFactors;
290  if (bufFactors.getFirst().inCoeffDomain())
291  {
292    if (i.hasItem())
293      i++;
294  }
295  for (; i.hasItem(); i++)
296  {
297    tryInvert (Lc (i.getItem()), M, inv ,fail);
298    if (fail)
299      return;
300    i.getItem()= reduce (i.getItem()*inv, M);
301  }
302  bufFactors= productsNTL (bufFactors, M);
303
304  CanonicalForm buf1, buf2, buf3, S, T;
305  i= bufFactors;
306  if (i.hasItem())
307    i++;
308  buf1= bufFactors.getFirst();
309  buf2= i.getItem();
310  tryExtgcd (buf1, buf2, M, buf3, S, T, fail);
311  if (fail)
312    return;
313  result.append (S);
314  result.append (T);
315  if (i.hasItem())
316    i++;
317  for (; i.hasItem(); i++)
318  {
319    buf1= i.getItem();
320    tryExtgcd (buf3, buf1, M, buf3, S, T, fail);
321    if (fail)
322      return;
323    CFListIterator k= factors;
324    for (CFListIterator j= result; j.hasItem(); j++, k++)
325    {
326      j.getItem() *= S;
327      j.getItem()= mod (j.getItem(), k.getItem());
328      j.getItem()= reduce (j.getItem(), M);
329    }
330    result.append (T);
331  }
332}
333
334static inline
335CFList mapinto (const CFList& L)
336{
337  CFList result;
338  for (CFListIterator i= L; i.hasItem(); i++)
339    result.append (mapinto (i.getItem()));
340  return result;
341}
342
343static inline
344int mod (const CFList& L, const CanonicalForm& p)
345{
346  for (CFListIterator i= L; i.hasItem(); i++)
347  {
348    if (mod (i.getItem(), p) == 0)
349      return 0;
350  }
351  return 1;
352}
353
354
355static inline void
356chineseRemainder (const CFList & x1, const CanonicalForm & q1,
357                  const CFList & x2, const CanonicalForm & q2,
358                  CFList & xnew, CanonicalForm & qnew)
359{
360  ASSERT (x1.length() == x2.length(), "expected lists of equal length");
361  CanonicalForm tmp1, tmp2;
362  CFListIterator j= x2;
363  for (CFListIterator i= x1; i.hasItem() && j.hasItem(); i++, j++)
364  {
365    chineseRemainder (i.getItem(), q1, j.getItem(), q2, tmp1, tmp2);
366    xnew.append (tmp1);
367  }
368  qnew= tmp2;
369}
370
371static inline
372CFList Farey (const CFList& L, const CanonicalForm& q)
373{
374  CFList result;
375  for (CFListIterator i= L; i.hasItem(); i++)
376    result.append (Farey (i.getItem(), q));
377  return result;
378}
379
380static inline
381CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
382{
383  CFList result;
384  for (CFListIterator i= L; i.hasItem(); i++)
385    result.append (replacevar (i.getItem(), a, b));
386  return result;
387}
388
389CFList
390modularDiophant (const CanonicalForm& f, const CFList& factors,
391                 const CanonicalForm& M)
392{
393  bool save_rat=!isOn (SW_RATIONAL);
394  On (SW_RATIONAL);
395  CanonicalForm F= f*bCommonDen (f);
396  CFList products= factors;
397  for (CFListIterator i= products; i.hasItem(); i++)
398  {
399    if (products.getFirst().level() == 1)
400      i.getItem() /= Lc (i.getItem());
401    i.getItem() *= bCommonDen (i.getItem());
402  }
403  if (products.getFirst().level() == 1)
404    products.insert (Lc (F));
405  CanonicalForm bound= maxNorm (F);
406  CFList leadingCoeffs;
407  leadingCoeffs.append (lc (F));
408  CanonicalForm dummy;
409  for (CFListIterator i= products; i.hasItem(); i++)
410  {
411    leadingCoeffs.append (lc (i.getItem()));
412    dummy= maxNorm (i.getItem());
413    bound= (dummy > bound) ? dummy : bound;
414  }
415  bound *= maxNorm (Lc (F))*maxNorm (Lc(F))*bound;
416  bound *= bound*bound;
417  bound= power (bound, degree (M));
418  bound *= power (CanonicalForm (2),degree (f));
419  CanonicalForm bufBound= bound;
420  int i = cf_getNumBigPrimes() - 1;
421  int p;
422  CFList resultModP, result, newResult;
423  CanonicalForm q (0), newQ;
424  bool fail= false;
425  Variable a= M.mvar();
426  Variable b= Variable (2);
427  setReduce (M.mvar(), false);
428  CanonicalForm mipo= bCommonDen (M)*M;
429  Off (SW_RATIONAL);
430  CanonicalForm modMipo;
431  leadingCoeffs.append (lc (mipo));
432  CFList tmp1, tmp2;
433  bool equal= false;
434  int count= 0;
435  do
436  {
437    p = cf_getBigPrime( i );
438    i--;
439    while ( i >= 0 && mod( leadingCoeffs, p ) == 0)
440    {
441      p = cf_getBigPrime( i );
442      i--;
443    }
444
445    ASSERT (i >= 0, "ran out of primes"); //sic
446
447    setCharacteristic (p);
448    modMipo= mapinto (mipo);
449    modMipo /= lc (modMipo);
450    resultModP= CFList();
451    tryDiophantine (resultModP, mapinto (F), mapinto (products), modMipo, fail);
452    setCharacteristic (0);
453    if (fail)
454    {
455      fail= false;
456      continue;
457    }
458
459    if ( q.isZero() )
460    {
461      result= replacevar (mapinto(resultModP), a, b);
462      q= p;
463    }
464    else
465    {
466      result= replacevar (result, a, b);
467      newResult= CFList();
468      chineseRemainder( result, q, replacevar (mapinto (resultModP), a, b),
469                        p, newResult, newQ );
470      q= newQ;
471      result= newResult;
472      if (newQ > bound)
473      {
474        count++;
475        tmp1= replacevar (Farey (result, q), b, a);
476        if (tmp2.isEmpty())
477          tmp2= tmp1;
478        else
479        {
480          equal= true;
481          CFListIterator k= tmp1;
482          for (CFListIterator j= tmp2; j.hasItem(); j++, k++)
483          {
484            if (j.getItem() != k.getItem())
485              equal= false;
486          }
487          if (!equal)
488            tmp2= tmp1;
489        }
490        if (count > 2)
491        {
492          bound *= bufBound;
493          equal= false;
494          count= 0;
495        }
496      }
497      if (newQ > bound && equal)
498      {
499        On( SW_RATIONAL );
500        CFList bufResult= result;
501        result= tmp2;
502        setReduce (M.mvar(), true);
503        if (factors.getFirst().level() == 1)
504        {
505          result.removeFirst();
506          CFListIterator j= factors;
507          CanonicalForm denf= bCommonDen (f);
508          for (CFListIterator i= result; i.hasItem(); i++, j++)
509            i.getItem() *= Lc (j.getItem())*denf;
510        }
511        if (factors.getFirst().level() != 1 && 
512            !bCommonDen (factors.getFirst()).isOne())
513        {
514          CanonicalForm denFirst= bCommonDen (factors.getFirst());
515          for (CFListIterator i= result; i.hasItem(); i++)
516            i.getItem() *= denFirst;
517        }
518
519        CanonicalForm test= 0;
520        CFListIterator jj= factors;
521        for (CFListIterator ii= result; ii.hasItem(); ii++, jj++)
522          test += ii.getItem()*(f/jj.getItem());
523        if (!test.isOne())
524        {
525          bound *= bufBound;
526          equal= false;
527          count= 0;
528          setReduce (M.mvar(), false);
529          result= bufResult;
530          Off (SW_RATIONAL);
531        }
532        else
533          break;
534      }
535    }
536  } while (1);
537  if (save_rat) Off(SW_RATIONAL);
538  return result;
539}
540
541CanonicalForm
542mulNTL (const CanonicalForm& F, const CanonicalForm& G)
543{
544  if (F.inCoeffDomain() || G.inCoeffDomain() || getCharacteristic() == 0)
545  {
546    Variable alpha;
547#ifdef HAVE_FLINT
548    if ((!F.inCoeffDomain() && !G.inCoeffDomain()) &&
549        (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha)))
550    {
551      CanonicalForm result= mulFLINTQa (F, G, alpha);
552      return result;
553    }
554    else if (!F.inCoeffDomain() && !G.inCoeffDomain())
555      return mulFLINTQ (F, G);
556#endif
557    return F*G;
558  }
559  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
560  ASSERT (F.level() == G.level(), "expected polys of same level");
561  if (CFFactory::gettype() == GaloisFieldDomain)
562    return F*G;
563  zz_p::init (getCharacteristic());
564  Variable alpha;
565  CanonicalForm result;
566  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
567  {
568    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
569    zz_pE::init (NTLMipo);
570    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
571    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
572    mul (NTLF, NTLF, NTLG);
573    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
574  }
575  else
576  {
577#ifdef HAVE_FLINT
578    nmod_poly_t FLINTF, FLINTG;
579    convertFacCF2nmod_poly_t (FLINTF, F);
580    convertFacCF2nmod_poly_t (FLINTG, G);
581    nmod_poly_mul (FLINTF, FLINTF, FLINTG);
582    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
583    nmod_poly_clear (FLINTF);
584    nmod_poly_clear (FLINTG);
585#else
586    zz_pX NTLF= convertFacCF2NTLzzpX (F);
587    zz_pX NTLG= convertFacCF2NTLzzpX (G);
588    mul (NTLF, NTLF, NTLG);
589    result= convertNTLzzpX2CF(NTLF, F.mvar());
590#endif
591  }
592  return result;
593}
594
595CanonicalForm
596modNTL (const CanonicalForm& F, const CanonicalForm& G)
597{
598  if (F.inCoeffDomain() && G.isUnivariate())
599    return F;
600  else if (F.inCoeffDomain() && G.inCoeffDomain())
601    return mod (F, G);
602  else if (F.isUnivariate() && G.inCoeffDomain())
603    return mod (F,G);
604
605  if (getCharacteristic() == 0)
606  {
607#ifdef HAVE_FLINT
608    Variable alpha;
609    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
610      return modFLINTQ (F, G);
611    else
612      //TODO newtonDivrem
613     return mod (F, G);
614#else
615    return mod (F, G);
616#endif
617  }
618
619  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
620  ASSERT (F.level() == G.level(), "expected polys of same level");
621  if (CFFactory::gettype() == GaloisFieldDomain)
622    return mod (F, G);
623  zz_p::init (getCharacteristic());
624  Variable alpha;
625  CanonicalForm result;
626  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
627  {
628    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
629    zz_pE::init (NTLMipo);
630    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
631    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
632    rem (NTLF, NTLF, NTLG);
633    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
634  }
635  else
636  {
637#ifdef HAVE_FLINT
638    nmod_poly_t FLINTF, FLINTG;
639    convertFacCF2nmod_poly_t (FLINTF, F);
640    convertFacCF2nmod_poly_t (FLINTG, G);
641    nmod_poly_divrem (FLINTG, FLINTF, FLINTF, FLINTG);
642    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
643    nmod_poly_clear (FLINTF);
644    nmod_poly_clear (FLINTG);
645#else
646    zz_pX NTLF= convertFacCF2NTLzzpX (F);
647    zz_pX NTLG= convertFacCF2NTLzzpX (G);
648    rem (NTLF, NTLF, NTLG);
649    result= convertNTLzzpX2CF(NTLF, F.mvar());
650#endif
651  }
652  return result;
653}
654
655CanonicalForm
656divNTL (const CanonicalForm& F, const CanonicalForm& G)
657{
658  if (F.inCoeffDomain() && G.isUnivariate())
659    return F;
660  else if (F.inCoeffDomain() && G.inCoeffDomain())
661    return div (F, G);
662  else if (F.isUnivariate() && G.inCoeffDomain())
663    return div (F,G);
664
665  if (getCharacteristic() == 0)
666  {
667#ifdef HAVE_FLINT
668    Variable alpha;
669    if (!hasFirstAlgVar (F, alpha) && !hasFirstAlgVar (G, alpha))
670      return divFLINTQ (F,G);
671    else
672      //TODO newtonDivrem
673      return div (F, G);
674#else
675    return div (F, G);
676#endif
677  }
678
679  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
680  ASSERT (F.level() == G.level(), "expected polys of same level");
681  if (CFFactory::gettype() == GaloisFieldDomain)
682    return div (F, G);
683  zz_p::init (getCharacteristic());
684  Variable alpha;
685  CanonicalForm result;
686  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
687  {
688    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
689    zz_pE::init (NTLMipo);
690    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
691    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
692    div (NTLF, NTLF, NTLG);
693    result= convertNTLzz_pEX2CF(NTLF, F.mvar(), alpha);
694  }
695  else
696  {
697#ifdef HAVE_FLINT
698    nmod_poly_t FLINTF, FLINTG;
699    convertFacCF2nmod_poly_t (FLINTF, F);
700    convertFacCF2nmod_poly_t (FLINTG, G);
701    nmod_poly_div (FLINTF, FLINTF, FLINTG);
702    result= convertnmod_poly_t2FacCF (FLINTF, F.mvar());
703    nmod_poly_clear (FLINTF);
704    nmod_poly_clear (FLINTG);
705#else
706    zz_pX NTLF= convertFacCF2NTLzzpX (F);
707    zz_pX NTLG= convertFacCF2NTLzzpX (G);
708    div (NTLF, NTLF, NTLG);
709    result= convertNTLzzpX2CF(NTLF, F.mvar());
710#endif
711  }
712  return result;
713}
714
715/*
716void
717divremNTL (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
718           CanonicalForm& R)
719{
720  if (F.inCoeffDomain() && G.isUnivariate())
721  {
722    R= F;
723    Q= 0;
724  }
725  else if (F.inCoeffDomain() && G.inCoeffDomain())
726  {
727    divrem (F, G, Q, R);
728    return;
729  }
730  else if (F.isUnivariate() && G.inCoeffDomain())
731  {
732    divrem (F, G, Q, R);
733    return;
734  }
735
736  ASSERT (F.isUnivariate() && G.isUnivariate(), "expected univariate polys");
737  ASSERT (F.level() == G.level(), "expected polys of same level");
738  if (CFFactory::gettype() == GaloisFieldDomain)
739  {
740    divrem (F, G, Q, R);
741    return;
742  }
743  zz_p::init (getCharacteristic());
744  Variable alpha;
745  CanonicalForm result;
746  if (hasFirstAlgVar (F, alpha) || hasFirstAlgVar (G, alpha))
747  {
748    zz_pX NTLMipo= convertFacCF2NTLzzpX(getMipo (alpha));
749    zz_pE::init (NTLMipo);
750    zz_pEX NTLF= convertFacCF2NTLzz_pEX (F, NTLMipo);
751    zz_pEX NTLG= convertFacCF2NTLzz_pEX (G, NTLMipo);
752    zz_pEX NTLQ;
753    zz_pEX NTLR;
754    DivRem (NTLQ, NTLR, NTLF, NTLG);
755    Q= convertNTLzz_pEX2CF(NTLQ, F.mvar(), alpha);
756    R= convertNTLzz_pEX2CF(NTLR, F.mvar(), alpha);
757    return;
758  }
759  else
760  {
761    zz_pX NTLF= convertFacCF2NTLzzpX (F);
762    zz_pX NTLG= convertFacCF2NTLzzpX (G);
763    zz_pX NTLQ;
764    zz_pX NTLR;
765    DivRem (NTLQ, NTLR, NTLF, NTLG);
766    Q= convertNTLzzpX2CF(NTLQ, F.mvar());
767    R= convertNTLzzpX2CF(NTLR, F.mvar());
768    return;
769  }
770}*/
771
772CanonicalForm mod (const CanonicalForm& F, const CFList& M)
773{
774  CanonicalForm A= F;
775  for (CFListIterator i= M; i.hasItem(); i++)
776    A= mod (A, i.getItem());
777  return A;
778}
779
780#ifdef HAVE_FLINT
781void kronSubFp (nmod_poly_t result, const CanonicalForm& A, int d)
782{
783  int degAy= degree (A);
784  nmod_poly_init2 (result, getCharacteristic(), d*(degAy + 1));
785
786  nmod_poly_t buf;
787
788  int j, k, bufRepLength;
789  for (CFIterator i= A; i.hasTerms(); i++)
790  {
791    convertFacCF2nmod_poly_t (buf, i.coeff());
792
793    k= i.exp()*d;
794    bufRepLength= (int) nmod_poly_length (buf);
795    for (j= 0; j < bufRepLength; j++)
796      nmod_poly_set_coeff_ui (result, j + k, nmod_poly_get_coeff_ui (buf, j));
797    nmod_poly_clear (buf);
798  }
799  _nmod_poly_normalise (result);
800}
801
802/*void kronSubQ (fmpz_poly_t result, const CanonicalForm& A, int d)
803{
804  int degAy= degree (A);
805  fmpz_poly_init2 (result, d*(degAy + 1));
806  _fmpz_poly_set_length (result, d*(degAy+1));
807
808  CFIterator j;
809  for (CFIterator i= A; i.hasTerms(); i++)
810  {
811    if (i.coeff().inBas
812    convertFacCF2Fmpz_poly_t (buf, i.coeff());
813
814    int k= i.exp()*d;
815    int bufRepLength= (int) fmpz_poly_length (buf);
816    for (int j= 0; j < bufRepLength; j++)
817    {
818      fmpz_poly_get_coeff_fmpz (coeff, buf, j);
819      fmpz_poly_set_coeff_fmpz (result, j + k, coeff);
820    }
821    fmpz_poly_clear (buf);
822  }
823  fmpz_clear (coeff);
824  _fmpz_poly_normalise (result);
825}*/
826
827// A is a bivariate poly over Qa!!!!
828// d2= 2*deg_alpha + 1
829// d1= 2*deg_x*d2+1
830void kronSubQa (fmpq_poly_t result, const CanonicalForm& A, int d1, int d2)
831{
832  int degAy= degree (A);
833  fmpq_poly_init2 (result, d1*(degAy + 1));
834
835  fmpq_poly_t buf;
836  fmpq_t coeff;
837
838  int k, l, bufRepLength;
839  CFIterator j;
840  for (CFIterator i= A; i.hasTerms(); i++)
841  {
842    if (i.coeff().inCoeffDomain())
843    {
844      k= d1*i.exp();
845      convertFacCF2Fmpq_poly_t (buf, i.coeff());
846      bufRepLength= (int) fmpq_poly_length(buf);
847      for (l= 0; l < bufRepLength; l++)
848      {
849        fmpq_poly_get_coeff_fmpq (coeff, buf, l);
850        fmpq_poly_set_coeff_fmpq (result, l + k, coeff);
851      }
852      fmpq_poly_clear (buf);
853    }
854    else
855    {
856      for (j= i.coeff(); j.hasTerms(); j++)
857      {
858        k= d1*i.exp();
859        k += d2*j.exp();
860        convertFacCF2Fmpq_poly_t (buf, j.coeff());
861        bufRepLength= (int) fmpq_poly_length(buf);
862        for (l= 0; l < bufRepLength; l++)
863        {
864          fmpq_poly_get_coeff_fmpq (coeff, buf, l);
865          fmpq_poly_set_coeff_fmpq (result, k + l, coeff);
866        }
867        fmpq_poly_clear (buf);
868      }
869    }
870  }
871  fmpq_clear (coeff);
872  _fmpq_poly_normalise (result);
873}
874#endif
875
876zz_pX kronSubFp (const CanonicalForm& A, int d)
877{
878  int degAy= degree (A);
879  zz_pX result;
880  result.rep.SetLength (d*(degAy + 1));
881
882  zz_p *resultp;
883  resultp= result.rep.elts();
884  zz_pX buf;
885  zz_p *bufp;
886  int j, k, bufRepLength;
887
888  for (CFIterator i= A; i.hasTerms(); i++)
889  {
890    if (i.coeff().inCoeffDomain())
891      buf= convertFacCF2NTLzzpX (i.coeff());
892    else
893      buf= convertFacCF2NTLzzpX (i.coeff());
894
895    k= i.exp()*d;
896    bufp= buf.rep.elts();
897    bufRepLength= (int) buf.rep.length();
898    for (j= 0; j < bufRepLength; j++)
899      resultp [j + k]= bufp [j];
900  }
901  result.normalize();
902
903  return result;
904}
905
906zz_pEX kronSub (const CanonicalForm& A, int d, const Variable& alpha)
907{
908  int degAy= degree (A);
909  zz_pEX result;
910  result.rep.SetLength (d*(degAy + 1));
911
912  Variable v;
913  zz_pE *resultp;
914  resultp= result.rep.elts();
915  zz_pEX buf1;
916  zz_pE *buf1p;
917  zz_pX buf2;
918  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
919  int j, k, buf1RepLength;
920
921  for (CFIterator i= A; i.hasTerms(); i++)
922  {
923    if (i.coeff().inCoeffDomain())
924    {
925      buf2= convertFacCF2NTLzzpX (i.coeff());
926      buf1= to_zz_pEX (to_zz_pE (buf2));
927    }
928    else
929      buf1= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
930
931    k= i.exp()*d;
932    buf1p= buf1.rep.elts();
933    buf1RepLength= (int) buf1.rep.length();
934    for (j= 0; j < buf1RepLength; j++)
935      resultp [j + k]= buf1p [j];
936  }
937  result.normalize();
938
939  return result;
940}
941
942void
943kronSubRecipro (zz_pEX& subA1, zz_pEX& subA2,const CanonicalForm& A, int d,
944                const Variable& alpha)
945{
946  int degAy= degree (A);
947  subA1.rep.SetLength ((long) d*(degAy + 2));
948  subA2.rep.SetLength ((long) d*(degAy + 2));
949
950  Variable v;
951  zz_pE *subA1p;
952  zz_pE *subA2p;
953  subA1p= subA1.rep.elts();
954  subA2p= subA2.rep.elts();
955  zz_pEX buf;
956  zz_pE *bufp;
957  zz_pX buf2;
958  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
959  int j, k, kk, bufRepLength;
960
961  for (CFIterator i= A; i.hasTerms(); i++)
962  {
963    if (i.coeff().inCoeffDomain())
964    {
965      buf2= convertFacCF2NTLzzpX (i.coeff());
966      buf= to_zz_pEX (to_zz_pE (buf2));
967    }
968    else
969      buf= convertFacCF2NTLzz_pEX (i.coeff(), NTLMipo);
970
971    k= i.exp()*d;
972    kk= (degAy - i.exp())*d;
973    bufp= buf.rep.elts();
974    bufRepLength= (int) buf.rep.length();
975    for (j= 0; j < bufRepLength; j++)
976    {
977      subA1p [j + k] += bufp [j];
978      subA2p [j + kk] += bufp [j];
979    }
980  }
981  subA1.normalize();
982  subA2.normalize();
983}
984
985void
986kronSubRecipro (zz_pX& subA1, zz_pX& subA2,const CanonicalForm& A, int d)
987{
988  int degAy= degree (A);
989  subA1.rep.SetLength ((long) d*(degAy + 2));
990  subA2.rep.SetLength ((long) d*(degAy + 2));
991
992  zz_p *subA1p;
993  zz_p *subA2p;
994  subA1p= subA1.rep.elts();
995  subA2p= subA2.rep.elts();
996  zz_pX buf;
997  zz_p *bufp;
998  int j, k, kk, bufRepLength;
999
1000  for (CFIterator i= A; i.hasTerms(); i++)
1001  {
1002    buf= convertFacCF2NTLzzpX (i.coeff());
1003
1004    k= i.exp()*d;
1005    kk= (degAy - i.exp())*d;
1006    bufp= buf.rep.elts();
1007    bufRepLength= (int) buf.rep.length();
1008    for (j= 0; j < bufRepLength; j++)
1009    {
1010      subA1p [j + k] += bufp [j];
1011      subA2p [j + kk] += bufp [j];
1012    }
1013  }
1014  subA1.normalize();
1015  subA2.normalize();
1016}
1017
1018CanonicalForm
1019reverseSubst (const zz_pEX& F, const zz_pEX& G, int d, int k,
1020              const Variable& alpha)
1021{
1022  Variable y= Variable (2);
1023  Variable x= Variable (1);
1024
1025  zz_pEX f= F;
1026  zz_pEX g= G;
1027  int degf= deg(f);
1028  int degg= deg(g);
1029
1030  zz_pEX buf1;
1031  zz_pEX buf2;
1032  zz_pEX buf3;
1033
1034  zz_pE *buf1p;
1035  zz_pE *buf2p;
1036  zz_pE *buf3p;
1037  if (f.rep.length() < (long) d*(k+1)) //zero padding
1038    f.rep.SetLength ((long)d*(k+1));
1039
1040  zz_pE *gp= g.rep.elts();
1041  zz_pE *fp= f.rep.elts();
1042  CanonicalForm result= 0;
1043  int i= 0;
1044  int lf= 0;
1045  int lg= d*k;
1046  int degfSubLf= degf;
1047  int deggSubLg= degg-lg;
1048  int repLengthBuf2, repLengthBuf1, ind, tmp;
1049  zz_pE zzpEZero= zz_pE();
1050
1051  while (degf >= lf || lg >= 0)
1052  {
1053    if (degfSubLf >= d)
1054      repLengthBuf1= d;
1055    else if (degfSubLf < 0)
1056      repLengthBuf1= 0;
1057    else
1058      repLengthBuf1= degfSubLf + 1;
1059    buf1.rep.SetLength((long) repLengthBuf1);
1060
1061    buf1p= buf1.rep.elts();
1062    for (ind= 0; ind < repLengthBuf1; ind++)
1063      buf1p [ind]= fp [ind + lf];
1064    buf1.normalize();
1065
1066    repLengthBuf1= buf1.rep.length();
1067
1068    if (deggSubLg >= d - 1)
1069      repLengthBuf2= d - 1;
1070    else if (deggSubLg < 0)
1071      repLengthBuf2= 0;
1072    else
1073      repLengthBuf2= deggSubLg + 1;
1074
1075    buf2.rep.SetLength ((long) repLengthBuf2);
1076    buf2p= buf2.rep.elts();
1077    for (ind= 0; ind < repLengthBuf2; ind++)
1078      buf2p [ind]= gp [ind + lg];
1079    buf2.normalize();
1080
1081    repLengthBuf2= buf2.rep.length();
1082
1083    buf3.rep.SetLength((long) repLengthBuf2 + d);
1084    buf3p= buf3.rep.elts();
1085    buf2p= buf2.rep.elts();
1086    buf1p= buf1.rep.elts();
1087    for (ind= 0; ind < repLengthBuf1; ind++)
1088      buf3p [ind]= buf1p [ind];
1089    for (ind= repLengthBuf1; ind < d; ind++)
1090      buf3p [ind]= zzpEZero;
1091    for (ind= 0; ind < repLengthBuf2; ind++)
1092      buf3p [ind + d]= buf2p [ind];
1093    buf3.normalize();
1094
1095    result += convertNTLzz_pEX2CF (buf3, x, alpha)*power (y, i);
1096    i++;
1097
1098
1099    lf= i*d;
1100    degfSubLf= degf - lf;
1101
1102    lg= d*(k-i);
1103    deggSubLg= degg - lg;
1104
1105    buf1p= buf1.rep.elts();
1106
1107    if (lg >= 0 && deggSubLg > 0)
1108    {
1109      if (repLengthBuf2 > degfSubLf + 1)
1110        degfSubLf= repLengthBuf2 - 1;
1111      tmp= tmin (repLengthBuf1, deggSubLg + 1);
1112      for (ind= 0; ind < tmp; ind++)
1113        gp [ind + lg] -= buf1p [ind];
1114    }
1115
1116    if (lg < 0)
1117      break;
1118
1119    buf2p= buf2.rep.elts();
1120    if (degfSubLf >= 0)
1121    {
1122      for (ind= 0; ind < repLengthBuf2; ind++)
1123        fp [ind + lf] -= buf2p [ind];
1124    }
1125  }
1126
1127  return result;
1128}
1129
1130CanonicalForm
1131reverseSubst (const zz_pX& F, const zz_pX& G, int d, int k)
1132{
1133  Variable y= Variable (2);
1134  Variable x= Variable (1);
1135
1136  zz_pX f= F;
1137  zz_pX g= G;
1138  int degf= deg(f);
1139  int degg= deg(g);
1140
1141  zz_pX buf1;
1142  zz_pX buf2;
1143  zz_pX buf3;
1144
1145  zz_p *buf1p;
1146  zz_p *buf2p;
1147  zz_p *buf3p;
1148
1149  if (f.rep.length() < (long) d*(k+1)) //zero padding
1150    f.rep.SetLength ((long)d*(k+1));
1151
1152  zz_p *gp= g.rep.elts();
1153  zz_p *fp= f.rep.elts();
1154  CanonicalForm result= 0;
1155  int i= 0;
1156  int lf= 0;
1157  int lg= d*k;
1158  int degfSubLf= degf;
1159  int deggSubLg= degg-lg;
1160  int repLengthBuf2, repLengthBuf1, ind, tmp;
1161  zz_p zzpZero= zz_p();
1162  while (degf >= lf || lg >= 0)
1163  {
1164    if (degfSubLf >= d)
1165      repLengthBuf1= d;
1166    else if (degfSubLf < 0)
1167      repLengthBuf1= 0;
1168    else
1169      repLengthBuf1= degfSubLf + 1;
1170    buf1.rep.SetLength((long) repLengthBuf1);
1171
1172    buf1p= buf1.rep.elts();
1173    for (ind= 0; ind < repLengthBuf1; ind++)
1174      buf1p [ind]= fp [ind + lf];
1175    buf1.normalize();
1176
1177    repLengthBuf1= buf1.rep.length();
1178
1179    if (deggSubLg >= d - 1)
1180      repLengthBuf2= d - 1;
1181    else if (deggSubLg < 0)
1182      repLengthBuf2= 0;
1183    else
1184      repLengthBuf2= deggSubLg + 1;
1185
1186    buf2.rep.SetLength ((long) repLengthBuf2);
1187    buf2p= buf2.rep.elts();
1188    for (ind= 0; ind < repLengthBuf2; ind++)
1189      buf2p [ind]= gp [ind + lg];
1190
1191    buf2.normalize();
1192
1193    repLengthBuf2= buf2.rep.length();
1194
1195
1196    buf3.rep.SetLength((long) repLengthBuf2 + d);
1197    buf3p= buf3.rep.elts();
1198    buf2p= buf2.rep.elts();
1199    buf1p= buf1.rep.elts();
1200    for (ind= 0; ind < repLengthBuf1; ind++)
1201      buf3p [ind]= buf1p [ind];
1202    for (ind= repLengthBuf1; ind < d; ind++)
1203      buf3p [ind]= zzpZero;
1204    for (ind= 0; ind < repLengthBuf2; ind++)
1205      buf3p [ind + d]= buf2p [ind];
1206    buf3.normalize();
1207
1208    result += convertNTLzzpX2CF (buf3, x)*power (y, i);
1209    i++;
1210
1211
1212    lf= i*d;
1213    degfSubLf= degf - lf;
1214
1215    lg= d*(k-i);
1216    deggSubLg= degg - lg;
1217
1218    buf1p= buf1.rep.elts();
1219
1220    if (lg >= 0 && deggSubLg > 0)
1221    {
1222      if (repLengthBuf2 > degfSubLf + 1)
1223        degfSubLf= repLengthBuf2 - 1;
1224      tmp= tmin (repLengthBuf1, deggSubLg + 1);
1225      for (ind= 0; ind < tmp; ind++)
1226        gp [ind + lg] -= buf1p [ind];
1227    }
1228    if (lg < 0)
1229      break;
1230
1231    buf2p= buf2.rep.elts();
1232    if (degfSubLf >= 0)
1233    {
1234      for (ind= 0; ind < repLengthBuf2; ind++)
1235        fp [ind + lf] -= buf2p [ind];
1236    }
1237  }
1238
1239  return result;
1240}
1241
1242CanonicalForm reverseSubst (const zz_pEX& F, int d, const Variable& alpha)
1243{
1244  Variable y= Variable (2);
1245  Variable x= Variable (1);
1246
1247  zz_pEX f= F;
1248  zz_pE *fp= f.rep.elts();
1249
1250  zz_pEX buf;
1251  zz_pE *bufp;
1252  CanonicalForm result= 0;
1253  int i= 0;
1254  int degf= deg(f);
1255  int k= 0;
1256  int degfSubK, repLength, j;
1257  while (degf >= k)
1258  {
1259    degfSubK= degf - k;
1260    if (degfSubK >= d)
1261      repLength= d;
1262    else
1263      repLength= degfSubK + 1;
1264
1265    buf.rep.SetLength ((long) repLength);
1266    bufp= buf.rep.elts();
1267    for (j= 0; j < repLength; j++)
1268      bufp [j]= fp [j + k];
1269    buf.normalize();
1270
1271    result += convertNTLzz_pEX2CF (buf, x, alpha)*power (y, i);
1272    i++;
1273    k= d*i;
1274  }
1275
1276  return result;
1277}
1278
1279CanonicalForm reverseSubstFp (const zz_pX& F, int d)
1280{
1281  Variable y= Variable (2);
1282  Variable x= Variable (1);
1283
1284  zz_pX f= F;
1285  zz_p *fp= f.rep.elts();
1286
1287  zz_pX buf;
1288  zz_p *bufp;
1289  CanonicalForm result= 0;
1290  int i= 0;
1291  int degf= deg(f);
1292  int k= 0;
1293  int degfSubK, repLength, j;
1294  while (degf >= k)
1295  {
1296    degfSubK= degf - k;
1297    if (degfSubK >= d)
1298      repLength= d;
1299    else
1300      repLength= degfSubK + 1;
1301
1302    buf.rep.SetLength ((long) repLength);
1303    bufp= buf.rep.elts();
1304    for (j= 0; j < repLength; j++)
1305      bufp [j]= fp [j + k];
1306    buf.normalize();
1307
1308    result += convertNTLzzpX2CF (buf, x)*power (y, i);
1309    i++;
1310    k= d*i;
1311  }
1312
1313  return result;
1314}
1315
1316// assumes input to be reduced mod M and to be an element of Fq not Fp
1317CanonicalForm
1318mulMod2NTLFpReci (const CanonicalForm& F, const CanonicalForm& G, const
1319                  CanonicalForm& M)
1320{
1321  int d1= degree (F, 1) + degree (G, 1) + 1;
1322  d1 /= 2;
1323  d1 += 1;
1324
1325  zz_pX F1, F2;
1326  kronSubRecipro (F1, F2, F, d1);
1327  zz_pX G1, G2;
1328  kronSubRecipro (G1, G2, G, d1);
1329
1330  int k= d1*degree (M);
1331  MulTrunc (F1, F1, G1, (long) k);
1332
1333  int degtailF= degree (tailcoeff (F), 1);
1334  int degtailG= degree (tailcoeff (G), 1);
1335  int taildegF= taildegree (F);
1336  int taildegG= taildegree (G);
1337  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
1338
1339  reverse (F2, F2);
1340  reverse (G2, G2);
1341  MulTrunc (F2, F2, G2, b + 1);
1342  reverse (F2, F2, b);
1343
1344  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
1345  return reverseSubst (F1, F2, d1, d2);
1346}
1347
1348//Kronecker substitution
1349CanonicalForm
1350mulMod2NTLFp (const CanonicalForm& F, const CanonicalForm& G, const
1351              CanonicalForm& M)
1352{
1353  CanonicalForm A= F;
1354  CanonicalForm B= G;
1355
1356  int degAx= degree (A, 1);
1357  int degAy= degree (A, 2);
1358  int degBx= degree (B, 1);
1359  int degBy= degree (B, 2);
1360  int d1= degAx + 1 + degBx;
1361  int d2= tmax (degAy, degBy);
1362
1363  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
1364    return mulMod2NTLFpReci (A, B, M);
1365
1366  zz_pX NTLA= kronSubFp (A, d1);
1367  zz_pX NTLB= kronSubFp (B, d1);
1368
1369  int k= d1*degree (M);
1370  MulTrunc (NTLA, NTLA, NTLB, (long) k);
1371
1372  A= reverseSubstFp (NTLA, d1);
1373
1374  return A;
1375}
1376
1377// assumes input to be reduced mod M and to be an element of Fq not Fp
1378CanonicalForm
1379mulMod2NTLFqReci (const CanonicalForm& F, const CanonicalForm& G, const
1380                  CanonicalForm& M, const Variable& alpha)
1381{
1382  int d1= degree (F, 1) + degree (G, 1) + 1;
1383  d1 /= 2;
1384  d1 += 1;
1385
1386  zz_pEX F1, F2;
1387  kronSubRecipro (F1, F2, F, d1, alpha);
1388  zz_pEX G1, G2;
1389  kronSubRecipro (G1, G2, G, d1, alpha);
1390
1391  int k= d1*degree (M);
1392  MulTrunc (F1, F1, G1, (long) k);
1393
1394  int degtailF= degree (tailcoeff (F), 1);
1395  int degtailG= degree (tailcoeff (G), 1);
1396  int taildegF= taildegree (F);
1397  int taildegG= taildegree (G);
1398  int b= k + degtailF + degtailG - d1*(2+taildegF+taildegG);
1399
1400  reverse (F2, F2);
1401  reverse (G2, G2);
1402  MulTrunc (F2, F2, G2, b + 1);
1403  reverse (F2, F2, b);
1404
1405  int d2= tmax (deg (F2)/d1, deg (F1)/d1);
1406  return reverseSubst (F1, F2, d1, d2, alpha);
1407}
1408
1409#ifdef HAVE_FLINT
1410CanonicalForm
1411mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
1412                CanonicalForm& M);
1413#endif
1414
1415CanonicalForm
1416mulMod2NTLFq (const CanonicalForm& F, const CanonicalForm& G, const
1417              CanonicalForm& M)
1418{
1419  Variable alpha;
1420  CanonicalForm A= F;
1421  CanonicalForm B= G;
1422
1423  if (hasFirstAlgVar (A, alpha) || hasFirstAlgVar (B, alpha))
1424  {
1425    int degAx= degree (A, 1);
1426    int degAy= degree (A, 2);
1427    int degBx= degree (B, 1);
1428    int degBy= degree (B, 2);
1429    int d1= degAx + degBx + 1;
1430    int d2= tmax (degAy, degBy);
1431    zz_p::init (getCharacteristic());
1432    zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1433    zz_pE::init (NTLMipo);
1434
1435    int degMipo= degree (getMipo (alpha));
1436    if ((d1 > 128/degMipo) && (d2 > 160/degMipo) && (degAy == degBy) &&
1437        (2*degAy > degree (M)))
1438      return mulMod2NTLFqReci (A, B, M, alpha);
1439
1440    zz_pEX NTLA= kronSub (A, d1, alpha);
1441    zz_pEX NTLB= kronSub (B, d1, alpha);
1442
1443    int k= d1*degree (M);
1444
1445    MulTrunc (NTLA, NTLA, NTLB, (long) k);
1446
1447    A= reverseSubst (NTLA, d1, alpha);
1448
1449    return A;
1450  }
1451  else
1452#ifdef HAVE_FLINT
1453    return mulMod2FLINTFp (A, B, M);
1454#else
1455    return mulMod2NTLFp (A, B, M);
1456#endif
1457}
1458
1459#ifdef HAVE_FLINT
1460void
1461kronSubRecipro (nmod_poly_t subA1, nmod_poly_t subA2, const CanonicalForm& A,
1462                int d)
1463{
1464  int degAy= degree (A);
1465  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1466  nmod_poly_init2_preinv (subA1, getCharacteristic(), ninv, d*(degAy + 2));
1467  nmod_poly_init2_preinv (subA2, getCharacteristic(), ninv, d*(degAy + 2));
1468
1469  nmod_poly_t buf;
1470
1471  int k, kk, j, bufRepLength;
1472  for (CFIterator i= A; i.hasTerms(); i++)
1473  {
1474    convertFacCF2nmod_poly_t (buf, i.coeff());
1475
1476    k= i.exp()*d;
1477    kk= (degAy - i.exp())*d;
1478    bufRepLength= (int) nmod_poly_length (buf);
1479    for (j= 0; j < bufRepLength; j++)
1480    {
1481      nmod_poly_set_coeff_ui (subA1, j + k,
1482                              n_addmod (nmod_poly_get_coeff_ui (subA1, j+k),
1483                                        nmod_poly_get_coeff_ui (buf, j),
1484                                        getCharacteristic()
1485                                       )
1486                             );
1487      nmod_poly_set_coeff_ui (subA2, j + kk,
1488                              n_addmod (nmod_poly_get_coeff_ui (subA2, j + kk),
1489                                        nmod_poly_get_coeff_ui (buf, j),
1490                                        getCharacteristic()
1491                                       )
1492                             );
1493    }
1494    nmod_poly_clear (buf);
1495  }
1496  _nmod_poly_normalise (subA1);
1497  _nmod_poly_normalise (subA2);
1498}
1499
1500void
1501kronSubRecipro (fmpz_poly_t subA1, fmpz_poly_t subA2, const CanonicalForm& A,
1502                int d)
1503{
1504  int degAy= degree (A);
1505  fmpz_poly_init2 (subA1, d*(degAy + 2));
1506  fmpz_poly_init2 (subA2, d*(degAy + 2));
1507
1508  fmpz_poly_t buf;
1509  fmpz_t coeff1, coeff2;
1510
1511  int k, kk, j, bufRepLength;
1512  for (CFIterator i= A; i.hasTerms(); i++)
1513  {
1514    convertFacCF2Fmpz_poly_t (buf, i.coeff());
1515
1516    k= i.exp()*d;
1517    kk= (degAy - i.exp())*d;
1518    bufRepLength= (int) fmpz_poly_length (buf);
1519    for (j= 0; j < bufRepLength; j++)
1520    {
1521      fmpz_poly_get_coeff_fmpz (coeff1, subA1, j+k);
1522      fmpz_poly_get_coeff_fmpz (coeff2, buf, j);
1523      fmpz_add (coeff1, coeff1, coeff2);
1524      fmpz_poly_set_coeff_fmpz (subA1, j + k, coeff1);
1525      fmpz_poly_get_coeff_fmpz (coeff1, subA2, j + kk);
1526      fmpz_add (coeff1, coeff1, coeff2);
1527      fmpz_poly_set_coeff_fmpz (subA2, j + kk, coeff1);
1528    }
1529    fmpz_poly_clear (buf);
1530  }
1531  fmpz_clear (coeff1);
1532  fmpz_clear (coeff2);
1533  _fmpz_poly_normalise (subA1);
1534  _fmpz_poly_normalise (subA2);
1535}
1536
1537CanonicalForm reverseSubstQ (const fmpz_poly_t F, int d)
1538{
1539  Variable y= Variable (2);
1540  Variable x= Variable (1);
1541
1542  fmpz_poly_t f;
1543  fmpz_poly_init (f);
1544  fmpz_poly_set (f, F);
1545
1546  fmpz_poly_t buf;
1547  CanonicalForm result= 0;
1548  int i= 0;
1549  int degf= fmpz_poly_degree(f);
1550  int k= 0;
1551  int degfSubK, repLength, j;
1552  fmpz_t coeff;
1553  while (degf >= k)
1554  {
1555    degfSubK= degf - k;
1556    if (degfSubK >= d)
1557      repLength= d;
1558    else
1559      repLength= degfSubK + 1;
1560
1561    fmpz_poly_init2 (buf, repLength);
1562    fmpz_init (coeff);
1563    for (j= 0; j < repLength; j++)
1564    {
1565      fmpz_poly_get_coeff_fmpz (coeff, f, j + k);
1566      fmpz_poly_set_coeff_fmpz (buf, j, coeff);
1567    }
1568    _fmpz_poly_normalise (buf);
1569
1570    result += convertFmpz_poly_t2FacCF (buf, x)*power (y, i);
1571    i++;
1572    k= d*i;
1573    fmpz_poly_clear (buf);
1574    fmpz_clear (coeff);
1575  }
1576  fmpz_poly_clear (f);
1577
1578  return result;
1579}
1580
1581CanonicalForm
1582reverseSubst (const nmod_poly_t F, const nmod_poly_t G, int d, int k)
1583{
1584  Variable y= Variable (2);
1585  Variable x= Variable (1);
1586
1587  nmod_poly_t f, g;
1588  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1589  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
1590  nmod_poly_init_preinv (g, getCharacteristic(), ninv);
1591  nmod_poly_set (f, F);
1592  nmod_poly_set (g, G);
1593  int degf= nmod_poly_degree(f);
1594  int degg= nmod_poly_degree(g);
1595
1596
1597  nmod_poly_t buf1,buf2, buf3;
1598
1599  if (nmod_poly_length (f) < (long) d*(k+1)) //zero padding
1600    nmod_poly_fit_length (f,(long)d*(k+1));
1601
1602  CanonicalForm result= 0;
1603  int i= 0;
1604  int lf= 0;
1605  int lg= d*k;
1606  int degfSubLf= degf;
1607  int deggSubLg= degg-lg;
1608  int repLengthBuf2, repLengthBuf1, ind, tmp;
1609  while (degf >= lf || lg >= 0)
1610  {
1611    if (degfSubLf >= d)
1612      repLengthBuf1= d;
1613    else if (degfSubLf < 0)
1614      repLengthBuf1= 0;
1615    else
1616      repLengthBuf1= degfSubLf + 1;
1617    nmod_poly_init2_preinv (buf1, getCharacteristic(), ninv, repLengthBuf1);
1618
1619    for (ind= 0; ind < repLengthBuf1; ind++)
1620      nmod_poly_set_coeff_ui (buf1, ind, nmod_poly_get_coeff_ui (f, ind+lf));
1621    _nmod_poly_normalise (buf1);
1622
1623    repLengthBuf1= nmod_poly_length (buf1);
1624
1625    if (deggSubLg >= d - 1)
1626      repLengthBuf2= d - 1;
1627    else if (deggSubLg < 0)
1628      repLengthBuf2= 0;
1629    else
1630      repLengthBuf2= deggSubLg + 1;
1631
1632    nmod_poly_init2_preinv (buf2, getCharacteristic(), ninv, repLengthBuf2);
1633    for (ind= 0; ind < repLengthBuf2; ind++)
1634      nmod_poly_set_coeff_ui (buf2, ind, nmod_poly_get_coeff_ui (g, ind + lg));
1635
1636    _nmod_poly_normalise (buf2);
1637    repLengthBuf2= nmod_poly_length (buf2);
1638
1639    nmod_poly_init2_preinv (buf3, getCharacteristic(), ninv, repLengthBuf2 + d);
1640    for (ind= 0; ind < repLengthBuf1; ind++)
1641      nmod_poly_set_coeff_ui (buf3, ind, nmod_poly_get_coeff_ui (buf1, ind));
1642    for (ind= repLengthBuf1; ind < d; ind++)
1643      nmod_poly_set_coeff_ui (buf3, ind, 0);
1644    for (ind= 0; ind < repLengthBuf2; ind++)
1645      nmod_poly_set_coeff_ui (buf3, ind+d, nmod_poly_get_coeff_ui (buf2, ind));
1646    _nmod_poly_normalise (buf3);
1647
1648    result += convertnmod_poly_t2FacCF (buf3, x)*power (y, i);
1649    i++;
1650
1651
1652    lf= i*d;
1653    degfSubLf= degf - lf;
1654
1655    lg= d*(k-i);
1656    deggSubLg= degg - lg;
1657
1658    if (lg >= 0 && deggSubLg > 0)
1659    {
1660      if (repLengthBuf2 > degfSubLf + 1)
1661        degfSubLf= repLengthBuf2 - 1;
1662      tmp= tmin (repLengthBuf1, deggSubLg + 1);
1663      for (ind= 0; ind < tmp; ind++)
1664        nmod_poly_set_coeff_ui (g, ind + lg,
1665                                n_submod (nmod_poly_get_coeff_ui (g, ind + lg),
1666                                          nmod_poly_get_coeff_ui (buf1, ind),
1667                                          getCharacteristic()
1668                                         )
1669                               );
1670    }
1671    if (lg < 0)
1672    {
1673      nmod_poly_clear (buf1);
1674      nmod_poly_clear (buf2);
1675      nmod_poly_clear (buf3);
1676      break;
1677    }
1678    if (degfSubLf >= 0)
1679    {
1680      for (ind= 0; ind < repLengthBuf2; ind++)
1681        nmod_poly_set_coeff_ui (f, ind + lf,
1682                                n_submod (nmod_poly_get_coeff_ui (f, ind + lf),
1683                                          nmod_poly_get_coeff_ui (buf2, ind),
1684                                          getCharacteristic()
1685                                         )
1686                               );
1687    }
1688    nmod_poly_clear (buf1);
1689    nmod_poly_clear (buf2);
1690    nmod_poly_clear (buf3);
1691  }
1692
1693  nmod_poly_clear (f);
1694  nmod_poly_clear (g);
1695
1696  return result;
1697}
1698
1699CanonicalForm
1700reverseSubst (const fmpz_poly_t F, const fmpz_poly_t G, int d, int k)
1701{
1702  Variable y= Variable (2);
1703  Variable x= Variable (1);
1704
1705  fmpz_poly_t f, g;
1706  fmpz_poly_init (f);
1707  fmpz_poly_init (g);
1708  fmpz_poly_set (f, F);
1709  fmpz_poly_set (g, G);
1710  int degf= fmpz_poly_degree(f);
1711  int degg= fmpz_poly_degree(g);
1712
1713
1714  fmpz_poly_t buf1,buf2, buf3;
1715
1716  if (fmpz_poly_length (f) < (long) d*(k+1)) //zero padding
1717    fmpz_poly_fit_length (f,(long)d*(k+1));
1718
1719  CanonicalForm result= 0;
1720  int i= 0;
1721  int lf= 0;
1722  int lg= d*k;
1723  int degfSubLf= degf;
1724  int deggSubLg= degg-lg;
1725  int repLengthBuf2, repLengthBuf1, ind, tmp;
1726  fmpz_t tmp1, tmp2;
1727  while (degf >= lf || lg >= 0)
1728  {
1729    if (degfSubLf >= d)
1730      repLengthBuf1= d;
1731    else if (degfSubLf < 0)
1732      repLengthBuf1= 0;
1733    else
1734      repLengthBuf1= degfSubLf + 1;
1735
1736    fmpz_poly_init2 (buf1, repLengthBuf1);
1737
1738    for (ind= 0; ind < repLengthBuf1; ind++)
1739    {
1740      fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1741      fmpz_poly_set_coeff_fmpz (buf1, ind, tmp1);
1742    }
1743    _fmpz_poly_normalise (buf1);
1744
1745    repLengthBuf1= fmpz_poly_length (buf1);
1746
1747    if (deggSubLg >= d - 1)
1748      repLengthBuf2= d - 1;
1749    else if (deggSubLg < 0)
1750      repLengthBuf2= 0;
1751    else
1752      repLengthBuf2= deggSubLg + 1;
1753
1754    fmpz_poly_init2 (buf2, repLengthBuf2);
1755
1756    for (ind= 0; ind < repLengthBuf2; ind++)
1757    {
1758      fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1759      fmpz_poly_set_coeff_fmpz (buf2, ind, tmp1);
1760    }
1761
1762    _fmpz_poly_normalise (buf2);
1763    repLengthBuf2= fmpz_poly_length (buf2);
1764
1765    fmpz_poly_init2 (buf3, repLengthBuf2 + d);
1766    for (ind= 0; ind < repLengthBuf1; ind++)
1767    {
1768      fmpz_poly_get_coeff_fmpz (tmp1, buf1, ind);  //oder fmpz_set (fmpz_poly_get_coeff_ptr (buf3, ind),fmpz_poly_get_coeff_ptr (buf1, ind))
1769      fmpz_poly_set_coeff_fmpz (buf3, ind, tmp1);
1770    }
1771    for (ind= repLengthBuf1; ind < d; ind++)
1772      fmpz_poly_set_coeff_ui (buf3, ind, 0);
1773    for (ind= 0; ind < repLengthBuf2; ind++)
1774    {
1775      fmpz_poly_get_coeff_fmpz (tmp1, buf2, ind);
1776      fmpz_poly_set_coeff_fmpz (buf3, ind + d, tmp1);
1777    }
1778    _fmpz_poly_normalise (buf3);
1779
1780    result += convertFmpz_poly_t2FacCF (buf3, x)*power (y, i);
1781    i++;
1782
1783
1784    lf= i*d;
1785    degfSubLf= degf - lf;
1786
1787    lg= d*(k-i);
1788    deggSubLg= degg - lg;
1789
1790    if (lg >= 0 && deggSubLg > 0)
1791    {
1792      if (repLengthBuf2 > degfSubLf + 1)
1793        degfSubLf= repLengthBuf2 - 1;
1794      tmp= tmin (repLengthBuf1, deggSubLg + 1);
1795      for (ind= 0; ind < tmp; ind++)
1796      {
1797        fmpz_poly_get_coeff_fmpz (tmp1, g, ind + lg);
1798        fmpz_poly_get_coeff_fmpz (tmp2, buf1, ind);
1799        fmpz_sub (tmp1, tmp1, tmp2);
1800        fmpz_poly_set_coeff_fmpz (g, ind + lg, tmp1);
1801      }
1802    }
1803    if (lg < 0)
1804    {
1805      fmpz_poly_clear (buf1);
1806      fmpz_poly_clear (buf2);
1807      fmpz_poly_clear (buf3);
1808      break;
1809    }
1810    if (degfSubLf >= 0)
1811    {
1812      for (ind= 0; ind < repLengthBuf2; ind++)
1813      {
1814        fmpz_poly_get_coeff_fmpz (tmp1, f, ind + lf);
1815        fmpz_poly_get_coeff_fmpz (tmp2, buf2, ind);
1816        fmpz_sub (tmp1, tmp1, tmp2);
1817        fmpz_poly_set_coeff_fmpz (f, ind + lf, tmp1);
1818      }
1819    }
1820    fmpz_poly_clear (buf1);
1821    fmpz_poly_clear (buf2);
1822    fmpz_poly_clear (buf3);
1823  }
1824
1825  fmpz_poly_clear (f);
1826  fmpz_poly_clear (g);
1827  fmpz_clear (tmp1);
1828  fmpz_clear (tmp2);
1829
1830  return result;
1831}
1832
1833CanonicalForm reverseSubstFp (const nmod_poly_t F, int d)
1834{
1835  Variable y= Variable (2);
1836  Variable x= Variable (1);
1837
1838  nmod_poly_t f;
1839  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1840  nmod_poly_init_preinv (f, getCharacteristic(), ninv);
1841  nmod_poly_set (f, F);
1842
1843  nmod_poly_t buf;
1844  CanonicalForm result= 0;
1845  int i= 0;
1846  int degf= nmod_poly_degree(f);
1847  int k= 0;
1848  int degfSubK, repLength, j;
1849  while (degf >= k)
1850  {
1851    degfSubK= degf - k;
1852    if (degfSubK >= d)
1853      repLength= d;
1854    else
1855      repLength= degfSubK + 1;
1856
1857    nmod_poly_init2_preinv (buf, getCharacteristic(), ninv, repLength);
1858    for (j= 0; j < repLength; j++)
1859      nmod_poly_set_coeff_ui (buf, j, nmod_poly_get_coeff_ui (f, j + k));
1860    _nmod_poly_normalise (buf);
1861
1862    result += convertnmod_poly_t2FacCF (buf, x)*power (y, i);
1863    i++;
1864    k= d*i;
1865    nmod_poly_clear (buf);
1866  }
1867  nmod_poly_clear (f);
1868
1869  return result;
1870}
1871
1872CanonicalForm
1873mulMod2FLINTFpReci (const CanonicalForm& F, const CanonicalForm& G, const
1874                    CanonicalForm& M)
1875{
1876  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
1877  d1 /= 2;
1878  d1 += 1;
1879
1880  nmod_poly_t F1, F2;
1881  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1882  nmod_poly_init_preinv (F1, getCharacteristic(), ninv);
1883  nmod_poly_init_preinv (F2, getCharacteristic(), ninv);
1884  kronSubRecipro (F1, F2, F, d1);
1885
1886  nmod_poly_t G1, G2;
1887  nmod_poly_init_preinv (G1, getCharacteristic(), ninv);
1888  nmod_poly_init_preinv (G2, getCharacteristic(), ninv);
1889  kronSubRecipro (G1, G2, G, d1);
1890
1891  int k= d1*degree (M);
1892  nmod_poly_mullow (F1, F1, G1, (long) k);
1893
1894  int degtailF= degree (tailcoeff (F), 1);;
1895  int degtailG= degree (tailcoeff (G), 1);
1896  int taildegF= taildegree (F);
1897  int taildegG= taildegree (G);
1898
1899  int b= nmod_poly_degree (F2) + nmod_poly_degree (G2) - k - degtailF - degtailG
1900         + d1*(2+taildegF + taildegG);
1901  nmod_poly_mulhigh (F2, F2, G2, b);
1902  nmod_poly_shift_right (F2, F2, b);
1903  int d2= tmax (nmod_poly_degree (F2)/d1, nmod_poly_degree (F1)/d1);
1904
1905
1906  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
1907
1908  nmod_poly_clear (F1);
1909  nmod_poly_clear (F2);
1910  nmod_poly_clear (G1);
1911  nmod_poly_clear (G2);
1912  return result;
1913}
1914
1915CanonicalForm
1916mulMod2FLINTFp (const CanonicalForm& F, const CanonicalForm& G, const
1917                CanonicalForm& M)
1918{
1919  CanonicalForm A= F;
1920  CanonicalForm B= G;
1921
1922  int degAx= degree (A, 1);
1923  int degAy= degree (A, 2);
1924  int degBx= degree (B, 1);
1925  int degBy= degree (B, 2);
1926  int d1= degAx + 1 + degBx;
1927  int d2= tmax (degAy, degBy);
1928
1929  if (d1 > 128 && d2 > 160 && (degAy == degBy) && (2*degAy > degree (M)))
1930    return mulMod2FLINTFpReci (A, B, M);
1931
1932  nmod_poly_t FLINTA, FLINTB;
1933  mp_limb_t ninv= n_preinvert_limb (getCharacteristic());
1934  nmod_poly_init_preinv (FLINTA, getCharacteristic(), ninv);
1935  nmod_poly_init_preinv (FLINTB, getCharacteristic(), ninv);
1936  kronSubFp (FLINTA, A, d1);
1937  kronSubFp (FLINTB, B, d1);
1938
1939  int k= d1*degree (M);
1940  nmod_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
1941
1942  A= reverseSubstFp (FLINTA, d1);
1943
1944  nmod_poly_clear (FLINTA);
1945  nmod_poly_clear (FLINTB);
1946  return A;
1947}
1948
1949CanonicalForm
1950mulMod2FLINTQReci (const CanonicalForm& F, const CanonicalForm& G, const
1951                    CanonicalForm& M)
1952{
1953  int d1= tmax (degree (F, 1), degree (G, 1)) + 1;
1954  d1 /= 2;
1955  d1 += 1;
1956
1957  fmpz_poly_t F1, F2;
1958  fmpz_poly_init (F1);
1959  fmpz_poly_init (F2);
1960  kronSubRecipro (F1, F2, F, d1);
1961
1962  fmpz_poly_t G1, G2;
1963  fmpz_poly_init (G1);
1964  fmpz_poly_init (G2);
1965  kronSubRecipro (G1, G2, G, d1);
1966
1967  int k= d1*degree (M);
1968  fmpz_poly_mullow (F1, F1, G1, (long) k);
1969
1970  int degtailF= degree (tailcoeff (F), 1);;
1971  int degtailG= degree (tailcoeff (G), 1);
1972  int taildegF= taildegree (F);
1973  int taildegG= taildegree (G);
1974
1975  int b= fmpz_poly_degree (F2) + fmpz_poly_degree (G2) - k - degtailF - degtailG
1976         + d1*(2+taildegF + taildegG);
1977  fmpz_poly_mulhigh_n (F2, F2, G2, b);
1978  fmpz_poly_shift_right (F2, F2, b);
1979  int d2= tmax (fmpz_poly_degree (F2)/d1, fmpz_poly_degree (F1)/d1);
1980
1981  CanonicalForm result= reverseSubst (F1, F2, d1, d2);
1982
1983  fmpz_poly_clear (F1);
1984  fmpz_poly_clear (F2);
1985  fmpz_poly_clear (G1);
1986  fmpz_poly_clear (G2);
1987  return result;
1988}
1989
1990CanonicalForm
1991mulMod2FLINTQ (const CanonicalForm& F, const CanonicalForm& G, const
1992                CanonicalForm& M)
1993{
1994  CanonicalForm A= F;
1995  CanonicalForm B= G;
1996
1997  int degAx= degree (A, 1);
1998  int degBx= degree (B, 1);
1999  int d1= degAx + 1 + degBx;
2000
2001  CanonicalForm f= bCommonDen (F);
2002  CanonicalForm g= bCommonDen (G);
2003  A *= f;
2004  B *= g;
2005
2006  fmpz_poly_t FLINTA, FLINTB;
2007  fmpz_poly_init (FLINTA);
2008  fmpz_poly_init (FLINTB);
2009  kronSub (FLINTA, A, d1);
2010  kronSub (FLINTB, B, d1);
2011  int k= d1*degree (M);
2012
2013  fmpz_poly_mullow (FLINTA, FLINTA, FLINTB, (long) k);
2014  A= reverseSubstQ (FLINTA, d1);
2015  fmpz_poly_clear (FLINTA);
2016  fmpz_poly_clear (FLINTB);
2017  return A/(f*g);
2018}
2019
2020#endif
2021
2022CanonicalForm mulMod2 (const CanonicalForm& A, const CanonicalForm& B,
2023                       const CanonicalForm& M)
2024{
2025  if (A.isZero() || B.isZero())
2026    return 0;
2027
2028  ASSERT (M.isUnivariate(), "M must be univariate");
2029
2030  CanonicalForm F= mod (A, M);
2031  CanonicalForm G= mod (B, M);
2032  if (F.inCoeffDomain() || G.inCoeffDomain())
2033    return F*G;
2034  Variable y= M.mvar();
2035  int degF= degree (F, y);
2036  int degG= degree (G, y);
2037
2038  if ((degF < 1 && degG < 1) && (F.isUnivariate() && G.isUnivariate()) &&
2039      (F.level() == G.level()))
2040  {
2041    CanonicalForm result= mulNTL (F, G);
2042    return mod (result, M);
2043  }
2044  else if (degF <= 1 && degG <= 1)
2045  {
2046    CanonicalForm result= F*G;
2047    return mod (result, M);
2048  }
2049
2050  int sizeF= size (F);
2051  int sizeG= size (G);
2052
2053  int fallBackToNaive= 50;
2054  if (sizeF < fallBackToNaive || sizeG < fallBackToNaive)
2055    return mod (F*G, M);
2056
2057#ifdef HAVE_FLINT
2058  Variable alpha;
2059  if (getCharacteristic() == 0 && !hasFirstAlgVar (F, alpha) && ! hasFirstAlgVar (G, alpha))
2060    return mulMod2FLINTQ (F, G, M);
2061#endif
2062
2063  if (getCharacteristic() > 0 && CFFactory::gettype() != GaloisFieldDomain &&
2064      (((degF-degG) < 50 && degF > degG) || ((degG-degF) < 50 && degF <= degG)))
2065    return mulMod2NTLFq (F, G, M);
2066
2067  int m= (int) ceil (degree (M)/2.0);
2068  if (degF >= m || degG >= m)
2069  {
2070    CanonicalForm MLo= power (y, m);
2071    CanonicalForm MHi= power (y, degree (M) - m);
2072    CanonicalForm F0= mod (F, MLo);
2073    CanonicalForm F1= div (F, MLo);
2074    CanonicalForm G0= mod (G, MLo);
2075    CanonicalForm G1= div (G, MLo);
2076    CanonicalForm F0G1= mulMod2 (F0, G1, MHi);
2077    CanonicalForm F1G0= mulMod2 (F1, G0, MHi);
2078    CanonicalForm F0G0= mulMod2 (F0, G0, M);
2079    return F0G0 + MLo*(F0G1 + F1G0);
2080  }
2081  else
2082  {
2083    m= (int) ceil (tmax (degF, degG)/2.0);
2084    CanonicalForm yToM= power (y, m);
2085    CanonicalForm F0= mod (F, yToM);
2086    CanonicalForm F1= div (F, yToM);
2087    CanonicalForm G0= mod (G, yToM);
2088    CanonicalForm G1= div (G, yToM);
2089    CanonicalForm H00= mulMod2 (F0, G0, M);
2090    CanonicalForm H11= mulMod2 (F1, G1, M);
2091    CanonicalForm H01= mulMod2 (F0 + F1, G0 + G1, M);
2092    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
2093  }
2094  DEBOUTLN (cerr, "fatal end in mulMod2");
2095}
2096
2097CanonicalForm mulMod (const CanonicalForm& A, const CanonicalForm& B,
2098                      const CFList& MOD)
2099{
2100  if (A.isZero() || B.isZero())
2101    return 0;
2102
2103  if (MOD.length() == 1)
2104    return mulMod2 (A, B, MOD.getLast());
2105
2106  CanonicalForm M= MOD.getLast();
2107  CanonicalForm F= mod (A, M);
2108  CanonicalForm G= mod (B, M);
2109  if (F.inCoeffDomain() || G.inCoeffDomain())
2110    return F*G;
2111  Variable y= M.mvar();
2112  int degF= degree (F, y);
2113  int degG= degree (G, y);
2114
2115  if ((degF <= 1 && F.level() <= M.level()) &&
2116      (degG <= 1 && G.level() <= M.level()))
2117  {
2118    CFList buf= MOD;
2119    buf.removeLast();
2120    if (degF == 1 && degG == 1)
2121    {
2122      CanonicalForm F0= mod (F, y);
2123      CanonicalForm F1= div (F, y);
2124      CanonicalForm G0= mod (G, y);
2125      CanonicalForm G1= div (G, y);
2126      if (degree (M) > 2)
2127      {
2128        CanonicalForm H00= mulMod (F0, G0, buf);
2129        CanonicalForm H11= mulMod (F1, G1, buf);
2130        CanonicalForm H01= mulMod (F0 + F1, G0 + G1, buf);
2131        return H11*y*y + (H01 - H00 - H11)*y + H00;
2132      }
2133      else //here degree (M) == 2
2134      {
2135        buf.append (y);
2136        CanonicalForm F0G1= mulMod (F0, G1, buf);
2137        CanonicalForm F1G0= mulMod (F1, G0, buf);
2138        CanonicalForm F0G0= mulMod (F0, G0, MOD);
2139        CanonicalForm result= F0G0 + y*(F0G1 + F1G0);
2140        return result;
2141      }
2142    }
2143    else if (degF == 1 && degG == 0)
2144      return mulMod (div (F, y), G, buf)*y + mulMod (mod (F, y), G, buf);
2145    else if (degF == 0 && degG == 1)
2146      return mulMod (div (G, y), F, buf)*y + mulMod (mod (G, y), F, buf);
2147    else
2148      return mulMod (F, G, buf);
2149  }
2150  int m= (int) ceil (degree (M)/2.0);
2151  if (degF >= m || degG >= m)
2152  {
2153    CanonicalForm MLo= power (y, m);
2154    CanonicalForm MHi= power (y, degree (M) - m);
2155    CanonicalForm F0= mod (F, MLo);
2156    CanonicalForm F1= div (F, MLo);
2157    CanonicalForm G0= mod (G, MLo);
2158    CanonicalForm G1= div (G, MLo);
2159    CFList buf= MOD;
2160    buf.removeLast();
2161    buf.append (MHi);
2162    CanonicalForm F0G1= mulMod (F0, G1, buf);
2163    CanonicalForm F1G0= mulMod (F1, G0, buf);
2164    CanonicalForm F0G0= mulMod (F0, G0, MOD);
2165    return F0G0 + MLo*(F0G1 + F1G0);
2166  }
2167  else
2168  {
2169    m= (int) ceil (tmax (degF, degG)/2.0);
2170    CanonicalForm yToM= power (y, m);
2171    CanonicalForm F0= mod (F, yToM);
2172    CanonicalForm F1= div (F, yToM);
2173    CanonicalForm G0= mod (G, yToM);
2174    CanonicalForm G1= div (G, yToM);
2175    CanonicalForm H00= mulMod (F0, G0, MOD);
2176    CanonicalForm H11= mulMod (F1, G1, MOD);
2177    CanonicalForm H01= mulMod (F0 + F1, G0 + G1, MOD);
2178    return H11*yToM*yToM + (H01 - H11 - H00)*yToM + H00;
2179  }
2180  DEBOUTLN (cerr, "fatal end in mulMod");
2181}
2182
2183CanonicalForm prodMod (const CFList& L, const CanonicalForm& M)
2184{
2185  if (L.isEmpty())
2186    return 1;
2187  int l= L.length();
2188  if (l == 1)
2189    return mod (L.getFirst(), M);
2190  else if (l == 2) {
2191    CanonicalForm result= mulMod2 (L.getFirst(), L.getLast(), M);
2192    return result;
2193  }
2194  else
2195  {
2196    l /= 2;
2197    CFList tmp1, tmp2;
2198    CFListIterator i= L;
2199    CanonicalForm buf1, buf2;
2200    for (int j= 1; j <= l; j++, i++)
2201      tmp1.append (i.getItem());
2202    tmp2= Difference (L, tmp1);
2203    buf1= prodMod (tmp1, M);
2204    buf2= prodMod (tmp2, M);
2205    CanonicalForm result= mulMod2 (buf1, buf2, M);
2206    return result;
2207  }
2208}
2209
2210CanonicalForm prodMod (const CFList& L, const CFList& M)
2211{
2212  if (L.isEmpty())
2213    return 1;
2214  else if (L.length() == 1)
2215    return L.getFirst();
2216  else if (L.length() == 2)
2217    return mulMod (L.getFirst(), L.getLast(), M);
2218  else
2219  {
2220    int l= L.length()/2;
2221    CFListIterator i= L;
2222    CFList tmp1, tmp2;
2223    CanonicalForm buf1, buf2;
2224    for (int j= 1; j <= l; j++, i++)
2225      tmp1.append (i.getItem());
2226    tmp2= Difference (L, tmp1);
2227    buf1= prodMod (tmp1, M);
2228    buf2= prodMod (tmp2, M);
2229    return mulMod (buf1, buf2, M);
2230  }
2231}
2232
2233
2234CanonicalForm reverse (const CanonicalForm& F, int d)
2235{
2236  if (d == 0)
2237    return F;
2238  CanonicalForm A= F;
2239  Variable y= Variable (2);
2240  Variable x= Variable (1);
2241  if (degree (A, x) > 0)
2242  {
2243    A= swapvar (A, x, y);
2244    CanonicalForm result= 0;
2245    CFIterator i= A;
2246    while (d - i.exp() < 0)
2247      i++;
2248
2249    for (; i.hasTerms() && (d - i.exp() >= 0); i++)
2250      result += swapvar (i.coeff(),x,y)*power (x, d - i.exp());
2251    return result;
2252  }
2253  else
2254    return A*power (x, d);
2255}
2256
2257CanonicalForm
2258newtonInverse (const CanonicalForm& F, const int n, const CanonicalForm& M)
2259{
2260  int l= ilog2(n);
2261
2262  CanonicalForm g= mod (F, M)[0] [0];
2263
2264  ASSERT (!g.isZero(), "expected a unit");
2265
2266  Variable alpha;
2267
2268  if (!g.isOne())
2269    g = 1/g;
2270  Variable x= Variable (1);
2271  CanonicalForm result;
2272  int exp= 0;
2273  if (n & 1)
2274  {
2275    result= g;
2276    exp= 1;
2277  }
2278  CanonicalForm h;
2279
2280  for (int i= 1; i <= l; i++)
2281  {
2282    h= mulMod2 (g, mod (F, power (x, (1 << i))), M);
2283    h= mod (h, power (x, (1 << i)) - 1);
2284    h= div (h, power (x, (1 << (i - 1))));
2285    h= mod (h, M);
2286    g -= power (x, (1 << (i - 1)))*
2287         mod (mulMod2 (g, h, M), power (x, (1 << (i - 1))));
2288
2289    if (n & (1 << i))
2290    {
2291      if (exp)
2292      {
2293        h= mulMod2 (result, mod (F, power (x, exp + (1 << i))), M);
2294        h= mod (h, power (x, exp + (1 << i)) - 1);
2295        h= div (h, power (x, exp));
2296        h= mod (h, M);
2297        result -= power(x, exp)*mod (mulMod2 (g, h, M),
2298                                       power (x, (1 << i)));
2299        exp += (1 << i);
2300      }
2301      else
2302      {
2303        exp= (1 << i);
2304        result= g;
2305      }
2306    }
2307  }
2308
2309  return result;
2310}
2311
2312CanonicalForm
2313newtonDiv (const CanonicalForm& F, const CanonicalForm& G, const CanonicalForm&
2314           M)
2315{
2316  ASSERT (getCharacteristic() > 0, "positive characteristic expected");
2317  ASSERT (CFFactory::gettype() != GaloisFieldDomain, "no GF expected");
2318
2319  CanonicalForm A= mod (F, M);
2320  CanonicalForm B= mod (G, M);
2321
2322  Variable x= Variable (1);
2323  int degA= degree (A, x);
2324  int degB= degree (B, x);
2325  int m= degA - degB;
2326  if (m < 0)
2327    return 0;
2328
2329  Variable v;
2330  CanonicalForm Q;
2331  if (degB < 1 || CFFactory::gettype() == GaloisFieldDomain)
2332  {
2333    CanonicalForm R;
2334    divrem2 (A, B, Q, R, M);
2335  }
2336  else
2337  {
2338    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
2339    {
2340      CanonicalForm R= reverse (A, degA);
2341      CanonicalForm revB= reverse (B, degB);
2342      revB= newtonInverse (revB, m + 1, M);
2343      Q= mulMod2 (R, revB, M);
2344      Q= mod (Q, power (x, m + 1));
2345      Q= reverse (Q, m);
2346    }
2347    else
2348    {
2349      zz_pX mipo= convertFacCF2NTLzzpX (M);
2350      Variable y= Variable (2);
2351      zz_pEX NTLA, NTLB;
2352      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
2353      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
2354      div (NTLA, NTLA, NTLB);
2355      Q= convertNTLzz_pEX2CF (NTLA, x, y);
2356    }
2357  }
2358
2359  return Q;
2360}
2361
2362void
2363newtonDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2364              CanonicalForm& R, const CanonicalForm& M)
2365{
2366  CanonicalForm A= mod (F, M);
2367  CanonicalForm B= mod (G, M);
2368  Variable x= Variable (1);
2369  int degA= degree (A, x);
2370  int degB= degree (B, x);
2371  int m= degA - degB;
2372
2373  if (m < 0)
2374  {
2375    R= A;
2376    Q= 0;
2377    return;
2378  }
2379
2380  Variable v;
2381  if (degB <= 1 || CFFactory::gettype() == GaloisFieldDomain)
2382  {
2383     divrem2 (A, B, Q, R, M);
2384  }
2385  else
2386  {
2387    if (hasFirstAlgVar (A, v) || hasFirstAlgVar (B, v))
2388    {
2389      R= reverse (A, degA);
2390
2391      CanonicalForm revB= reverse (B, degB);
2392      revB= newtonInverse (revB, m + 1, M);
2393      Q= mulMod2 (R, revB, M);
2394
2395      Q= mod (Q, power (x, m + 1));
2396      Q= reverse (Q, m);
2397
2398      R= A - mulMod2 (Q, B, M);
2399    }
2400    else
2401    {
2402      zz_pX mipo= convertFacCF2NTLzzpX (M);
2403      Variable y= Variable (2);
2404      zz_pEX NTLA, NTLB;
2405      NTLA= convertFacCF2NTLzz_pEX (swapvar (A, x, y), mipo);
2406      NTLB= convertFacCF2NTLzz_pEX (swapvar (B, x, y), mipo);
2407      zz_pEX NTLQ, NTLR;
2408      DivRem (NTLQ, NTLR, NTLA, NTLB);
2409      Q= convertNTLzz_pEX2CF (NTLQ, x, y);
2410      R= convertNTLzz_pEX2CF (NTLR, x, y);
2411    }
2412  }
2413}
2414
2415static inline
2416CFList split (const CanonicalForm& F, const int m, const Variable& x)
2417{
2418  CanonicalForm A= F;
2419  CanonicalForm buf= 0;
2420  bool swap= false;
2421  if (degree (A, x) <= 0)
2422    return CFList(A);
2423  else if (x.level() != A.level())
2424  {
2425    swap= true;
2426    A= swapvar (A, x, A.mvar());
2427  }
2428
2429  int j= (int) floor ((double) degree (A)/ m);
2430  CFList result;
2431  CFIterator i= A;
2432  for (; j >= 0; j--)
2433  {
2434    while (i.hasTerms() && i.exp() - j*m >= 0)
2435    {
2436      if (swap)
2437        buf += i.coeff()*power (A.mvar(), i.exp() - j*m);
2438      else
2439        buf += i.coeff()*power (x, i.exp() - j*m);
2440      i++;
2441    }
2442    if (swap)
2443      result.append (swapvar (buf, x, F.mvar()));
2444    else
2445      result.append (buf);
2446    buf= 0;
2447  }
2448  return result;
2449}
2450
2451static inline
2452void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2453               CanonicalForm& R, const CFList& M);
2454
2455static inline
2456void divrem21 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2457               CanonicalForm& R, const CFList& M)
2458{
2459  CanonicalForm A= mod (F, M);
2460  CanonicalForm B= mod (G, M);
2461  Variable x= Variable (1);
2462  int degB= degree (B, x);
2463  int degA= degree (A, x);
2464  if (degA < degB)
2465  {
2466    Q= 0;
2467    R= A;
2468    return;
2469  }
2470  ASSERT (2*degB > degA, "expected degree (F, 1) < 2*degree (G, 1)");
2471  if (degB < 1)
2472  {
2473    divrem (A, B, Q, R);
2474    Q= mod (Q, M);
2475    R= mod (R, M);
2476    return;
2477  }
2478
2479  int m= (int) ceil ((double) (degB + 1)/2.0) + 1;
2480  CFList splitA= split (A, m, x);
2481  if (splitA.length() == 3)
2482    splitA.insert (0);
2483  if (splitA.length() == 2)
2484  {
2485    splitA.insert (0);
2486    splitA.insert (0);
2487  }
2488  if (splitA.length() == 1)
2489  {
2490    splitA.insert (0);
2491    splitA.insert (0);
2492    splitA.insert (0);
2493  }
2494
2495  CanonicalForm xToM= power (x, m);
2496
2497  CFListIterator i= splitA;
2498  CanonicalForm H= i.getItem();
2499  i++;
2500  H *= xToM;
2501  H += i.getItem();
2502  i++;
2503  H *= xToM;
2504  H += i.getItem();
2505  i++;
2506
2507  divrem32 (H, B, Q, R, M);
2508
2509  CFList splitR= split (R, m, x);
2510  if (splitR.length() == 1)
2511    splitR.insert (0);
2512
2513  H= splitR.getFirst();
2514  H *= xToM;
2515  H += splitR.getLast();
2516  H *= xToM;
2517  H += i.getItem();
2518
2519  CanonicalForm bufQ;
2520  divrem32 (H, B, bufQ, R, M);
2521
2522  Q *= xToM;
2523  Q += bufQ;
2524  return;
2525}
2526
2527static inline
2528void divrem32 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2529               CanonicalForm& R, const CFList& M)
2530{
2531  CanonicalForm A= mod (F, M);
2532  CanonicalForm B= mod (G, M);
2533  Variable x= Variable (1);
2534  int degB= degree (B, x);
2535  int degA= degree (A, x);
2536  if (degA < degB)
2537  {
2538    Q= 0;
2539    R= A;
2540    return;
2541  }
2542  ASSERT (3*(degB/2) > degA, "expected degree (F, 1) < 3*(degree (G, 1)/2)");
2543  if (degB < 1)
2544  {
2545    divrem (A, B, Q, R);
2546    Q= mod (Q, M);
2547    R= mod (R, M);
2548    return;
2549  }
2550  int m= (int) ceil ((double) (degB + 1)/ 2.0);
2551
2552  CFList splitA= split (A, m, x);
2553  CFList splitB= split (B, m, x);
2554
2555  if (splitA.length() == 2)
2556  {
2557    splitA.insert (0);
2558  }
2559  if (splitA.length() == 1)
2560  {
2561    splitA.insert (0);
2562    splitA.insert (0);
2563  }
2564  CanonicalForm xToM= power (x, m);
2565
2566  CanonicalForm H;
2567  CFListIterator i= splitA;
2568  i++;
2569
2570  if (degree (splitA.getFirst(), x) < degree (splitB.getFirst(), x))
2571  {
2572    H= splitA.getFirst()*xToM + i.getItem();
2573    divrem21 (H, splitB.getFirst(), Q, R, M);
2574  }
2575  else
2576  {
2577    R= splitA.getFirst()*xToM + i.getItem() + splitB.getFirst() -
2578       splitB.getFirst()*xToM;
2579    Q= xToM - 1;
2580  }
2581
2582  H= mulMod (Q, splitB.getLast(), M);
2583
2584  R= R*xToM + splitA.getLast() - H;
2585
2586  while (degree (R, x) >= degB)
2587  {
2588    xToM= power (x, degree (R, x) - degB);
2589    Q += LC (R, x)*xToM;
2590    R -= mulMod (LC (R, x), B, M)*xToM;
2591    Q= mod (Q, M);
2592    R= mod (R, M);
2593  }
2594
2595  return;
2596}
2597
2598void divrem2 (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2599              CanonicalForm& R, const CanonicalForm& M)
2600{
2601  CanonicalForm A= mod (F, M);
2602  CanonicalForm B= mod (G, M);
2603
2604  if (B.inCoeffDomain())
2605  {
2606    divrem (A, B, Q, R);
2607    return;
2608  }
2609  if (A.inCoeffDomain() && !B.inCoeffDomain())
2610  {
2611    Q= 0;
2612    R= A;
2613    return;
2614  }
2615
2616  if (B.level() < A.level())
2617  {
2618    divrem (A, B, Q, R);
2619    return;
2620  }
2621  if (A.level() > B.level())
2622  {
2623    R= A;
2624    Q= 0;
2625    return;
2626  }
2627  if (B.level() == 1 && B.isUnivariate())
2628  {
2629    divrem (A, B, Q, R);
2630    return;
2631  }
2632  if (!(B.level() == 1 && B.isUnivariate()) && (A.level() == 1 && A.isUnivariate()))
2633  {
2634    Q= 0;
2635    R= A;
2636    return;
2637  }
2638
2639  Variable x= Variable (1);
2640  int degB= degree (B, x);
2641  if (degB > degree (A, x))
2642  {
2643    Q= 0;
2644    R= A;
2645    return;
2646  }
2647
2648  CFList splitA= split (A, degB, x);
2649
2650  CanonicalForm xToDegB= power (x, degB);
2651  CanonicalForm H, bufQ;
2652  Q= 0;
2653  CFListIterator i= splitA;
2654  H= i.getItem()*xToDegB;
2655  i++;
2656  H += i.getItem();
2657  CFList buf;
2658  while (i.hasItem())
2659  {
2660    buf= CFList (M);
2661    divrem21 (H, B, bufQ, R, buf);
2662    i++;
2663    if (i.hasItem())
2664      H= R*xToDegB + i.getItem();
2665    Q *= xToDegB;
2666    Q += bufQ;
2667  }
2668  return;
2669}
2670
2671void divrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
2672             CanonicalForm& R, const CFList& MOD)
2673{
2674  CanonicalForm A= mod (F, MOD);
2675  CanonicalForm B= mod (G, MOD);
2676  Variable x= Variable (1);
2677  int degB= degree (B, x);
2678  if (degB > degree (A, x))
2679  {
2680    Q= 0;
2681    R= A;
2682    return;
2683  }
2684
2685  if (degB <= 0)
2686  {
2687    divrem (A, B, Q, R);
2688    Q= mod (Q, MOD);
2689    R= mod (R, MOD);
2690    return;
2691  }
2692  CFList splitA= split (A, degB, x);
2693
2694  CanonicalForm xToDegB= power (x, degB);
2695  CanonicalForm H, bufQ;
2696  Q= 0;
2697  CFListIterator i= splitA;
2698  H= i.getItem()*xToDegB;
2699  i++;
2700  H += i.getItem();
2701  while (i.hasItem())
2702  {
2703    divrem21 (H, B, bufQ, R, MOD);
2704    i++;
2705    if (i.hasItem())
2706      H= R*xToDegB + i.getItem();
2707    Q *= xToDegB;
2708    Q += bufQ;
2709  }
2710  return;
2711}
2712
2713void sortList (CFList& list, const Variable& x)
2714{
2715  int l= 1;
2716  int k= 1;
2717  CanonicalForm buf;
2718  CFListIterator m;
2719  for (CFListIterator i= list; l <= list.length(); i++, l++)
2720  {
2721    for (CFListIterator j= list; k <= list.length() - l; k++)
2722    {
2723      m= j;
2724      m++;
2725      if (degree (j.getItem(), x) > degree (m.getItem(), x))
2726      {
2727        buf= m.getItem();
2728        m.getItem()= j.getItem();
2729        j.getItem()= buf;
2730        j++;
2731        j.getItem()= m.getItem();
2732      }
2733      else
2734        j++;
2735    }
2736    k= 1;
2737  }
2738}
2739
2740static inline
2741CFList diophantine (const CanonicalForm& F, const CFList& factors)
2742{
2743  if (getCharacteristic() == 0)
2744  {
2745    Variable v;
2746    bool hasAlgVar= hasFirstAlgVar (F, v);
2747    for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++)
2748      hasAlgVar= hasFirstAlgVar (i.getItem(), v);
2749    if (hasAlgVar)
2750    {
2751      CFList result= modularDiophant (F, factors, getMipo (v));
2752      return result;
2753    }
2754  }
2755
2756  CanonicalForm buf1, buf2, buf3, S, T;
2757  CFListIterator i= factors;
2758  CFList result;
2759  if (i.hasItem())
2760    i++;
2761  buf1= F/factors.getFirst();
2762  buf2= divNTL (F, i.getItem());
2763  buf3= extgcd (buf1, buf2, S, T);
2764  result.append (S);
2765  result.append (T);
2766  if (i.hasItem())
2767    i++;
2768  for (; i.hasItem(); i++)
2769  {
2770    buf1= divNTL (F, i.getItem());
2771    buf3= extgcd (buf3, buf1, S, T);
2772    CFListIterator k= factors;
2773    for (CFListIterator j= result; j.hasItem(); j++, k++)
2774    {
2775      j.getItem()= mulNTL (j.getItem(), S);
2776      j.getItem()= modNTL (j.getItem(), k.getItem());
2777    }
2778    result.append (T);
2779  }
2780  return result;
2781}
2782
2783void
2784henselStep12 (const CanonicalForm& F, const CFList& factors,
2785              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
2786              CFArray& Pi, int j)
2787{
2788  CanonicalForm E;
2789  CanonicalForm xToJ= power (F.mvar(), j);
2790  Variable x= F.mvar();
2791  // compute the error
2792  if (j == 1)
2793    E= F[j];
2794  else
2795  {
2796    if (degree (Pi [factors.length() - 2], x) > 0)
2797      E= F[j] - Pi [factors.length() - 2] [j];
2798    else
2799      E= F[j];
2800  }
2801
2802  CFArray buf= CFArray (diophant.length());
2803  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
2804  int k= 0;
2805  CanonicalForm remainder;
2806  // actual lifting
2807  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
2808  {
2809    if (degree (bufFactors[k], x) > 0)
2810    {
2811      if (k > 0)
2812        remainder= modNTL (E, bufFactors[k] [0]);
2813      else
2814        remainder= E;
2815    }
2816    else
2817      remainder= modNTL (E, bufFactors[k]);
2818
2819    buf[k]= mulNTL (i.getItem(), remainder);
2820    if (degree (bufFactors[k], x) > 0)
2821      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
2822    else
2823      buf[k]= modNTL (buf[k], bufFactors[k]);
2824  }
2825  for (k= 1; k < factors.length(); k++)
2826    bufFactors[k] += xToJ*buf[k];
2827
2828  // update Pi [0]
2829  int degBuf0= degree (bufFactors[0], x);
2830  int degBuf1= degree (bufFactors[1], x);
2831  if (degBuf0 > 0 && degBuf1 > 0)
2832    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
2833  CanonicalForm uIZeroJ;
2834  if (j == 1)
2835  {
2836    if (degBuf0 > 0 && degBuf1 > 0)
2837      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
2838                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
2839    else if (degBuf0 > 0)
2840      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
2841    else if (degBuf1 > 0)
2842      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
2843    else
2844      uIZeroJ= 0;
2845    Pi [0] += xToJ*uIZeroJ;
2846  }
2847  else
2848  {
2849    if (degBuf0 > 0 && degBuf1 > 0)
2850      uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
2851                  (bufFactors[1] [0] + buf[1])) - M(1, 1) - M(j + 1, 1);
2852    else if (degBuf0 > 0)
2853      uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1]);
2854    else if (degBuf1 > 0)
2855      uIZeroJ= mulNTL (bufFactors[0], buf[1]);
2856    else
2857      uIZeroJ= 0;
2858    Pi [0] += xToJ*uIZeroJ;
2859  }
2860  CFArray tmp= CFArray (factors.length() - 1);
2861  for (k= 0; k < factors.length() - 1; k++)
2862    tmp[k]= 0;
2863  CFIterator one, two;
2864  one= bufFactors [0];
2865  two= bufFactors [1];
2866  if (degBuf0 > 0 && degBuf1 > 0)
2867  {
2868    for (k= 1; k <= (int) ceil (j/2.0); k++)
2869    {
2870      if (k != j - k + 1)
2871      {
2872        if ((one.hasTerms() && one.exp() == j - k + 1) && (two.hasTerms() && two.exp() == j - k + 1))
2873        {
2874          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
2875                     two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
2876          one++;
2877          two++;
2878        }
2879        else if (one.hasTerms() && one.exp() == j - k + 1)
2880        {
2881          tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) -
2882                     M (k + 1, 1);
2883          one++;
2884        }
2885        else if (two.hasTerms() && two.exp() == j - k + 1)
2886        {
2887          tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) -
2888                    M (k + 1, 1);
2889          two++;
2890        }
2891      }
2892      else
2893      {
2894        tmp[0] += M (k + 1, 1);
2895      }
2896    }
2897  }
2898  Pi [0] += tmp[0]*xToJ*F.mvar();
2899
2900  // update Pi [l]
2901  int degPi, degBuf;
2902  for (int l= 1; l < factors.length() - 1; l++)
2903  {
2904    degPi= degree (Pi [l - 1], x);
2905    degBuf= degree (bufFactors[l + 1], x);
2906    if (degPi > 0 && degBuf > 0)
2907      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
2908    if (j == 1)
2909    {
2910      if (degPi > 0 && degBuf > 0)
2911        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
2912                  bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
2913                  M (1, l + 1));
2914      else if (degPi > 0)
2915        Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
2916      else if (degBuf > 0)
2917        Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
2918    }
2919    else
2920    {
2921      if (degPi > 0 && degBuf > 0)
2922      {
2923        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
2924        uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
2925      }
2926      else if (degPi > 0)
2927        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
2928      else if (degBuf > 0)
2929      {
2930        uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
2931        uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
2932      }
2933      Pi[l] += xToJ*uIZeroJ;
2934    }
2935    one= bufFactors [l + 1];
2936    two= Pi [l - 1];
2937    if (two.hasTerms() && two.exp() == j + 1)
2938    {
2939      if (degBuf > 0 && degPi > 0)
2940      {
2941          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
2942          two++;
2943      }
2944      else if (degPi > 0)
2945      {
2946          tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
2947          two++;
2948      }
2949    }
2950    if (degBuf > 0 && degPi > 0)
2951    {
2952      for (k= 1; k <= (int) ceil (j/2.0); k++)
2953      {
2954        if (k != j - k + 1)
2955        {
2956          if ((one.hasTerms() && one.exp() == j - k + 1) &&
2957              (two.hasTerms() && two.exp() == j - k + 1))
2958          {
2959            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
2960                       two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
2961            one++;
2962            two++;
2963          }
2964          else if (one.hasTerms() && one.exp() == j - k + 1)
2965          {
2966            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) -
2967                       M (k + 1, l + 1);
2968            one++;
2969          }
2970          else if (two.hasTerms() && two.exp() == j - k + 1)
2971          {
2972            tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) -
2973                      M (k + 1, l + 1);
2974            two++;
2975          }
2976        }
2977        else
2978          tmp[l] += M (k + 1, l + 1);
2979      }
2980    }
2981    Pi[l] += tmp[l]*xToJ*F.mvar();
2982  }
2983  return;
2984}
2985
2986void
2987henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
2988              CFList& diophant, CFMatrix& M, bool sort)
2989{
2990  if (sort)
2991    sortList (factors, Variable (1));
2992  Pi= CFArray (factors.length() - 1);
2993  CFListIterator j= factors;
2994  diophant= diophantine (F[0], factors);
2995  DEBOUTLN (cerr, "diophant= " << diophant);
2996  j++;
2997  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()));
2998  M (1, 1)= Pi [0];
2999  int i= 1;
3000  if (j.hasItem())
3001    j++;
3002  for (; j.hasItem(); j++, i++)
3003  {
3004    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
3005    M (1, i + 1)= Pi [i];
3006  }
3007  CFArray bufFactors= CFArray (factors.length());
3008  i= 0;
3009  for (CFListIterator k= factors; k.hasItem(); i++, k++)
3010  {
3011    if (i == 0)
3012      bufFactors[i]= mod (k.getItem(), F.mvar());
3013    else
3014      bufFactors[i]= k.getItem();
3015  }
3016  for (i= 1; i < l; i++)
3017    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
3018
3019  CFListIterator k= factors;
3020  for (i= 0; i < factors.length (); i++, k++)
3021    k.getItem()= bufFactors[i];
3022  factors.removeFirst();
3023  return;
3024}
3025
3026void
3027henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
3028                    end, CFArray& Pi, const CFList& diophant, CFMatrix& M)
3029{
3030  CFArray bufFactors= CFArray (factors.length());
3031  int i= 0;
3032  CanonicalForm xToStart= power (F.mvar(), start);
3033  for (CFListIterator k= factors; k.hasItem(); k++, i++)
3034  {
3035    if (i == 0)
3036      bufFactors[i]= mod (k.getItem(), xToStart);
3037    else
3038      bufFactors[i]= k.getItem();
3039  }
3040  for (i= start; i < end; i++)
3041    henselStep12 (F, factors, bufFactors, diophant, M, Pi, i);
3042
3043  CFListIterator k= factors;
3044  for (i= 0; i < factors.length(); k++, i++)
3045    k.getItem()= bufFactors [i];
3046  factors.removeFirst();
3047  return;
3048}
3049
3050static inline
3051CFList
3052biDiophantine (const CanonicalForm& F, const CFList& factors, const int d)
3053{
3054  Variable y= F.mvar();
3055  CFList result;
3056  if (y.level() == 1)
3057  {
3058    result= diophantine (F, factors);
3059    return result;
3060  }
3061  else
3062  {
3063    CFList buf= factors;
3064    for (CFListIterator i= buf; i.hasItem(); i++)
3065      i.getItem()= mod (i.getItem(), y);
3066    CanonicalForm A= mod (F, y);
3067    int bufD= 1;
3068    CFList recResult= biDiophantine (A, buf, bufD);
3069    CanonicalForm e= 1;
3070    CFList p;
3071    CFArray bufFactors= CFArray (factors.length());
3072    CanonicalForm yToD= power (y, d);
3073    int k= 0;
3074    for (CFListIterator i= factors; i.hasItem(); i++, k++)
3075    {
3076      bufFactors [k]= i.getItem();
3077    }
3078    CanonicalForm b, quot;
3079    for (k= 0; k < factors.length(); k++) //TODO compute b's faster
3080    {
3081      b= 1;
3082      if (fdivides (bufFactors[k], F, quot))
3083        b= quot;
3084      else
3085      {
3086        for (int l= 0; l < factors.length(); l++)
3087        {
3088          if (l == k)
3089            continue;
3090          else
3091          {
3092            b= mulMod2 (b, bufFactors[l], yToD);
3093          }
3094        }
3095      }
3096      p.append (b);
3097    }
3098
3099    CFListIterator j= p;
3100    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
3101      e -= i.getItem()*j.getItem();
3102
3103    if (e.isZero())
3104      return recResult;
3105    CanonicalForm coeffE;
3106    CFList s;
3107    result= recResult;
3108    CanonicalForm g;
3109    for (int i= 1; i < d; i++)
3110    {
3111      if (degree (e, y) > 0)
3112        coeffE= e[i];
3113      else
3114        coeffE= 0;
3115      if (!coeffE.isZero())
3116      {
3117        CFListIterator k= result;
3118        CFListIterator l= p;
3119        int ii= 0;
3120        j= recResult;
3121        for (; j.hasItem(); j++, k++, l++, ii++)
3122        {
3123          g= coeffE*j.getItem();
3124          if (degree (bufFactors[ii], y) <= 0)
3125            g= mod (g, bufFactors[ii]);
3126          else
3127            g= mod (g, bufFactors[ii][0]);
3128          k.getItem() += g*power (y, i);
3129          e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
3130          DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
3131                    mod (e, power (y, i + 1)));
3132        }
3133      }
3134      if (e.isZero())
3135        break;
3136    }
3137
3138    DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
3139
3140#ifdef DEBUGOUTPUT
3141    CanonicalForm test= 0;
3142    j= p;
3143    for (CFListIterator i= result; i.hasItem(); i++, j++)
3144      test += mod (i.getItem()*j.getItem(), power (y, d));
3145    DEBOUTLN (cerr, "test= " << test);
3146#endif
3147    return result;
3148  }
3149}
3150
3151static inline
3152CFList
3153multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
3154                     const CFList& recResult, const CFList& M, const int d)
3155{
3156  Variable y= F.mvar();
3157  CFList result;
3158  CFListIterator i;
3159  CanonicalForm e= 1;
3160  CFListIterator j= factors;
3161  CFList p;
3162  CFArray bufFactors= CFArray (factors.length());
3163  CanonicalForm yToD= power (y, d);
3164  int k= 0;
3165  for (CFListIterator i= factors; i.hasItem(); i++, k++)
3166    bufFactors [k]= i.getItem();
3167  CanonicalForm b, quot;
3168  CFList buf= M;
3169  buf.removeLast();
3170  buf.append (yToD);
3171  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
3172  {
3173    b= 1;
3174    if (fdivides (bufFactors[k], F, quot))
3175      b= quot;
3176    else
3177    {
3178      for (int l= 0; l < factors.length(); l++)
3179      {
3180        if (l == k)
3181          continue;
3182        else
3183        {
3184          b= mulMod (b, bufFactors[l], buf);
3185        }
3186      }
3187    }
3188    p.append (b);
3189  }
3190  j= p;
3191  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
3192    e -= mulMod (i.getItem(), j.getItem(), M);
3193
3194  if (e.isZero())
3195    return recResult;
3196  CanonicalForm coeffE;
3197  CFList s;
3198  result= recResult;
3199  CanonicalForm g;
3200  for (int i= 1; i < d; i++)
3201  {
3202    if (degree (e, y) > 0)
3203      coeffE= e[i];
3204    else
3205      coeffE= 0;
3206    if (!coeffE.isZero())
3207    {
3208      CFListIterator k= result;
3209      CFListIterator l= p;
3210      j= recResult;
3211      int ii= 0;
3212      CanonicalForm dummy;
3213      for (; j.hasItem(); j++, k++, l++, ii++)
3214      {
3215        g= mulMod (coeffE, j.getItem(), M);
3216        if (degree (bufFactors[ii], y) <= 0)
3217          divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
3218                  g, M);
3219        else
3220          divrem (g, bufFactors[ii][0], dummy, g, M);
3221        k.getItem() += g*power (y, i);
3222        e -= mulMod (g*power (y, i), l.getItem(), M);
3223      }
3224    }
3225
3226    if (e.isZero())
3227      break;
3228  }
3229
3230#ifdef DEBUGOUTPUT
3231    CanonicalForm test= 0;
3232    j= p;
3233    for (CFListIterator i= result; i.hasItem(); i++, j++)
3234      test += mod (i.getItem()*j.getItem(), power (y, d));
3235    DEBOUTLN (cerr, "test= " << test);
3236#endif
3237  return result;
3238}
3239
3240static inline
3241void
3242henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
3243            const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
3244            const CFList& MOD)
3245{
3246  CanonicalForm E;
3247  CanonicalForm xToJ= power (F.mvar(), j);
3248  Variable x= F.mvar();
3249  // compute the error
3250  if (j == 1)
3251  {
3252    E= F[j];
3253#ifdef DEBUGOUTPUT
3254    CanonicalForm test= 1;
3255    for (int i= 0; i < factors.length(); i++)
3256    {
3257      if (i == 0)
3258        test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
3259      else
3260        test= mulMod (test, bufFactors[i], MOD);
3261    }
3262    CanonicalForm test2= mod (F-test, xToJ);
3263
3264    test2= mod (test2, MOD);
3265    DEBOUTLN (cerr, "test= " << test2);
3266#endif
3267  }
3268  else
3269  {
3270#ifdef DEBUGOUTPUT
3271    CanonicalForm test= 1;
3272    for (int i= 0; i < factors.length(); i++)
3273    {
3274      if (i == 0)
3275        test *= mod (bufFactors [i], power (x, j));
3276      else
3277        test *= bufFactors[i];
3278    }
3279    test= mod (test, power (x, j));
3280    test= mod (test, MOD);
3281    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
3282    DEBOUTLN (cerr, "test= " << test2);
3283#endif
3284
3285    if (degree (Pi [factors.length() - 2], x) > 0)
3286      E= F[j] - Pi [factors.length() - 2] [j];
3287    else
3288      E= F[j];
3289  }
3290
3291  CFArray buf= CFArray (diophant.length());
3292  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
3293  int k= 0;
3294  // actual lifting
3295  CanonicalForm dummy, rest1;
3296  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
3297  {
3298    if (degree (bufFactors[k], x) > 0)
3299    {
3300      if (k > 0)
3301        divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
3302      else
3303        rest1= E;
3304    }
3305    else
3306      divrem (E, bufFactors[k], dummy, rest1, MOD);
3307
3308    buf[k]= mulMod (i.getItem(), rest1, MOD);
3309
3310    if (degree (bufFactors[k], x) > 0)
3311      divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
3312    else
3313      divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
3314  }
3315  for (k= 1; k < factors.length(); k++)
3316    bufFactors[k] += xToJ*buf[k];
3317
3318  // update Pi [0]
3319  int degBuf0= degree (bufFactors[0], x);
3320  int degBuf1= degree (bufFactors[1], x);
3321  if (degBuf0 > 0 && degBuf1 > 0)
3322    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
3323  CanonicalForm uIZeroJ;
3324  if (j == 1)
3325  {
3326    if (degBuf0 > 0 && degBuf1 > 0)
3327      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
3328                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
3329    else if (degBuf0 > 0)
3330      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
3331    else if (degBuf1 > 0)
3332      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
3333    else
3334      uIZeroJ= 0;
3335    Pi [0] += xToJ*uIZeroJ;
3336  }
3337  else
3338  {
3339    if (degBuf0 > 0 && degBuf1 > 0)
3340      uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
3341                  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
3342    else if (degBuf0 > 0)
3343      uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
3344    else if (degBuf1 > 0)
3345      uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
3346    else
3347      uIZeroJ= 0;
3348    Pi [0] += xToJ*uIZeroJ;
3349  }
3350
3351  CFArray tmp= CFArray (factors.length() - 1);
3352  for (k= 0; k < factors.length() - 1; k++)
3353    tmp[k]= 0;
3354  CFIterator one, two;
3355  one= bufFactors [0];
3356  two= bufFactors [1];
3357  if (degBuf0 > 0 && degBuf1 > 0)
3358  {
3359    for (k= 1; k <= (int) ceil (j/2.0); k++)
3360    {
3361      if (k != j - k + 1)
3362      {
3363        if ((one.hasTerms() && one.exp() == j - k + 1) &&
3364            (two.hasTerms() && two.exp() == j - k + 1))
3365        {
3366          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
3367                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
3368                    M (j - k + 2, 1);
3369          one++;
3370          two++;
3371        }
3372        else if (one.hasTerms() && one.exp() == j - k + 1)
3373        {
3374          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
3375                    bufFactors[1] [k], MOD) - M (k + 1, 1);
3376          one++;
3377        }
3378        else if (two.hasTerms() && two.exp() == j - k + 1)
3379        {
3380          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
3381                    two.coeff()), MOD) - M (k + 1, 1);
3382          two++;
3383        }
3384      }
3385      else
3386      {
3387        tmp[0] += M (k + 1, 1);
3388      }
3389    }
3390  }
3391  Pi [0] += tmp[0]*xToJ*F.mvar();
3392
3393  // update Pi [l]
3394  int degPi, degBuf;
3395  for (int l= 1; l < factors.length() - 1; l++)
3396  {
3397    degPi= degree (Pi [l - 1], x);
3398    degBuf= degree (bufFactors[l + 1], x);
3399    if (degPi > 0 && degBuf > 0)
3400      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
3401    if (j == 1)
3402    {
3403      if (degPi > 0 && degBuf > 0)
3404        Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
3405                  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
3406                  M (1, l + 1));
3407      else if (degPi > 0)
3408        Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
3409      else if (degBuf > 0)
3410        Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
3411    }
3412    else
3413    {
3414      if (degPi > 0 && degBuf > 0)
3415      {
3416        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
3417        uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
3418      }
3419      else if (degPi > 0)
3420        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
3421      else if (degBuf > 0)
3422      {
3423        uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
3424        uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
3425      }
3426      Pi[l] += xToJ*uIZeroJ;
3427    }
3428    one= bufFactors [l + 1];
3429    two= Pi [l - 1];
3430    if (two.hasTerms() && two.exp() == j + 1)
3431    {
3432      if (degBuf > 0 && degPi > 0)
3433      {
3434          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
3435          two++;
3436      }
3437      else if (degPi > 0)
3438      {
3439          tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
3440          two++;
3441      }
3442    }
3443    if (degBuf > 0 && degPi > 0)
3444    {
3445      for (k= 1; k <= (int) ceil (j/2.0); k++)
3446      {
3447        if (k != j - k + 1)
3448        {
3449          if ((one.hasTerms() && one.exp() == j - k + 1) &&
3450              (two.hasTerms() && two.exp() == j - k + 1))
3451          {
3452            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
3453                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
3454                      M (j - k + 2, l + 1);
3455            one++;
3456            two++;
3457          }
3458          else if (one.hasTerms() && one.exp() == j - k + 1)
3459          {
3460            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
3461                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
3462            one++;
3463          }
3464          else if (two.hasTerms() && two.exp() == j - k + 1)
3465          {
3466            tmp[l] += mulMod (bufFactors[l + 1] [k],
3467                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
3468            two++;
3469          }
3470        }
3471        else
3472          tmp[l] += M (k + 1, l + 1);
3473      }
3474    }
3475    Pi[l] += tmp[l]*xToJ*F.mvar();
3476  }
3477
3478  return;
3479}
3480
3481CFList
3482henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
3483              diophant, CFArray& Pi, CFMatrix& M)
3484{
3485  CFList buf= factors;
3486  int k= 0;
3487  int liftBoundBivar= l[k];
3488  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
3489  CFList MOD;
3490  MOD.append (power (Variable (2), liftBoundBivar));
3491  CFArray bufFactors= CFArray (factors.length());
3492  k= 0;
3493  CFListIterator j= eval;
3494  j++;
3495  buf.removeFirst();
3496  buf.insert (LC (j.getItem(), 1));
3497  for (CFListIterator i= buf; i.hasItem(); i++, k++)
3498    bufFactors[k]= i.getItem();
3499  Pi= CFArray (factors.length() - 1);
3500  CFListIterator i= buf;
3501  i++;
3502  Variable y= j.getItem().mvar();
3503  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
3504  M (1, 1)= Pi [0];
3505  k= 1;
3506  if (i.hasItem())
3507    i++;
3508  for (; i.hasItem(); i++, k++)
3509  {
3510    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
3511    M (1, k + 1)= Pi [k];
3512  }
3513
3514  for (int d= 1; d < l[1]; d++)
3515    henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
3516  CFList result;
3517  for (k= 1; k < factors.length(); k++)
3518    result.append (bufFactors[k]);
3519  return result;
3520}
3521
3522void
3523henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
3524                  CFArray& Pi, const CFList& diophant, CFMatrix& M,
3525                  const CFList& MOD)
3526{
3527  CFArray bufFactors= CFArray (factors.length());
3528  int i= 0;
3529  CanonicalForm xToStart= power (F.mvar(), start);
3530  for (CFListIterator k= factors; k.hasItem(); k++, i++)
3531  {
3532    if (i == 0)
3533      bufFactors[i]= mod (k.getItem(), xToStart);
3534    else
3535      bufFactors[i]= k.getItem();
3536  }
3537  for (i= start; i < end; i++)
3538    henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD);
3539
3540  CFListIterator k= factors;
3541  for (i= 0; i < factors.length(); k++, i++)
3542    k.getItem()= bufFactors [i];
3543  factors.removeFirst();
3544  return;
3545}
3546
3547CFList
3548henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
3549            diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
3550            lNew)
3551{
3552  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
3553  int k= 0;
3554  CFArray bufFactors= CFArray (factors.length());
3555  for (CFListIterator i= factors; i.hasItem(); i++, k++)
3556  {
3557    if (k == 0)
3558      bufFactors[k]= LC (F.getLast(), 1);
3559    else
3560      bufFactors[k]= i.getItem();
3561  }
3562  CFList buf= factors;
3563  buf.removeFirst();
3564  buf.insert (LC (F.getLast(), 1));
3565  CFListIterator i= buf;
3566  i++;
3567  Variable y= F.getLast().mvar();
3568  Variable x= F.getFirst().mvar();
3569  CanonicalForm xToLOld= power (x, lOld);
3570  Pi [0]= mod (Pi[0], xToLOld);
3571  M (1, 1)= Pi [0];
3572  k= 1;
3573  if (i.hasItem())
3574    i++;
3575  for (; i.hasItem(); i++, k++)
3576  {
3577    Pi [k]= mod (Pi [k], xToLOld);
3578    M (1, k + 1)= Pi [k];
3579  }
3580
3581  for (int d= 1; d < lNew; d++)
3582    henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
3583  CFList result;
3584  for (k= 1; k < factors.length(); k++)
3585    result.append (bufFactors[k]);
3586  return result;
3587}
3588
3589CFList
3590henselLift (const CFList& eval, const CFList& factors, const int* l, const int
3591            lLength, bool sort)
3592{
3593  CFList diophant;
3594  CFList buf= factors;
3595  buf.insert (LC (eval.getFirst(), 1));
3596  if (sort)
3597    sortList (buf, Variable (1));
3598  CFArray Pi;
3599  CFMatrix M= CFMatrix (l[1], factors.length());
3600  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
3601  if (eval.length() == 2)
3602    return result;
3603  CFList MOD;
3604  for (int i= 0; i < 2; i++)
3605    MOD.append (power (Variable (i + 2), l[i]));
3606  CFListIterator j= eval;
3607  j++;
3608  CFList bufEval;
3609  bufEval.append (j.getItem());
3610  j++;
3611
3612  for (int i= 2; i < lLength && j.hasItem(); i++, j++)
3613  {
3614    result.insert (LC (bufEval.getFirst(), 1));
3615    bufEval.append (j.getItem());
3616    M= CFMatrix (l[i], factors.length());
3617    result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
3618    MOD.append (power (Variable (i + 2), l[i]));
3619    bufEval.removeFirst();
3620  }
3621  return result;
3622}
3623
3624void
3625henselStep122 (const CanonicalForm& F, const CFList& factors,
3626              CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
3627              CFArray& Pi, int j, const CFArray& /*LCs*/)
3628{
3629  Variable x= F.mvar();
3630  CanonicalForm xToJ= power (x, j);
3631
3632  CanonicalForm E;
3633  // compute the error
3634  if (degree (Pi [factors.length() - 2], x) > 0)
3635    E= F[j] - Pi [factors.length() - 2] [j];
3636  else
3637    E= F[j];
3638
3639  CFArray buf= CFArray (diophant.length());
3640
3641  int k= 0;
3642  CanonicalForm remainder;
3643  // actual lifting
3644  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
3645  {
3646    if (degree (bufFactors[k], x) > 0)
3647      remainder= modNTL (E, bufFactors[k] [0]);
3648    else
3649      remainder= modNTL (E, bufFactors[k]);
3650    buf[k]= mulNTL (i.getItem(), remainder);
3651    if (degree (bufFactors[k], x) > 0)
3652      buf[k]= modNTL (buf[k], bufFactors[k] [0]);
3653    else
3654      buf[k]= modNTL (buf[k], bufFactors[k]);
3655  }
3656
3657  for (k= 0; k < factors.length(); k++)
3658    bufFactors[k] += xToJ*buf[k];
3659
3660  // update Pi [0]
3661  int degBuf0= degree (bufFactors[0], x);
3662  int degBuf1= degree (bufFactors[1], x);
3663  if (degBuf0 > 0 && degBuf1 > 0)
3664  {
3665    M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
3666    if (j + 2 <= M.rows())
3667      M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
3668  }
3669  CanonicalForm uIZeroJ;
3670  if (degBuf0 > 0 && degBuf1 > 0)
3671    uIZeroJ= mulNTL(bufFactors[0][0],buf[1])+mulNTL (bufFactors[1][0], buf[0]);
3672  else if (degBuf0 > 0)
3673    uIZeroJ= mulNTL (buf[0], bufFactors[1]);
3674  else if (degBuf1 > 0)
3675    uIZeroJ= mulNTL (bufFactors[0], buf [1]);
3676  else
3677    uIZeroJ= 0;
3678  Pi [0] += xToJ*uIZeroJ;
3679
3680  CFArray tmp= CFArray (factors.length() - 1);
3681  for (k= 0; k < factors.length() - 1; k++)
3682    tmp[k]= 0;
3683  CFIterator one, two;
3684  one= bufFactors [0];
3685  two= bufFactors [1];
3686  if (degBuf0 > 0 && degBuf1 > 0)
3687  {
3688    while (one.hasTerms() && one.exp() > j) one++;
3689    while (two.hasTerms() && two.exp() > j) two++;
3690    for (k= 1; k <= (int) ceil (j/2.0); k++)
3691    {
3692      if (one.hasTerms() && two.hasTerms())
3693      {
3694        if (k != j - k + 1)
3695        {
3696          if ((one.hasTerms() && one.exp() == j - k + 1) && +
3697              (two.hasTerms() && two.exp() == j - k + 1))
3698        {
3699          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
3700                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
3701          one++;
3702          two++;
3703        }
3704        else if (one.hasTerms() && one.exp() == j - k + 1)
3705        {
3706          tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
3707                      M (k + 1, 1);
3708          one++;
3709        }
3710        else if (two.hasTerms() && two.exp() == j - k + 1)
3711        {
3712          tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
3713                    M (k + 1, 1);
3714          two++;
3715        }
3716      }
3717      else
3718        tmp[0] += M (k + 1, 1);
3719      }
3720    }
3721  }
3722
3723  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
3724  {
3725    if (j + 2 <= M.rows())
3726      tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
3727                         (bufFactors [1] [j + 1] + bufFactors [1] [0]))
3728                         - M(1,1) - M (j + 2,1);
3729  }
3730  else if (degBuf0 >= j + 1)
3731  {
3732    if (degBuf1 > 0)
3733      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
3734    else
3735      tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]);
3736  }
3737  else if (degBuf1 >= j + 1)
3738  {
3739    if (degBuf0 > 0)
3740      tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
3741    else
3742      tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]);
3743  }
3744
3745  Pi [0] += tmp[0]*xToJ*F.mvar();
3746
3747  int degPi, degBuf;
3748  for (int l= 1; l < factors.length() - 1; l++)
3749  {
3750    degPi= degree (Pi [l - 1], x);
3751    degBuf= degree (bufFactors[l + 1], x);
3752    if (degPi > 0 && degBuf > 0)
3753    {
3754      M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
3755      if (j + 2 <= M.rows())
3756        M (j + 2, l + 1)= mulNTL (Pi [l - 1][j + 1], bufFactors[l + 1] [j + 1]);
3757    }
3758
3759    if (degPi > 0 && degBuf > 0)
3760      uIZeroJ= mulNTL (Pi[l -1] [0], buf[l + 1]) +
3761               mulNTL (uIZeroJ, bufFactors[l+1] [0]);
3762    else if (degPi > 0)
3763      uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]);
3764    else if (degBuf > 0)
3765      uIZeroJ= mulNTL (Pi[l - 1], buf[1]);
3766    else
3767      uIZeroJ= 0;
3768
3769    Pi [l] += xToJ*uIZeroJ;
3770
3771    one= bufFactors [l + 1];
3772    two= Pi [l - 1];
3773    if (degBuf > 0 && degPi > 0)
3774    {
3775      while (one.hasTerms() && one.exp() > j) one++;
3776      while (two.hasTerms() && two.exp() > j) two++;
3777      for (k= 1; k <= (int) ceil (j/2.0); k++)
3778      {
3779        if (k != j - k + 1)
3780        {
3781          if ((one.hasTerms() && one.exp() == j - k + 1) &&
3782              (two.hasTerms() && two.exp() == j - k + 1))
3783          {
3784            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
3785                      (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1) -
3786                      M (j - k + 2, l + 1);
3787            one++;
3788            two++;
3789          }
3790          else if (one.hasTerms() && one.exp() == j - k + 1)
3791          {
3792            tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
3793                               Pi[l - 1] [k]) - M (k + 1, l + 1);
3794            one++;
3795          }
3796          else if (two.hasTerms() && two.exp() == j - k + 1)
3797          {
3798            tmp[l] += mulNTL (bufFactors[l + 1] [k],
3799                      (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1);
3800            two++;
3801           }
3802        }
3803        else
3804          tmp[l] += M (k + 1, l + 1);
3805      }
3806    }
3807
3808    if (degPi >= j + 1 && degBuf >= j + 1)
3809    {
3810      if (j + 2 <= M.rows())
3811        tmp [l] += mulNTL ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
3812                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
3813                          ) - M(1,l+1) - M (j + 2,l+1);
3814    }
3815    else if (degPi >= j + 1)
3816    {
3817      if (degBuf > 0)
3818        tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1] [0]);
3819      else
3820        tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1]);
3821    }
3822    else if (degBuf >= j + 1)
3823    {
3824      if (degPi > 0)
3825        tmp[l] += mulNTL (Pi [l - 1] [0], bufFactors [l + 1] [j + 1]);
3826      else
3827        tmp[l] += mulNTL (Pi [l - 1], bufFactors [l + 1] [j + 1]);
3828    }
3829
3830    Pi[l] += tmp[l]*xToJ*F.mvar();
3831  }
3832  return;
3833}
3834
3835void
3836henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
3837              CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
3838{
3839  if (sort)
3840    sortList (factors, Variable (1));
3841  Pi= CFArray (factors.length() - 2);
3842  CFList bufFactors2= factors;
3843  bufFactors2.removeFirst();
3844  diophant= diophantine (F[0], bufFactors2);
3845  DEBOUTLN (cerr, "diophant= " << diophant);
3846
3847  CFArray bufFactors= CFArray (bufFactors2.length());
3848  int i= 0;
3849  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
3850    bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
3851
3852  Variable x= F.mvar();
3853  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
3854  {
3855    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
3856    Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
3857                        mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
3858  }
3859  else if (degree (bufFactors[0], x) > 0)
3860  {
3861    M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
3862    Pi [0]= M (1, 1) +
3863            mulNTL (bufFactors [0] [1], bufFactors[1])*x;
3864  }
3865  else if (degree (bufFactors[1], x) > 0)
3866  {
3867    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
3868    Pi [0]= M (1, 1) +
3869            mulNTL (bufFactors [0], bufFactors[1] [1])*x;
3870  }
3871  else
3872  {
3873    M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
3874    Pi [0]= M (1, 1);
3875  }
3876
3877  for (i= 1; i < Pi.size(); i++)
3878  {
3879    if (degree (Pi[i-1], x) > 0 && degree (bufFactors [i+1], x) > 0)
3880    {
3881      M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors[i+1] [0]);
3882      Pi [i]= M (1,i+1) + (mulNTL (Pi[i-1] [1], bufFactors[i+1] [0]) +
3883                       mulNTL (Pi[i-1] [0], bufFactors [i+1] [1]))*x;
3884    }
3885    else if (degree (Pi[i-1], x) > 0)
3886    {
3887      M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors [i+1]);
3888      Pi [i]=  M(1,i+1) + mulNTL (Pi[i-1] [1], bufFactors[i+1])*x;
3889    }
3890    else if (degree (bufFactors[i+1], x) > 0)
3891    {
3892      M (1,i+1)= mulNTL (Pi[i-1], bufFactors [i+1] [0]);
3893      Pi [i]= M (1,i+1) + mulNTL (Pi[i-1], bufFactors[i+1] [1])*x;
3894    }
3895    else
3896    {
3897      M (1,i+1)= mulNTL (Pi [i-1], bufFactors [i+1]);
3898      Pi [i]= M (1,i+1);
3899    }
3900  }
3901
3902  for (i= 1; i < l; i++)
3903    henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
3904
3905  factors= CFList();
3906  for (i= 0; i < bufFactors.size(); i++)
3907    factors.append (bufFactors[i]);
3908  return;
3909}
3910
3911
3912/// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
3913/// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
3914static inline
3915CFList
3916diophantine (const CFList& recResult, const CFList& factors,
3917             const CFList& products, const CFList& M, const CanonicalForm& E,
3918             bool& bad)
3919{
3920  if (M.isEmpty())
3921  {
3922    CFList result;
3923    CFListIterator j= factors;
3924    CanonicalForm buf;
3925    for (CFListIterator i= recResult; i.hasItem(); i++, j++)
3926    {
3927      ASSERT (E.isUnivariate() || E.inCoeffDomain(),
3928              "constant or univariate poly expected");
3929      ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
3930              "constant or univariate poly expected");
3931      ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
3932              "constant or univariate poly expected");
3933      buf= mulNTL (E, i.getItem());
3934      result.append (modNTL (buf, j.getItem()));
3935    }
3936    return result;
3937  }
3938  Variable y= M.getLast().mvar();
3939  CFList bufFactors= factors;
3940  for (CFListIterator i= bufFactors; i.hasItem(); i++)
3941    i.getItem()= mod (i.getItem(), y);
3942  CFList bufProducts= products;
3943  for (CFListIterator i= bufProducts; i.hasItem(); i++)
3944    i.getItem()= mod (i.getItem(), y);
3945  CFList buf= M;
3946  buf.removeLast();
3947  CanonicalForm bufE= mod (E, y);
3948  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
3949                                      bufE, bad);
3950
3951  if (bad)
3952    return CFList();
3953
3954  CanonicalForm e= E;
3955  CFListIterator j= products;
3956  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
3957    e -= i.getItem()*j.getItem();
3958
3959  CFList result= recDiophantine;
3960  int d= degree (M.getLast());
3961  CanonicalForm coeffE;
3962  for (int i= 1; i < d; i++)
3963  {
3964    if (degree (e, y) > 0)
3965      coeffE= e[i];
3966    else
3967      coeffE= 0;
3968    if (!coeffE.isZero())
3969    {
3970      CFListIterator k= result;
3971      recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
3972                                   coeffE, bad);
3973      if (bad)
3974        return CFList();
3975      CFListIterator l= products;
3976      for (j= recDiophantine; j.hasItem(); j++, k++, l++)
3977      {
3978        k.getItem() += j.getItem()*power (y, i);
3979        e -= j.getItem()*power (y, i)*l.getItem();
3980      }
3981    }
3982    if (e.isZero())
3983      break;
3984  }
3985  if (!e.isZero())
3986  {
3987    bad= true;
3988    return CFList();
3989  }
3990  return result;
3991}
3992
3993static inline
3994void
3995henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
3996            const CFList& diophant, CFMatrix& M, CFArray& Pi,
3997            const CFList& products, int j, const CFList& MOD, bool& noOneToOne)
3998{
3999  CanonicalForm E;
4000  CanonicalForm xToJ= power (F.mvar(), j);
4001  Variable x= F.mvar();
4002
4003  // compute the error
4004#ifdef DEBUGOUTPUT
4005    CanonicalForm test= 1;
4006    for (int i= 0; i < factors.length(); i++)
4007    {
4008      if (i == 0)
4009        test *= mod (bufFactors [i], power (x, j));
4010      else
4011        test *= bufFactors[i];
4012    }
4013    test= mod (test, power (x, j));
4014    test= mod (test, MOD);
4015    CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
4016    DEBOUTLN (cerr, "test= " << test2);
4017#endif
4018
4019  if (degree (Pi [factors.length() - 2], x) > 0)
4020    E= F[j] - Pi [factors.length() - 2] [j];
4021  else
4022    E= F[j];
4023
4024  CFArray buf= CFArray (diophant.length());
4025
4026  // actual lifting
4027  CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
4028                                    noOneToOne);
4029
4030  if (noOneToOne)
4031    return;
4032
4033  int k= 0;
4034  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
4035  {
4036    buf[k]= i.getItem();
4037    bufFactors[k] += xToJ*i.getItem();
4038  }
4039
4040  // update Pi [0]
4041  int degBuf0= degree (bufFactors[0], x);
4042  int degBuf1= degree (bufFactors[1], x);
4043  if (degBuf0 > 0 && degBuf1 > 0)
4044  {
4045    M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
4046    if (j + 2 <= M.rows())
4047      M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
4048  }
4049  CanonicalForm uIZeroJ;
4050
4051  if (degBuf0 > 0 && degBuf1 > 0)
4052    uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
4053             mulMod (bufFactors[1] [0], buf[0], MOD);
4054  else if (degBuf0 > 0)
4055    uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
4056  else if (degBuf1 > 0)
4057    uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
4058  else
4059    uIZeroJ= 0;
4060  Pi [0] += xToJ*uIZeroJ;
4061
4062  CFArray tmp= CFArray (factors.length() - 1);
4063  for (k= 0; k < factors.length() - 1; k++)
4064    tmp[k]= 0;
4065  CFIterator one, two;
4066  one= bufFactors [0];
4067  two= bufFactors [1];
4068  if (degBuf0 > 0 && degBuf1 > 0)
4069  {
4070    while (one.hasTerms() && one.exp() > j) one++;
4071    while (two.hasTerms() && two.exp() > j) two++;
4072    for (k= 1; k <= (int) ceil (j/2.0); k++)
4073    {
4074      if (k != j - k + 1)
4075      {
4076        if ((one.hasTerms() && one.exp() == j - k + 1) &&
4077            (two.hasTerms() && two.exp() == j - k + 1))
4078        {
4079          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
4080                    (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
4081                    M (j - k + 2, 1);
4082          one++;
4083          two++;
4084        }
4085        else if (one.hasTerms() && one.exp() == j - k + 1)
4086        {
4087          tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
4088                    bufFactors[1] [k], MOD) - M (k + 1, 1);
4089          one++;
4090        }
4091        else if (two.hasTerms() && two.exp() == j - k + 1)
4092        {
4093          tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
4094                    two.coeff()), MOD) - M (k + 1, 1);
4095          two++;
4096        }
4097      }
4098      else
4099      {
4100        tmp[0] += M (k + 1, 1);
4101      }
4102    }
4103  }
4104
4105  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
4106  {
4107    if (j + 2 <= M.rows())
4108      tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
4109                         (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
4110                         - M(1,1) - M (j + 2,1);
4111  }
4112  else if (degBuf0 >= j + 1)
4113  {
4114    if (degBuf1 > 0)
4115      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
4116    else
4117      tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
4118  }
4119  else if (degBuf1 >= j + 1)
4120  {
4121    if (degBuf0 > 0)
4122      tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
4123    else
4124      tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
4125  }
4126  Pi [0] += tmp[0]*xToJ*F.mvar();
4127
4128  // update Pi [l]
4129  int degPi, degBuf;
4130  for (int l= 1; l < factors.length() - 1; l++)
4131  {
4132    degPi= degree (Pi [l - 1], x);
4133    degBuf= degree (bufFactors[l + 1], x);
4134    if (degPi > 0 && degBuf > 0)
4135    {
4136      M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
4137      if (j + 2 <= M.rows())
4138        M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
4139                                  MOD);
4140    }
4141
4142    if (degPi > 0 && degBuf > 0)
4143      uIZeroJ= mulMod (Pi[l -1] [0], buf[l + 1], MOD) +
4144               mulMod (uIZeroJ, bufFactors[l+1] [0], MOD);
4145    else if (degPi > 0)
4146      uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD);
4147    else if (degBuf > 0)
4148      uIZeroJ= mulMod (Pi[l - 1], buf[1], MOD);
4149    else
4150      uIZeroJ= 0;
4151
4152    Pi [l] += xToJ*uIZeroJ;
4153
4154    one= bufFactors [l + 1];
4155    two= Pi [l - 1];
4156    if (degBuf > 0 && degPi > 0)
4157    {
4158      while (one.hasTerms() && one.exp() > j) one++;
4159      while (two.hasTerms() && two.exp() > j) two++;
4160      for (k= 1; k <= (int) ceil (j/2.0); k++)
4161      {
4162        if (k != j - k + 1)
4163        {
4164          if ((one.hasTerms() && one.exp() == j - k + 1) &&
4165              (two.hasTerms() && two.exp() == j - k + 1))
4166          {
4167            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
4168                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
4169                      M (j - k + 2, l + 1);
4170            one++;
4171            two++;
4172          }
4173          else if (one.hasTerms() && one.exp() == j - k + 1)
4174          {
4175            tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
4176                      Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
4177            one++;
4178          }
4179          else if (two.hasTerms() && two.exp() == j - k + 1)
4180          {
4181            tmp[l] += mulMod (bufFactors[l + 1] [k],
4182                      (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
4183            two++;
4184           }
4185        }
4186        else
4187          tmp[l] += M (k + 1, l + 1);
4188      }
4189    }
4190
4191    if (degPi >= j + 1 && degBuf >= j + 1)
4192    {
4193      if (j + 2 <= M.rows())
4194        tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
4195                           (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
4196                            , MOD) - M(1,l+1) - M (j + 2,l+1);
4197    }
4198    else if (degPi >= j + 1)
4199    {
4200      if (degBuf > 0)
4201        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
4202      else
4203        tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
4204    }
4205    else if (degBuf >= j + 1)
4206    {
4207      if (degPi > 0)
4208        tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
4209      else
4210        tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
4211    }
4212
4213    Pi[l] += tmp[l]*xToJ*F.mvar();
4214  }
4215  return;
4216}
4217
4218// wrt. Variable (1)
4219CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
4220{
4221  if (degree (F, 1) <= 0)
4222   return c;
4223  else
4224  {
4225    CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
4226    result += (swapvar (c, Variable (F.level() + 1), Variable (1))
4227              - LC (result))*power (result.mvar(), degree (result));
4228    return swapvar (result, Variable (F.level() + 1), Variable (1));
4229  }
4230}
4231
4232CFList
4233henselLift232 (const CFList& eval, const CFList& factors, int* l, CFList&
4234              diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
4235              const CFList& LCs2, bool& bad)
4236{
4237  CFList buf= factors;
4238  int k= 0;
4239  int liftBoundBivar= l[k];
4240  CFList bufbuf= factors;
4241  Variable v= Variable (2);
4242
4243  CFList MOD;
4244  MOD.append (power (Variable (2), liftBoundBivar));
4245  CFArray bufFactors= CFArray (factors.length());
4246  k= 0;
4247  CFListIterator j= eval;
4248  j++;
4249  CFListIterator iter1= LCs1;
4250  CFListIterator iter2= LCs2;
4251  iter1++;
4252  iter2++;
4253  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
4254  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
4255
4256  CFListIterator i= buf;
4257  i++;
4258  Variable y= j.getItem().mvar();
4259  if (y.level() != 3)
4260    y= Variable (3);
4261
4262  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
4263  M (1, 1)= Pi[0];
4264  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
4265    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
4266                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
4267  else if (degree (bufFactors[0], y) > 0)
4268    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
4269  else if (degree (bufFactors[1], y) > 0)
4270    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
4271
4272  CFList products;
4273  for (int i= 0; i < bufFactors.size(); i++)
4274  {
4275    if (degree (bufFactors[i], y) > 0)
4276      products.append (eval.getFirst()/bufFactors[i] [0]);
4277    else
4278      products.append (eval.getFirst()/bufFactors[i]);
4279  }
4280
4281  for (int d= 1; d < l[1]; d++)
4282  {
4283    henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
4284    if (bad)
4285      return CFList();
4286  }
4287  CFList result;
4288  for (k= 0; k < factors.length(); k++)
4289    result.append (bufFactors[k]);
4290  return result;
4291}
4292
4293
4294CFList
4295henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
4296            diophant, CFArray& Pi, CFMatrix& M, const int lOld, int&
4297            lNew, const CFList& LCs1, const CFList& LCs2, bool& bad)
4298{
4299  int k= 0;
4300  CFArray bufFactors= CFArray (factors.length());
4301  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
4302  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
4303  CFList buf= factors;
4304  Variable y= F.getLast().mvar();
4305  Variable x= F.getFirst().mvar();
4306  CanonicalForm xToLOld= power (x, lOld);
4307  Pi [0]= mod (Pi[0], xToLOld);
4308  M (1, 1)= Pi [0];
4309
4310  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
4311    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
4312                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
4313  else if (degree (bufFactors[0], y) > 0)
4314    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
4315  else if (degree (bufFactors[1], y) > 0)
4316    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
4317
4318  CFList products;
4319  CanonicalForm quot;
4320  for (int i= 0; i < bufFactors.size(); i++)
4321  {
4322    if (degree (bufFactors[i], y) > 0)
4323    {
4324      if (!fdivides (bufFactors[i] [0], F.getFirst(), quot))
4325      {
4326        bad= true;
4327        return CFList();
4328      }
4329      products.append (quot);
4330    }
4331    else
4332    {
4333      if (!fdivides (bufFactors[i], F.getFirst(), quot))
4334      {
4335        bad= true;
4336        return CFList();
4337      }
4338      products.append (quot);
4339    }
4340  }
4341
4342  for (int d= 1; d < lNew; d++)
4343  {
4344    henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD, bad);
4345    if (bad)
4346      return CFList();
4347  }
4348
4349  CFList result;
4350  for (k= 0; k < factors.length(); k++)
4351    result.append (bufFactors[k]);
4352  return result;
4353}
4354
4355CFList
4356henselLift2 (const CFList& eval, const CFList& factors, int* l, const int
4357             lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
4358             const CFArray& Pi, const CFList& diophant, bool& bad)
4359{
4360  CFList bufDiophant= diophant;
4361  CFList buf= factors;
4362  if (sort)
4363    sortList (buf, Variable (1));
4364  CFArray bufPi= Pi;
4365  CFMatrix M= CFMatrix (l[1], factors.length());
4366  CFList result= henselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2,
4367                               bad);
4368  if (bad)
4369    return CFList();
4370
4371  if (eval.length() == 2)
4372    return result;
4373  CFList MOD;
4374  for (int i= 0; i < 2; i++)
4375    MOD.append (power (Variable (i + 2), l[i]));
4376  CFListIterator j= eval;
4377  j++;
4378  CFList bufEval;
4379  bufEval.append (j.getItem());
4380  j++;
4381  CFListIterator jj= LCs1;
4382  CFListIterator jjj= LCs2;
4383  CFList bufLCs1, bufLCs2;
4384  jj++, jjj++;
4385  bufLCs1.append (jj.getItem());
4386  bufLCs2.append (jjj.getItem());
4387  jj++, jjj++;
4388
4389  for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
4390  {
4391    bufEval.append (j.getItem());
4392    bufLCs1.append (jj.getItem());
4393    bufLCs2.append (jjj.getItem());
4394    M= CFMatrix (l[i], factors.length());
4395    result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
4396                         l[i], bufLCs1, bufLCs2, bad);
4397    if (bad)
4398      return CFList();
4399    MOD.append (power (Variable (i + 2), l[i]));
4400    bufEval.removeFirst();
4401    bufLCs1.removeFirst();
4402    bufLCs2.removeFirst();
4403  }
4404  return result;
4405}
4406
4407CFList
4408nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
4409                      CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
4410                      int bivarLiftBound, bool& bad)
4411{
4412  CFList bufFactors2= factors;
4413
4414  Variable y= Variable (2);
4415  for (CFListIterator i= bufFactors2; i.hasItem(); i++)
4416    i.getItem()= mod (i.getItem(), y);
4417
4418  CanonicalForm bufF= F;
4419  bufF= mod (bufF, y);
4420  bufF= mod (bufF, Variable (3));
4421
4422  diophant= diophantine (bufF, bufFactors2);
4423
4424  CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
4425
4426  Pi= CFArray (bufFactors2.length() - 1);
4427
4428  CFArray bufFactors= CFArray (bufFactors2.length());
4429  CFListIterator j= LCs;
4430  int i= 0;
4431  for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
4432    bufFactors[i]= replaceLC (k.getItem(), j.getItem());
4433
4434  //initialise Pi
4435  Variable v= Variable (3);
4436  CanonicalForm yToL= power (y, bivarLiftBound);
4437  if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
4438  {
4439    M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
4440    Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
4441                       mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
4442  }
4443  else if (degree (bufFactors[0], v) > 0)
4444  {
4445    M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
4446    Pi [0]=  M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
4447  }
4448  else if (degree (bufFactors[1], v) > 0)
4449  {
4450    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
4451    Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
4452  }
4453  else
4454  {
4455    M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
4456    Pi [0]= M (1,1);
4457  }
4458
4459  for (i= 1; i < Pi.size(); i++)
4460  {
4461    if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
4462    {
4463      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
4464      Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
4465                       mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
4466    }
4467    else if (degree (Pi[i-1], v) > 0)
4468    {
4469      M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
4470      Pi [i]=  M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
4471    }
4472    else if (degree (bufFactors[i+1], v) > 0)
4473    {
4474      M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
4475      Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
4476    }
4477    else
4478    {
4479      M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
4480      Pi [i]= M (1,i+1);
4481    }
4482  }
4483
4484  CFList products;
4485  bufF= mod (F, Variable (3));
4486  for (CFListIterator k= factors; k.hasItem(); k++)
4487    products.append (bufF/k.getItem());
4488
4489  CFList MOD= CFList (power (v, liftBound));
4490  MOD.insert (yToL);
4491  for (int d= 1; d < liftBound; d++)
4492  {
4493    henselStep2 (F, factors, bufFactors, diophant, M, Pi, products, d, MOD, bad);
4494    if (bad)
4495      return CFList();
4496  }
4497
4498  CFList result;
4499  for (i= 0; i < factors.length(); i++)
4500    result.append (bufFactors[i]);
4501  return result;
4502}
4503
4504CFList
4505nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
4506                    CFList& diophant, CFArray& Pi, CFMatrix& M, const int lOld,
4507                    int& lNew, const CFList& MOD, bool& noOneToOne
4508                   )
4509{
4510
4511  int k= 0;
4512  CFArray bufFactors= CFArray (factors.length());
4513  CFListIterator j= LCs;
4514  for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
4515    bufFactors [k]= replaceLC (i.getItem(), j.getItem());
4516
4517  Variable y= F.getLast().mvar();
4518  Variable x= F.getFirst().mvar();
4519  CanonicalForm xToLOld= power (x, lOld);
4520
4521  Pi [0]= mod (Pi[0], xToLOld);
4522  M (1, 1)= Pi [0];
4523
4524  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
4525    Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
4526                        mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
4527  else if (degree (bufFactors[0], y) > 0)
4528    Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
4529  else if (degree (bufFactors[1], y) > 0)
4530    Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
4531
4532  for (int i= 1; i < Pi.size(); i++)
4533  {
4534    Pi [i]= mod (Pi [i], xToLOld);
4535    M (1, i + 1)= Pi [i];
4536
4537    if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
4538      Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
4539                 mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
4540    else if (degree (Pi[i-1], y) > 0)
4541      Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
4542    else if (degree (bufFactors[i+1], y) > 0)
4543      Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
4544  }
4545
4546  CFList products;
4547  CanonicalForm quot, bufF= F.getFirst();
4548
4549  for (int i= 0; i < bufFactors.size(); i++)
4550  {
4551    if (degree (bufFactors[i], y) > 0)
4552    {
4553      if (!fdivides (bufFactors[i] [0], bufF, quot))
4554      {
4555        noOneToOne= true;
4556        return factors;
4557      }
4558      products.append (quot);
4559    }
4560    else
4561    {
4562      if (!fdivides (bufFactors[i], bufF, quot))
4563      {
4564        noOneToOne= true;
4565        return factors;
4566      }
4567      products.append (quot);
4568    }
4569  }
4570
4571  for (int d= 1; d < lNew; d++)
4572  {
4573    henselStep2 (F.getLast(), factors, bufFactors, diophant, M, Pi, products, d,
4574                 MOD, noOneToOne);
4575    if (noOneToOne)
4576      return CFList();
4577  }
4578
4579  CFList result;
4580  for (k= 0; k < factors.length(); k++)
4581    result.append (bufFactors[k]);
4582  return result;
4583}
4584
4585CFList
4586nonMonicHenselLift (const CFList& eval, const CFList& factors,
4587                    CFList* const& LCs, CFList& diophant, CFArray& Pi,
4588                    int* liftBound, int length, bool& noOneToOne
4589                   )
4590{
4591  CFList bufDiophant= diophant;
4592  CFList buf= factors;
4593  CFArray bufPi= Pi;
4594  CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
4595  int k= 0;
4596
4597  CFList result=
4598  nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
4599                        liftBound[1], liftBound[0], noOneToOne);
4600
4601  if (noOneToOne)
4602    return CFList();
4603
4604  if (eval.length() == 1)
4605    return result;
4606
4607  k++;
4608  CFList MOD;
4609  for (int i= 0; i < 2; i++)
4610    MOD.append (power (Variable (i + 2), liftBound[i]));
4611
4612  CFListIterator j= eval;
4613  CFList bufEval;
4614  bufEval.append (j.getItem());
4615  j++;
4616
4617  for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
4618  {
4619    bufEval.append (j.getItem());
4620    M= CFMatrix (liftBound[i], factors.length() - 1);
4621    result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
4622                                liftBound[i-1], liftBound[i], MOD, noOneToOne);
4623    if (noOneToOne)
4624      return result;
4625    MOD.append (power (Variable (i + 2), liftBound[i]));
4626    bufEval.removeFirst();
4627  }
4628
4629  return result;
4630}
4631
4632#endif
4633/* HAVE_NTL */
4634
Note: See TracBrowser for help on using the repository browser.