source: git/factory/facHensel.cc @ 3ef2d6

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