facHensel.cc
Go to the documentation of this file.
1 /*****************************************************************************\
2  * Computer Algebra System SINGULAR
3 \*****************************************************************************/
4 /** @file facHensel.cc
5  *
6  * This file implements functions to lift factors via Hensel lifting.
7  *
8  * ABSTRACT: Hensel lifting is described in "Efficient Multivariate
9  * Factorization over Finite Fields" by L. Bernardin & M. Monagon and
10  * "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn
11  *
12  * @author Martin Lee
13  *
14  **/
15 /*****************************************************************************/
16 
17 
18 #include "config.h"
19 
20 
21 #include "cf_assert.h"
22 #include "debug.h"
23 #include "timing.h"
24 
25 #include "cfGcdAlgExt.h"
26 #include "facHensel.h"
27 #include "facMul.h"
28 #include "fac_util.h"
29 #include "cf_algorithm.h"
30 #include "cf_primes.h"
31 #include "facBivar.h"
32 #include "cfNTLzzpEXGCD.h"
33 #include "cfUnivarGcd.h"
34 
35 #ifdef HAVE_NTL
36 #include <NTL/lzz_pEX.h>
37 #include "NTLconvert.h"
38 
39 #ifdef HAVE_FLINT
40 #include "FLINTconvert.h"
41 #endif
42 
43 TIMING_DEFINE_PRINT (diotime)
44 TIMING_DEFINE_PRINT (product1)
45 TIMING_DEFINE_PRINT (product2)
46 TIMING_DEFINE_PRINT (hensel23)
47 TIMING_DEFINE_PRINT (hensel)
48 
49 #if (!(HAVE_FLINT && __FLINT_RELEASE >= 20400))
50 static
51 CFList productsNTL (const CFList& factors, const CanonicalForm& M)
52 {
54  {
56  zz_p::init (getCharacteristic());
57  }
58  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
59  zz_pE::init (NTLMipo);
60  zz_pEX prod;
61  vec_zz_pEX v;
62  v.SetLength (factors.length());
63  int j= 0;
64  for (CFListIterator i= factors; i.hasItem(); i++, j++)
65  {
66  if (i.getItem().inCoeffDomain())
67  v[j]= to_zz_pEX (to_zz_pE (convertFacCF2NTLzzpX (i.getItem())));
68  else
69  v[j]= convertFacCF2NTLzz_pEX (i.getItem(), NTLMipo);
70  }
71  CFList result;
72  Variable x= Variable (1);
73  for (int j= 0; j < factors.length(); j++)
74  {
75  int k= 0;
76  set(prod);
77  for (int i= 0; i < factors.length(); i++, k++)
78  {
79  if (k == j)
80  continue;
81  prod *= v[i];
82  }
83  result.append (convertNTLzz_pEX2CF (prod, x, M.mvar()));
84  }
85  return result;
86 }
87 #endif
88 
89 #if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
90 static
91 CFList productsFLINT (const CFList& factors, const CanonicalForm& M)
92 {
93  nmod_poly_t FLINTmipo;
94  fq_nmod_ctx_t fq_con;
95  fq_nmod_poly_t prod;
96  fq_nmod_t buf;
97 
98  nmod_poly_init (FLINTmipo, getCharacteristic());
99  convertFacCF2nmod_poly_t (FLINTmipo, M);
100 
101  fq_nmod_ctx_init_modulus (fq_con, FLINTmipo, "Z");
102 
103  fq_nmod_poly_t * vec=new fq_nmod_poly_t [factors.length()];
104 
105  int j= 0;
106 
107  for (CFListIterator i= factors; i.hasItem(); i++, j++)
108  {
109  if (i.getItem().inCoeffDomain())
110  {
111  fq_nmod_poly_init (vec[j], fq_con);
112  fq_nmod_init2 (buf, fq_con);
113  convertFacCF2Fq_nmod_t (buf, i.getItem(), fq_con);
114  fq_nmod_poly_set_coeff (vec[j], 0, buf, fq_con);
115  fq_nmod_clear (buf, fq_con);
116  }
117  else
118  convertFacCF2Fq_nmod_poly_t (vec[j], i.getItem(), fq_con);
119  }
120 
123  fq_nmod_poly_init (prod, fq_con);
124  for (j= 0; j < factors.length(); j++)
125  {
126  int k= 0;
127  fq_nmod_poly_one (prod, fq_con);
128  for (int i= 0; i < factors.length(); i++, k++)
129  {
130  if (k == j)
131  continue;
132  fq_nmod_poly_mul (prod, prod, vec[i], fq_con);
133  }
134  result.append (convertFq_nmod_poly_t2FacCF (prod, x, M.mvar(), fq_con));
135  }
136  for (j= 0; j < factors.length(); j++)
137  fq_nmod_poly_clear (vec[j], fq_con);
138 
139  nmod_poly_clear (FLINTmipo);
140  fq_nmod_poly_clear (prod, fq_con);
141  fq_nmod_ctx_clear (fq_con);
142  delete [] vec;
143  return result;
144 }
145 #endif
146 
147 static
148 void tryDiophantine (CFList& result, const CanonicalForm& F,
149  const CFList& factors, const CanonicalForm& M, bool& fail)
150 {
151  ASSERT (M.isUnivariate(), "expected univariate poly");
152 
153  CFList bufFactors= factors;
154  bufFactors.removeFirst();
155  bufFactors.insert (factors.getFirst () (0,2));
156  CanonicalForm inv, leadingCoeff= Lc (F);
157  CFListIterator i= bufFactors;
158  if (bufFactors.getFirst().inCoeffDomain())
159  {
160  if (i.hasItem())
161  i++;
162  }
163  for (; i.hasItem(); i++)
164  {
165  tryInvert (Lc (i.getItem()), M, inv ,fail);
166  if (fail)
167  return;
168  i.getItem()= reduce (i.getItem()*inv, M);
169  }
170 #if (HAVE_FLINT && __FLINT_RELEASE >= 20400)
171  bufFactors= productsFLINT (bufFactors, M);
172 #else
173  bufFactors= productsNTL (bufFactors, M);
174 #endif
175 
176  CanonicalForm buf1, buf2, buf3, S, T;
177  i= bufFactors;
178  if (i.hasItem())
179  i++;
180  buf1= bufFactors.getFirst();
181  buf2= i.getItem();
182 #ifdef HAVE_NTL
183  Variable x= Variable (1);
185  {
187  zz_p::init (getCharacteristic());
188  }
189  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
190  zz_pE::init (NTLMipo);
191  zz_pEX NTLbuf1, NTLbuf2, NTLbuf3, NTLS, NTLT;
192  NTLbuf1= convertFacCF2NTLzz_pEX (buf1, NTLMipo);
193  NTLbuf2= convertFacCF2NTLzz_pEX (buf2, NTLMipo);
194  tryNTLXGCD (NTLbuf3, NTLS, NTLT, NTLbuf1, NTLbuf2, fail);
195  if (fail)
196  return;
197  S= convertNTLzz_pEX2CF (NTLS, x, M.mvar());
198  T= convertNTLzz_pEX2CF (NTLT, x, M.mvar());
199 #else
200  tryExtgcd (buf1, buf2, M, buf3, S, T, fail);
201  if (fail)
202  return;
203 #endif
204  result.append (S);
205  result.append (T);
206  if (i.hasItem())
207  i++;
208  for (; i.hasItem(); i++)
209  {
210 #ifdef HAVE_NTL
211  NTLbuf1= convertFacCF2NTLzz_pEX (i.getItem(), NTLMipo);
212  tryNTLXGCD (NTLbuf3, NTLS, NTLT, NTLbuf3, NTLbuf1, fail);
213  if (fail)
214  return;
215  S= convertNTLzz_pEX2CF (NTLS, x, M.mvar());
216  T= convertNTLzz_pEX2CF (NTLT, x, M.mvar());
217 #else
218  buf1= i.getItem();
219  tryExtgcd (buf3, buf1, M, buf3, S, T, fail);
220  if (fail)
221  return;
222 #endif
223  CFListIterator k= factors;
224  for (CFListIterator j= result; j.hasItem(); j++, k++)
225  {
226  j.getItem() *= S;
227  j.getItem()= mod (j.getItem(), k.getItem());
228  j.getItem()= reduce (j.getItem(), M);
229  }
230  result.append (T);
231  }
232 }
233 
234 static
236 {
237  CFList result;
238  for (CFListIterator i= L; i.hasItem(); i++)
239  result.append (mapinto (i.getItem()));
240  return result;
241 }
242 
243 static
244 int mod (const CFList& L, const CanonicalForm& p)
245 {
246  for (CFListIterator i= L; i.hasItem(); i++)
247  {
248  if (mod (i.getItem(), p) == 0)
249  return 0;
250  }
251  return 1;
252 }
253 
254 
255 static void
256 chineseRemainder (const CFList & x1, const CanonicalForm & q1,
257  const CFList & x2, const CanonicalForm & q2,
258  CFList & xnew, CanonicalForm & qnew)
259 {
260  ASSERT (x1.length() == x2.length(), "expected lists of equal length");
262  CFListIterator j= x2;
263  for (CFListIterator i= x1; i.hasItem() && j.hasItem(); i++, j++)
264  {
265  chineseRemainder (i.getItem(), q1, j.getItem(), q2, tmp1, tmp2);
266  xnew.append (tmp1);
267  }
268  qnew= tmp2;
269 }
270 
271 static
272 CFList Farey (const CFList& L, const CanonicalForm& q)
273 {
274  CFList result;
275  for (CFListIterator i= L; i.hasItem(); i++)
276  result.append (Farey (i.getItem(), q));
277  return result;
278 }
279 
280 static
281 CFList replacevar (const CFList& L, const Variable& a, const Variable& b)
282 {
283  CFList result;
284  for (CFListIterator i= L; i.hasItem(); i++)
285  result.append (replacevar (i.getItem(), a, b));
286  return result;
287 }
288 
289 CFList
290 modularDiophant (const CanonicalForm& f, const CFList& factors,
291  const CanonicalForm& M)
292 {
293  bool save_rat=!isOn (SW_RATIONAL);
294  On (SW_RATIONAL);
295  CanonicalForm F= f*bCommonDen (f);
296  CFList products= factors;
297  for (CFListIterator i= products; i.hasItem(); i++)
298  {
299  if (products.getFirst().level() == 1)
300  i.getItem() /= Lc (i.getItem());
301  i.getItem() *= bCommonDen (i.getItem());
302  }
303  if (products.getFirst().level() == 1)
304  products.insert (Lc (F));
306  CFList leadingCoeffs;
307  leadingCoeffs.append (lc (F));
308  CanonicalForm dummy;
309  for (CFListIterator i= products; i.hasItem(); i++)
310  {
311  leadingCoeffs.append (lc (i.getItem()));
312  dummy= maxNorm (i.getItem());
313  bound= (dummy > bound) ? dummy : bound;
314  }
315  bound *= maxNorm (Lc (F))*maxNorm (Lc(F))*bound;
316  bound *= bound*bound;
317  bound= power (bound, degree (M));
318  bound *= power (CanonicalForm (2),degree (f));
319  CanonicalForm bufBound= bound;
320  int i = cf_getNumBigPrimes() - 1;
321  int p;
322  CFList resultModP, result, newResult;
323  CanonicalForm q (0), newQ;
324  bool fail= false;
325  Variable a= M.mvar();
326  Variable b= Variable (2);
327  setReduce (M.mvar(), false);
328  CanonicalForm mipo= bCommonDen (M)*M;
329  Off (SW_RATIONAL);
330  CanonicalForm modMipo;
331  leadingCoeffs.append (lc (mipo));
332  CFList tmp1, tmp2;
333  bool equal= false;
334  int count= 0;
335  do
336  {
337  p = cf_getBigPrime( i );
338  i--;
339  while ( i >= 0 && mod( leadingCoeffs, p ) == 0)
340  {
341  p = cf_getBigPrime( i );
342  i--;
343  }
344 
345  ASSERT (i >= 0, "ran out of primes"); //sic
346 
347  setCharacteristic (p);
348  modMipo= mapinto (mipo);
349  modMipo /= lc (modMipo);
350  resultModP= CFList();
351  tryDiophantine (resultModP, mapinto (F), mapinto (products), modMipo, fail);
352  setCharacteristic (0);
353  if (fail)
354  {
355  fail= false;
356  continue;
357  }
358 
359  if ( q.isZero() )
360  {
361  result= replacevar (mapinto(resultModP), a, b);
362  q= p;
363  }
364  else
365  {
366  result= replacevar (result, a, b);
367  newResult= CFList();
368  chineseRemainder( result, q, replacevar (mapinto (resultModP), a, b),
369  p, newResult, newQ );
370  q= newQ;
371  result= newResult;
372  if (newQ > bound)
373  {
374  count++;
375  tmp1= replacevar (Farey (result, q), b, a);
376  if (tmp2.isEmpty())
377  tmp2= tmp1;
378  else
379  {
380  equal= true;
381  CFListIterator k= tmp1;
382  for (CFListIterator j= tmp2; j.hasItem(); j++, k++)
383  {
384  if (j.getItem() != k.getItem())
385  equal= false;
386  }
387  if (!equal)
388  tmp2= tmp1;
389  }
390  if (count > 2)
391  {
392  bound *= bufBound;
393  equal= false;
394  count= 0;
395  }
396  }
397  if (newQ > bound && equal)
398  {
399  On( SW_RATIONAL );
400  CFList bufResult= result;
401  result= tmp2;
402  setReduce (M.mvar(), true);
403  if (factors.getFirst().level() == 1)
404  {
405  result.removeFirst();
406  CFListIterator j= factors;
407  CanonicalForm denf= bCommonDen (f);
408  for (CFListIterator i= result; i.hasItem(); i++, j++)
409  i.getItem() *= Lc (j.getItem())*denf;
410  }
411  if (factors.getFirst().level() != 1 &&
412  !bCommonDen (factors.getFirst()).isOne())
413  {
414  CanonicalForm denFirst= bCommonDen (factors.getFirst());
415  for (CFListIterator i= result; i.hasItem(); i++)
416  i.getItem() *= denFirst;
417  }
418 
419  CanonicalForm test= 0;
420  CFListIterator jj= factors;
421  for (CFListIterator ii= result; ii.hasItem(); ii++, jj++)
422  test += ii.getItem()*(f/jj.getItem());
423  if (!test.isOne())
424  {
425  bound *= bufBound;
426  equal= false;
427  count= 0;
428  setReduce (M.mvar(), false);
429  result= bufResult;
430  Off (SW_RATIONAL);
431  }
432  else
433  break;
434  }
435  }
436  } while (1);
437  if (save_rat) Off(SW_RATIONAL);
438  return result;
439 }
440 
441 void sortList (CFList& list, const Variable& x)
442 {
443  int l= 1;
444  int k= 1;
447  for (CFListIterator i= list; l <= list.length(); i++, l++)
448  {
449  for (CFListIterator j= list; k <= list.length() - l; k++)
450  {
451  m= j;
452  m++;
453  if (degree (j.getItem(), x) > degree (m.getItem(), x))
454  {
455  buf= m.getItem();
456  m.getItem()= j.getItem();
457  j.getItem()= buf;
458  j++;
459  j.getItem()= m.getItem();
460  }
461  else
462  j++;
463  }
464  k= 1;
465  }
466 }
467 
468 CFList
469 diophantine (const CanonicalForm& F, const CFList& factors);
470 
471 CFList
472 diophantineHensel (const CanonicalForm & F, const CFList& factors,
473  const modpk& b)
474 {
475  int p= b.getp();
476  setCharacteristic (p);
477  CFList recResult= diophantine (mapinto (F), mapinto (factors));
478  setCharacteristic (0);
479  recResult= mapinto (recResult);
480  CanonicalForm e= 1;
481  CFList L;
482  CFArray bufFactors= CFArray (factors.length());
483  int k= 0;
484  for (CFListIterator i= factors; i.hasItem(); i++, k++)
485  {
486  if (k == 0)
487  bufFactors[k]= i.getItem() (0);
488  else
489  bufFactors [k]= i.getItem();
490  }
491  CanonicalForm tmp, quot;
492  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
493  {
494  tmp= 1;
495  for (int l= 0; l < factors.length(); l++)
496  {
497  if (l == k)
498  continue;
499  else
500  {
501  tmp= mulNTL (tmp, bufFactors[l]);
502  }
503  }
504  L.append (tmp);
505  }
506 
507  setCharacteristic (p);
508  for (k= 0; k < factors.length(); k++)
509  bufFactors [k]= bufFactors[k].mapinto();
511 
512  CFListIterator j= L;
513  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
514  e= b (e - mulNTL (i.getItem(),j.getItem(), b));
515 
516  if (e.isZero())
517  return recResult;
518  CanonicalForm coeffE;
519  CFList s;
520  CFList result= recResult;
521  setCharacteristic (p);
522  recResult= mapinto (recResult);
523  setCharacteristic (0);
525  CanonicalForm modulus= p;
526  int d= b.getk();
527  modpk b2;
528  for (int i= 1; i < d; i++)
529  {
530  coeffE= div (e, modulus);
531  setCharacteristic (p);
532  coeffE= coeffE.mapinto();
533  setCharacteristic (0);
534  b2= modpk (p, d - i);
535  if (!coeffE.isZero())
536  {
538  CFListIterator l= L;
539  int ii= 0;
540  j= recResult;
541  for (; j.hasItem(); j++, k++, l++, ii++)
542  {
543  setCharacteristic (p);
544  g= modNTL (coeffE, bufFactors[ii]);
545  g= mulNTL (g, j.getItem());
546  g= modNTL (g, bufFactors[ii]);
547  setCharacteristic (0);
548  k.getItem() += g.mapinto()*modulus;
549  e -= mulNTL (g.mapinto(), b2 (l.getItem()), b2)*modulus;
550  e= b(e);
551  }
552  }
553  modulus *= p;
554  if (e.isZero())
555  break;
556  }
557 
558  return result;
559 }
560 
561 /// solve \f$ 1=\sum_{i=1}^n{\delta_{i} \prod_{j\neq i}{f_j}} \f$ mod \f$p^k\f$
562 /// over \f$ Q(\alpha) \f$ by p-adic lifting
563 CFList
565  const CFList& factors, modpk& b, const Variable& alpha)
566 {
567  bool fail= false;
568  CFList recResult;
569  CanonicalForm modMipo, mipo;
570  //here SW_RATIONAL is off
571  On (SW_RATIONAL);
572  mipo= getMipo (alpha);
573  bool mipoHasDen= false;
574  if (!bCommonDen (mipo).isOne())
575  {
576  mipo *= bCommonDen (mipo);
577  mipoHasDen= true;
578  }
579  Off (SW_RATIONAL);
580  int p= b.getp();
581  setCharacteristic (p);
582  setReduce (alpha, false);
583  while (1)
584  {
585  setCharacteristic (p);
586  modMipo= mapinto (mipo);
587  modMipo /= lc (modMipo);
588  tryDiophantine (recResult, mapinto (F), mapinto (factors), modMipo, fail);
589  if (fail)
590  {
591  int i= 0;
592  while (cf_getBigPrime (i) < p)
593  i++;
594  findGoodPrime (F, i);
595  findGoodPrime (G, i);
596  p=cf_getBigPrime(i);
597  b = coeffBound( G, p, mipo );
598  modpk bb= coeffBound (F, p, mipo );
599  if (bb.getk() > b.getk() ) b=bb;
600  fail= false;
601  }
602  else
603  break;
604  }
605  setCharacteristic (0);
606  recResult= mapinto (recResult);
607  setReduce (alpha, true);
608  CanonicalForm e= 1;
609  CFList L;
610  CFArray bufFactors= CFArray (factors.length());
611  int k= 0;
612  for (CFListIterator i= factors; i.hasItem(); i++, k++)
613  {
614  if (k == 0)
615  bufFactors[k]= i.getItem() (0);
616  else
617  bufFactors [k]= i.getItem();
618  }
619  CanonicalForm tmp;
620  On (SW_RATIONAL);
621  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
622  {
623  tmp= 1;
624  for (int l= 0; l < factors.length(); l++)
625  {
626  if (l == k)
627  continue;
628  else
629  tmp= mulNTL (tmp, bufFactors[l]);
630  }
631  L.append (tmp*bCommonDen(tmp));
632  }
633 
634  Variable gamma;
636  if (mipoHasDen)
637  {
638  modMipo= getMipo (alpha);
639  den= bCommonDen (modMipo);
640  modMipo *= den;
641  Off (SW_RATIONAL);
642  setReduce (alpha, false);
643  gamma= rootOf (b (modMipo*b.inverse (den)));
644  setReduce (alpha, true);
645  }
646 
647  setCharacteristic (p);
648  Variable beta;
649  Off (SW_RATIONAL);
650  setReduce (alpha, false);
651  modMipo= modMipo.mapinto();
652  modMipo /= lc (modMipo);
653  beta= rootOf (modMipo);
654  setReduce (alpha, true);
655 
656  setReduce (alpha, false);
657  for (k= 0; k < factors.length(); k++)
658  {
659  bufFactors [k]= bufFactors[k].mapinto();
660  bufFactors [k]= replacevar (bufFactors[k], alpha, beta);
661  }
662  setReduce (alpha, true);
664 
665  CFListIterator j= L;
666  for (;j.hasItem(); j++)
667  {
668  if (mipoHasDen)
669  j.getItem()= replacevar (b(j.getItem()*b.inverse(lc(j.getItem()))),
670  alpha, gamma);
671  else
672  j.getItem()= b(j.getItem()*b.inverse(lc(j.getItem())));
673  }
674  j= L;
675  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
676  {
677  if (mipoHasDen)
678  e= b (e - mulNTL (replacevar (i.getItem(), alpha, gamma),j.getItem(), b));
679  else
680  e= b (e - mulNTL (i.getItem(), j.getItem(), b));
681  }
682 
683  if (e.isZero())
684  {
685  if (mipoHasDen)
686  {
687  for (CFListIterator i= recResult; i.hasItem(); i++)
688  i.getItem()= replacevar (i.getItem(), alpha, gamma);
689  }
690  return recResult;
691  }
692  CanonicalForm coeffE;
693  CFList result= recResult;
694  if (mipoHasDen)
695  {
696  for (CFListIterator i= result; i.hasItem(); i++)
697  i.getItem()= replacevar (i.getItem(), alpha, gamma);
698  }
699  setCharacteristic (p);
700  setReduce (alpha, false);
701  recResult= mapinto (recResult);
702  setReduce (alpha, true);
703 
704  for (CFListIterator i= recResult; i.hasItem(); i++)
705  i.getItem()= replacevar (i.getItem(), alpha, beta);
706 
707  setCharacteristic (0);
709  CanonicalForm modulus= p;
710  int d= b.getk();
711  modpk b2;
712  for (int i= 1; i < d; i++)
713  {
714  coeffE= div (e, modulus);
715  setCharacteristic (p);
716  if (mipoHasDen)
717  setReduce (gamma, false);
718  else
719  setReduce (alpha, false);
720  coeffE= coeffE.mapinto();
721  if (mipoHasDen)
722  setReduce (gamma, true);
723  else
724  setReduce (alpha, true);
725  if (mipoHasDen)
726  coeffE= replacevar (coeffE, gamma, beta);
727  else
728  coeffE= replacevar (coeffE, alpha, beta);
729  setCharacteristic (0);
730  b2= modpk (p, d - i);
731  if (!coeffE.isZero())
732  {
734  CFListIterator l= L;
735  int ii= 0;
736  j= recResult;
737  for (; j.hasItem(); j++, k++, l++, ii++)
738  {
739  setCharacteristic (p);
740  g= modNTL (coeffE, bufFactors[ii]);
741  g= mulNTL (g, j.getItem());
742  g= modNTL (g, bufFactors[ii]);
743  setCharacteristic (0);
744  if (mipoHasDen)
745  {
746  setReduce (beta, false);
747  k.getItem() += replacevar (g.mapinto()*modulus, beta, gamma);
748  e -= mulNTL (replacevar (g.mapinto(), beta, gamma),
749  b2 (l.getItem()), b2)*modulus;
750  setReduce (beta, true);
751  }
752  else
753  {
754  setReduce (beta, false);
755  k.getItem() += replacevar (g.mapinto()*modulus, beta, alpha);
756  e -= mulNTL (replacevar (g.mapinto(), beta, alpha),
757  b2 (l.getItem()), b2)*modulus;
758  setReduce (beta, true);
759  }
760  e= b(e);
761  }
762  }
763  modulus *= p;
764  if (e.isZero())
765  break;
766  }
767 
768  return result;
769 }
770 
771 
772 /// solve \f$ 1=\sum_{i=1}^n{\delta_{i} \prod_{j\neq i}{f_j}} \f$ mod \f$p^k\f$
773 /// over \f$ Q(\alpha) \f$ by first computing mod \f$p\f$ and if no zero divisor
774 /// occurred compute it mod \f$p^k\f$
775 CFList
777  const CFList& factors, modpk& b, const Variable& alpha)
778 {
779  bool fail= false;
780  CFList recResult;
781  CanonicalForm modMipo, mipo;
782  //here SW_RATIONAL is off
783  On (SW_RATIONAL);
784  mipo= getMipo (alpha);
785  bool mipoHasDen= false;
786  if (!bCommonDen (mipo).isOne())
787  {
788  mipo *= bCommonDen (mipo);
789  mipoHasDen= true;
790  }
791  Off (SW_RATIONAL);
792  int p= b.getp();
793  setCharacteristic (p);
794  setReduce (alpha, false);
795  while (1)
796  {
797  setCharacteristic (p);
798  modMipo= mapinto (mipo);
799  modMipo /= lc (modMipo);
800  tryDiophantine (recResult, mapinto (F), mapinto (factors), modMipo, fail);
801  if (fail)
802  {
803  int i= 0;
804  while (cf_getBigPrime (i) < p)
805  i++;
806  findGoodPrime (F, i);
807  findGoodPrime (G, i);
808  p=cf_getBigPrime(i);
809  b = coeffBound( G, p, mipo );
810  modpk bb= coeffBound (F, p, mipo );
811  if (bb.getk() > b.getk() ) b=bb;
812  fail= false;
813  }
814  else
815  break;
816  }
817  setReduce (alpha, true);
818  setCharacteristic (0);
819 
820  Variable gamma= alpha;
822  if (mipoHasDen)
823  {
824  On (SW_RATIONAL);
825  modMipo= getMipo (alpha);
826  den= bCommonDen (modMipo);
827  modMipo *= den;
828  Off (SW_RATIONAL);
829  setReduce (alpha, false);
830  gamma= rootOf (b (modMipo*b.inverse (den)));
831  setReduce (alpha, true);
832  }
833 
834  Variable x= Variable (1);
835  CanonicalForm buf1, buf2, buf3, S;
836  CFList bufFactors= factors;
837  CFListIterator i= bufFactors;
838  if (mipoHasDen)
839  {
840  for (; i.hasItem(); i++)
841  i.getItem()= replacevar (i.getItem(), alpha, gamma);
842  }
843  i= bufFactors;
844  CFList result;
845  if (i.hasItem())
846  i++;
847  buf1= 0;
848  CanonicalForm Freplaced;
849  if (mipoHasDen)
850  {
851  Freplaced= replacevar (F, alpha, gamma);
852  buf2= divNTL (Freplaced, replacevar (i.getItem(), alpha, gamma), b);
853  }
854  else
855  buf2= divNTL (F, i.getItem(), b);
856  ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
857  ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (gamma)));
858  ZZ_pE::init (NTLmipo);
859  ZZ_pEX NTLS, NTLT, NTLbuf3;
860  ZZ_pEX NTLbuf1= convertFacCF2NTLZZ_pEX (buf1, NTLmipo);
861  ZZ_pEX NTLbuf2= convertFacCF2NTLZZ_pEX (buf2, NTLmipo);
862  XGCD (NTLbuf3, NTLS, NTLT, NTLbuf1, NTLbuf2);
863  result.append (b (convertNTLZZ_pEX2CF (NTLS, x, gamma)));
864  result.append (b (convertNTLZZ_pEX2CF (NTLT, x, gamma)));
865  if (i.hasItem())
866  i++;
867  for (; i.hasItem(); i++)
868  {
869  if (mipoHasDen)
870  buf1= divNTL (Freplaced, i.getItem(), b);
871  else
872  buf1= divNTL (F, i.getItem(), b);
873  XGCD (NTLbuf3, NTLS, NTLT, NTLbuf3, convertFacCF2NTLZZ_pEX (buf1, NTLmipo));
874  CFListIterator k= bufFactors;
875  S= convertNTLZZ_pEX2CF (NTLS, x, gamma);
876  for (CFListIterator j= result; j.hasItem(); j++, k++)
877  {
878  j.getItem()= mulNTL (j.getItem(), S, b);
879  j.getItem()= modNTL (j.getItem(), k.getItem(), b);
880  }
881  result.append (b (convertNTLZZ_pEX2CF (NTLT, x, gamma)));
882  }
883  return result;
884 }
885 
886 CFList
888  const CFList& factors, modpk& b)
889 {
890  if (getCharacteristic() == 0)
891  {
892  Variable v;
893  bool hasAlgVar= hasFirstAlgVar (F, v);
894  for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++)
895  hasAlgVar= hasFirstAlgVar (i.getItem(), v);
896  if (hasAlgVar)
897  {
898  if (b.getp() != 0)
899  {
900  CFList result= diophantineQa (F, G, factors, b, v);
901  return result;
902  }
903  CFList result= modularDiophant (F, factors, getMipo (v));
904  return result;
905  }
906  if (b.getp() != 0)
907  return diophantineHensel (F, factors, b);
908  }
909 
910  CanonicalForm buf1, buf2, buf3, S, T;
911  CFListIterator i= factors;
912  CFList result;
913  if (i.hasItem())
914  i++;
915  buf1= F/factors.getFirst();
916  buf2= divNTL (F, i.getItem());
917  buf3= extgcd (buf1, buf2, S, T);
918  result.append (S);
919  result.append (T);
920  if (i.hasItem())
921  i++;
922  for (; i.hasItem(); i++)
923  {
924  buf1= divNTL (F, i.getItem());
925  buf3= extgcd (buf3, buf1, S, T);
926  CFListIterator k= factors;
927  for (CFListIterator j= result; j.hasItem(); j++, k++)
928  {
929  j.getItem()= mulNTL (j.getItem(), S);
930  j.getItem()= modNTL (j.getItem(), k.getItem());
931  }
932  result.append (T);
933  }
934  return result;
935 }
936 
937 CFList
938 diophantine (const CanonicalForm& F, const CFList& factors)
939 {
940  modpk b= modpk();
941  return diophantine (F, 1, factors, b);
942 }
943 
944 void
945 henselStep12 (const CanonicalForm& F, const CFList& factors,
946  CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
947  CFArray& Pi, int j, const modpk& b)
948 {
950  CanonicalForm xToJ= power (F.mvar(), j);
951  Variable x= F.mvar();
952  // compute the error
953  if (j == 1)
954  E= F[j];
955  else
956  {
957  if (degree (Pi [factors.length() - 2], x) > 0)
958  E= F[j] - Pi [factors.length() - 2] [j];
959  else
960  E= F[j];
961  }
962 
963  if (b.getp() != 0)
964  E= b(E);
965  CFArray buf= CFArray (diophant.length());
966  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
967  int k= 0;
968  CanonicalForm remainder;
969  // actual lifting
970  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
971  {
972  if (degree (bufFactors[k], x) > 0)
973  {
974  if (k > 0)
975  remainder= modNTL (E, bufFactors[k] [0], b); //TODO precompute inverses of bufFactors[k][0]
976  else
977  remainder= E;
978  }
979  else
980  remainder= modNTL (E, bufFactors[k], b);
981 
982  buf[k]= mulNTL (i.getItem(), remainder, b);
983  if (degree (bufFactors[k], x) > 0)
984  buf[k]= modNTL (buf[k], bufFactors[k] [0], b);
985  else
986  buf[k]= modNTL (buf[k], bufFactors[k], b);
987  }
988  for (k= 1; k < factors.length(); k++)
989  {
990  bufFactors[k] += xToJ*buf[k];
991  if (b.getp() != 0)
992  bufFactors[k]= b(bufFactors[k]);
993  }
994 
995  // update Pi [0]
996  int degBuf0= degree (bufFactors[0], x);
997  int degBuf1= degree (bufFactors[1], x);
998  if (degBuf0 > 0 && degBuf1 > 0)
999  M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j], b);
1000  CanonicalForm uIZeroJ;
1001 
1002  if (degBuf0 > 0 && degBuf1 > 0)
1003  uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
1004  (bufFactors[1] [0] + buf[1]), b) - M(1, 1) - M(j + 1, 1);
1005  else if (degBuf0 > 0)
1006  uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1], b);
1007  else if (degBuf1 > 0)
1008  uIZeroJ= mulNTL (bufFactors[0], buf[1], b);
1009  else
1010  uIZeroJ= 0;
1011  if (b.getp() != 0)
1012  uIZeroJ= b (uIZeroJ);
1013  Pi [0] += xToJ*uIZeroJ;
1014  if (b.getp() != 0)
1015  Pi [0]= b (Pi[0]);
1016 
1017  CFArray tmp= CFArray (factors.length() - 1);
1018  for (k= 0; k < factors.length() - 1; k++)
1019  tmp[k]= 0;
1020  CFIterator one, two;
1021  one= bufFactors [0];
1022  two= bufFactors [1];
1023  if (degBuf0 > 0 && degBuf1 > 0)
1024  {
1025  for (k= 1; k <= (int) ceil (j/2.0); k++)
1026  {
1027  if (k != j - k + 1)
1028  {
1029  if ((one.hasTerms() && one.exp() == j - k + 1)
1030  && (two.hasTerms() && two.exp() == j - k + 1))
1031  {
1032  tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), (bufFactors[1][k]+
1033  two.coeff()), b) - M (k + 1, 1) - M (j - k + 2, 1);
1034  one++;
1035  two++;
1036  }
1037  else if (one.hasTerms() && one.exp() == j - k + 1)
1038  {
1039  tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1][k], b)
1040  - M (k + 1, 1);
1041  one++;
1042  }
1043  else if (two.hasTerms() && two.exp() == j - k + 1)
1044  {
1045  tmp[0] += mulNTL (bufFactors[0][k], (bufFactors[1][k]+two.coeff()), b)
1046  - M (k + 1, 1);
1047  two++;
1048  }
1049  }
1050  else
1051  {
1052  tmp[0] += M (k + 1, 1);
1053  }
1054  }
1055  }
1056  if (b.getp() != 0)
1057  tmp[0]= b (tmp[0]);
1058  Pi [0] += tmp[0]*xToJ*F.mvar();
1059 
1060  // update Pi [l]
1061  int degPi, degBuf;
1062  for (int l= 1; l < factors.length() - 1; l++)
1063  {
1064  degPi= degree (Pi [l - 1], x);
1065  degBuf= degree (bufFactors[l + 1], x);
1066  if (degPi > 0 && degBuf > 0)
1067  M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j], b);
1068  if (j == 1)
1069  {
1070  if (degPi > 0 && degBuf > 0)
1071  Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
1072  bufFactors[l + 1] [0] + buf[l + 1], b) - M (j + 1, l +1) -
1073  M (1, l + 1));
1074  else if (degPi > 0)
1075  Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1], b));
1076  else if (degBuf > 0)
1077  Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1], b));
1078  }
1079  else
1080  {
1081  if (degPi > 0 && degBuf > 0)
1082  {
1083  uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0], b);
1084  uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1], b);
1085  }
1086  else if (degPi > 0)
1087  uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1], b);
1088  else if (degBuf > 0)
1089  {
1090  uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0], b);
1091  uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1], b);
1092  }
1093  Pi[l] += xToJ*uIZeroJ;
1094  }
1095  one= bufFactors [l + 1];
1096  two= Pi [l - 1];
1097  if (two.hasTerms() && two.exp() == j + 1)
1098  {
1099  if (degBuf > 0 && degPi > 0)
1100  {
1101  tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0], b);
1102  two++;
1103  }
1104  else if (degPi > 0)
1105  {
1106  tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1], b);
1107  two++;
1108  }
1109  }
1110  if (degBuf > 0 && degPi > 0)
1111  {
1112  for (k= 1; k <= (int) ceil (j/2.0); k++)
1113  {
1114  if (k != j - k + 1)
1115  {
1116  if ((one.hasTerms() && one.exp() == j - k + 1) &&
1117  (two.hasTerms() && two.exp() == j - k + 1))
1118  {
1119  tmp[l] += mulNTL ((bufFactors[l+1][k] + one.coeff()), (Pi[l-1][k] +
1120  two.coeff()),b) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
1121  one++;
1122  two++;
1123  }
1124  else if (one.hasTerms() && one.exp() == j - k + 1)
1125  {
1126  tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()), Pi[l-1][k], b) -
1127  M (k + 1, l + 1);
1128  one++;
1129  }
1130  else if (two.hasTerms() && two.exp() == j - k + 1)
1131  {
1132  tmp[l] += mulNTL (bufFactors[l+1][k], (Pi[l-1][k] + two.coeff()), b)
1133  - M (k + 1, l + 1);
1134  two++;
1135  }
1136  }
1137  else
1138  tmp[l] += M (k + 1, l + 1);
1139  }
1140  }
1141  if (b.getp() != 0)
1142  tmp[l]= b (tmp[l]);
1143  Pi[l] += tmp[l]*xToJ*F.mvar();
1144  }
1145 }
1146 
1147 void
1148 henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1149  CFList& diophant, CFMatrix& M, modpk& b, bool sort)
1150 {
1151  if (sort)
1152  sortList (factors, Variable (1));
1153  Pi= CFArray (factors.length() - 1);
1154  CFListIterator j= factors;
1155  diophant= diophantine (F[0], F, factors, b);
1156  CanonicalForm bufF= F;
1157  if (getCharacteristic() == 0 && b.getp() != 0)
1158  {
1159  Variable v;
1160  bool hasAlgVar= hasFirstAlgVar (F, v);
1161  for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++)
1162  hasAlgVar= hasFirstAlgVar (i.getItem(), v);
1163  Variable w;
1164  bool hasAlgVar2= false;
1165  for (CFListIterator i= diophant; i.hasItem() && !hasAlgVar2; i++)
1166  hasAlgVar2= hasFirstAlgVar (i.getItem(), w);
1167  if (hasAlgVar && hasAlgVar2 && v!=w)
1168  {
1169  bufF= replacevar (bufF, v, w);
1170  for (CFListIterator i= factors; i.hasItem(); i++)
1171  i.getItem()= replacevar (i.getItem(), v, w);
1172  }
1173  }
1174 
1175  DEBOUTLN (cerr, "diophant= " << diophant);
1176  j++;
1177  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()), b);
1178  M (1, 1)= Pi [0];
1179  int i= 1;
1180  if (j.hasItem())
1181  j++;
1182  for (; j.hasItem(); j++, i++)
1183  {
1184  Pi [i]= mulNTL (Pi [i - 1], j.getItem(), b);
1185  M (1, i + 1)= Pi [i];
1186  }
1187  CFArray bufFactors= CFArray (factors.length());
1188  i= 0;
1189  for (CFListIterator k= factors; k.hasItem(); i++, k++)
1190  {
1191  if (i == 0)
1192  bufFactors[i]= mod (k.getItem(), F.mvar());
1193  else
1194  bufFactors[i]= k.getItem();
1195  }
1196  for (i= 1; i < l; i++)
1197  henselStep12 (bufF, factors, bufFactors, diophant, M, Pi, i, b);
1198 
1199  CFListIterator k= factors;
1200  for (i= 0; i < factors.length (); i++, k++)
1201  k.getItem()= bufFactors[i];
1202  factors.removeFirst();
1203 }
1204 
1205 void
1206 henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1207  CFList& diophant, CFMatrix& M, bool sort)
1208 {
1209  modpk dummy= modpk();
1210  henselLift12 (F, factors, l, Pi, diophant, M, dummy, sort);
1211 }
1212 
1213 void
1214 henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
1215  end, CFArray& Pi, const CFList& diophant, CFMatrix& M,
1216  const modpk& b)
1217 {
1218  CFArray bufFactors= CFArray (factors.length());
1219  int i= 0;
1220  CanonicalForm xToStart= power (F.mvar(), start);
1221  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1222  {
1223  if (i == 0)
1224  bufFactors[i]= mod (k.getItem(), xToStart);
1225  else
1226  bufFactors[i]= k.getItem();
1227  }
1228  for (i= start; i < end; i++)
1229  henselStep12 (F, factors, bufFactors, diophant, M, Pi, i, b);
1230 
1231  CFListIterator k= factors;
1232  for (i= 0; i < factors.length(); k++, i++)
1233  k.getItem()= bufFactors [i];
1234  factors.removeFirst();
1235  return;
1236 }
1237 
1238 CFList
1239 biDiophantine (const CanonicalForm& F, const CFList& factors, int d)
1240 {
1241  Variable y= F.mvar();
1242  CFList result;
1243  if (y.level() == 1)
1244  {
1245  result= diophantine (F, factors);
1246  return result;
1247  }
1248  else
1249  {
1250  CFList buf= factors;
1251  for (CFListIterator i= buf; i.hasItem(); i++)
1252  i.getItem()= mod (i.getItem(), y);
1253  CanonicalForm A= mod (F, y);
1254  int bufD= 1;
1255  CFList recResult= biDiophantine (A, buf, bufD);
1256  CanonicalForm e= 1;
1257  CFList p;
1258  CFArray bufFactors= CFArray (factors.length());
1259  CanonicalForm yToD= power (y, d);
1260  int k= 0;
1261  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1262  {
1263  bufFactors [k]= i.getItem();
1264  }
1265  CanonicalForm b, quot;
1266  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1267  {
1268  b= 1;
1269  if (fdivides (bufFactors[k], F, quot))
1270  b= quot;
1271  else
1272  {
1273  for (int l= 0; l < factors.length(); l++)
1274  {
1275  if (l == k)
1276  continue;
1277  else
1278  {
1279  b= mulMod2 (b, bufFactors[l], yToD);
1280  }
1281  }
1282  }
1283  p.append (b);
1284  }
1285 
1286  CFListIterator j= p;
1287  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1288  e -= i.getItem()*j.getItem();
1289 
1290  if (e.isZero())
1291  return recResult;
1292  CanonicalForm coeffE;
1293  CFList s;
1294  result= recResult;
1295  CanonicalForm g;
1296  for (int i= 1; i < d; i++)
1297  {
1298  if (degree (e, y) > 0)
1299  coeffE= e[i];
1300  else
1301  coeffE= 0;
1302  if (!coeffE.isZero())
1303  {
1305  CFListIterator l= p;
1306  int ii= 0;
1307  j= recResult;
1308  for (; j.hasItem(); j++, k++, l++, ii++)
1309  {
1310  g= coeffE*j.getItem();
1311  if (degree (bufFactors[ii], y) <= 0)
1312  g= mod (g, bufFactors[ii]);
1313  else
1314  g= mod (g, bufFactors[ii][0]);
1315  k.getItem() += g*power (y, i);
1316  e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
1317  DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
1318  mod (e, power (y, i + 1)));
1319  }
1320  }
1321  if (e.isZero())
1322  break;
1323  }
1324 
1325  DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
1326 
1327 #ifdef DEBUGOUTPUT
1328  CanonicalForm test= 0;
1329  j= p;
1330  for (CFListIterator i= result; i.hasItem(); i++, j++)
1331  test += mod (i.getItem()*j.getItem(), power (y, d));
1332  DEBOUTLN (cerr, "test= " << test);
1333 #endif
1334  return result;
1335  }
1336 }
1337 
1338 CFList
1339 multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
1340  const CFList& recResult, const CFList& M, int d)
1341 {
1342  Variable y= F.mvar();
1343  CFList result;
1344  CFListIterator i;
1345  CanonicalForm e= 1;
1346  CFListIterator j= factors;
1347  CFList p;
1348  CFArray bufFactors= CFArray (factors.length());
1349  CanonicalForm yToD= power (y, d);
1350  int k= 0;
1351  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1352  bufFactors [k]= i.getItem();
1353  CanonicalForm b, quot;
1354  CFList buf= M;
1355  buf.removeLast();
1356  buf.append (yToD);
1357  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1358  {
1359  b= 1;
1360  if (fdivides (bufFactors[k], F, quot))
1361  b= quot;
1362  else
1363  {
1364  for (int l= 0; l < factors.length(); l++)
1365  {
1366  if (l == k)
1367  continue;
1368  else
1369  {
1370  b= mulMod (b, bufFactors[l], buf);
1371  }
1372  }
1373  }
1374  p.append (b);
1375  }
1376  j= p;
1377  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1378  e -= mulMod (i.getItem(), j.getItem(), M);
1379 
1380  if (e.isZero())
1381  return recResult;
1382  CanonicalForm coeffE;
1383  CFList s;
1384  result= recResult;
1385  CanonicalForm g;
1386  for (int i= 1; i < d; i++)
1387  {
1388  if (degree (e, y) > 0)
1389  coeffE= e[i];
1390  else
1391  coeffE= 0;
1392  if (!coeffE.isZero())
1393  {
1395  CFListIterator l= p;
1396  j= recResult;
1397  int ii= 0;
1398  CanonicalForm dummy;
1399  for (; j.hasItem(); j++, k++, l++, ii++)
1400  {
1401  g= mulMod (coeffE, j.getItem(), M);
1402  if (degree (bufFactors[ii], y) <= 0)
1403  divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
1404  g, M);
1405  else
1406  divrem (g, bufFactors[ii][0], dummy, g, M);
1407  k.getItem() += g*power (y, i);
1408  e -= mulMod (g*power (y, i), l.getItem(), M);
1409  }
1410  }
1411 
1412  if (e.isZero())
1413  break;
1414  }
1415 
1416 #ifdef DEBUGOUTPUT
1417  CanonicalForm test= 0;
1418  j= p;
1419  for (CFListIterator i= result; i.hasItem(); i++, j++)
1420  test += mod (i.getItem()*j.getItem(), power (y, d));
1421  DEBOUTLN (cerr, "test in multiRecDiophantine= " << test);
1422 #endif
1423  return result;
1424 }
1425 
1426 static
1427 void
1428 henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
1429  const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
1430  const CFList& MOD)
1431 {
1432  CanonicalForm E;
1433  CanonicalForm xToJ= power (F.mvar(), j);
1434  Variable x= F.mvar();
1435  // compute the error
1436  if (j == 1)
1437  {
1438  E= F[j];
1439 #ifdef DEBUGOUTPUT
1440  CanonicalForm test= 1;
1441  for (int i= 0; i < factors.length(); i++)
1442  {
1443  if (i == 0)
1444  test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
1445  else
1446  test= mulMod (test, bufFactors[i], MOD);
1447  }
1448  CanonicalForm test2= mod (F-test, xToJ);
1449 
1450  test2= mod (test2, MOD);
1451  DEBOUTLN (cerr, "test in henselStep= " << test2);
1452 #endif
1453  }
1454  else
1455  {
1456 #ifdef DEBUGOUTPUT
1457  CanonicalForm test= 1;
1458  for (int i= 0; i < factors.length(); i++)
1459  {
1460  if (i == 0)
1461  test *= mod (bufFactors [i], power (x, j));
1462  else
1463  test *= bufFactors[i];
1464  }
1465  test= mod (test, power (x, j));
1466  test= mod (test, MOD);
1467  CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
1468  DEBOUTLN (cerr, "test in henselStep= " << test2);
1469 #endif
1470 
1471  if (degree (Pi [factors.length() - 2], x) > 0)
1472  E= F[j] - Pi [factors.length() - 2] [j];
1473  else
1474  E= F[j];
1475  }
1476 
1477  CFArray buf= CFArray (diophant.length());
1478  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1479  int k= 0;
1480  // actual lifting
1481  CanonicalForm dummy, rest1;
1482  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1483  {
1484  if (degree (bufFactors[k], x) > 0)
1485  {
1486  if (k > 0)
1487  divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
1488  else
1489  rest1= E;
1490  }
1491  else
1492  divrem (E, bufFactors[k], dummy, rest1, MOD);
1493 
1494  buf[k]= mulMod (i.getItem(), rest1, MOD);
1495 
1496  if (degree (bufFactors[k], x) > 0)
1497  divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
1498  else
1499  divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
1500  }
1501  for (k= 1; k < factors.length(); k++)
1502  bufFactors[k] += xToJ*buf[k];
1503 
1504  // update Pi [0]
1505  int degBuf0= degree (bufFactors[0], x);
1506  int degBuf1= degree (bufFactors[1], x);
1507  if (degBuf0 > 0 && degBuf1 > 0)
1508  M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
1509  CanonicalForm uIZeroJ;
1510 
1511  if (degBuf0 > 0 && degBuf1 > 0)
1512  uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
1513  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
1514  else if (degBuf0 > 0)
1515  uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
1516  else if (degBuf1 > 0)
1517  uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
1518  else
1519  uIZeroJ= 0;
1520  Pi [0] += xToJ*uIZeroJ;
1521 
1522  CFArray tmp= CFArray (factors.length() - 1);
1523  for (k= 0; k < factors.length() - 1; k++)
1524  tmp[k]= 0;
1525  CFIterator one, two;
1526  one= bufFactors [0];
1527  two= bufFactors [1];
1528  if (degBuf0 > 0 && degBuf1 > 0)
1529  {
1530  for (k= 1; k <= (int) ceil (j/2.0); k++)
1531  {
1532  if (k != j - k + 1)
1533  {
1534  if ((one.hasTerms() && one.exp() == j - k + 1) &&
1535  (two.hasTerms() && two.exp() == j - k + 1))
1536  {
1537  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
1538  (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
1539  M (j - k + 2, 1);
1540  one++;
1541  two++;
1542  }
1543  else if (one.hasTerms() && one.exp() == j - k + 1)
1544  {
1545  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
1546  bufFactors[1] [k], MOD) - M (k + 1, 1);
1547  one++;
1548  }
1549  else if (two.hasTerms() && two.exp() == j - k + 1)
1550  {
1551  tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
1552  two.coeff()), MOD) - M (k + 1, 1);
1553  two++;
1554  }
1555  }
1556  else
1557  {
1558  tmp[0] += M (k + 1, 1);
1559  }
1560  }
1561  }
1562  Pi [0] += tmp[0]*xToJ*F.mvar();
1563 
1564  // update Pi [l]
1565  int degPi, degBuf;
1566  for (int l= 1; l < factors.length() - 1; l++)
1567  {
1568  degPi= degree (Pi [l - 1], x);
1569  degBuf= degree (bufFactors[l + 1], x);
1570  if (degPi > 0 && degBuf > 0)
1571  M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
1572  if (j == 1)
1573  {
1574  if (degPi > 0 && degBuf > 0)
1575  Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
1576  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
1577  M (1, l + 1));
1578  else if (degPi > 0)
1579  Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
1580  else if (degBuf > 0)
1581  Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
1582  }
1583  else
1584  {
1585  if (degPi > 0 && degBuf > 0)
1586  {
1587  uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
1588  uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
1589  }
1590  else if (degPi > 0)
1591  uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
1592  else if (degBuf > 0)
1593  {
1594  uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
1595  uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
1596  }
1597  Pi[l] += xToJ*uIZeroJ;
1598  }
1599  one= bufFactors [l + 1];
1600  two= Pi [l - 1];
1601  if (two.hasTerms() && two.exp() == j + 1)
1602  {
1603  if (degBuf > 0 && degPi > 0)
1604  {
1605  tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
1606  two++;
1607  }
1608  else if (degPi > 0)
1609  {
1610  tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
1611  two++;
1612  }
1613  }
1614  if (degBuf > 0 && degPi > 0)
1615  {
1616  for (k= 1; k <= (int) ceil (j/2.0); k++)
1617  {
1618  if (k != j - k + 1)
1619  {
1620  if ((one.hasTerms() && one.exp() == j - k + 1) &&
1621  (two.hasTerms() && two.exp() == j - k + 1))
1622  {
1623  tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1624  (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
1625  M (j - k + 2, l + 1);
1626  one++;
1627  two++;
1628  }
1629  else if (one.hasTerms() && one.exp() == j - k + 1)
1630  {
1631  tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1632  Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
1633  one++;
1634  }
1635  else if (two.hasTerms() && two.exp() == j - k + 1)
1636  {
1637  tmp[l] += mulMod (bufFactors[l + 1] [k],
1638  (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
1639  two++;
1640  }
1641  }
1642  else
1643  tmp[l] += M (k + 1, l + 1);
1644  }
1645  }
1646  Pi[l] += tmp[l]*xToJ*F.mvar();
1647  }
1648 
1649  return;
1650 }
1651 
1652 CFList
1653 henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList&
1654  diophant, CFArray& Pi, CFMatrix& M)
1655 {
1656  CFList buf= factors;
1657  int k= 0;
1658  int liftBoundBivar= l[k];
1659  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
1660  CFList MOD;
1661  MOD.append (power (Variable (2), liftBoundBivar));
1662  CFArray bufFactors= CFArray (factors.length());
1663  k= 0;
1664  CFListIterator j= eval;
1665  j++;
1666  buf.removeFirst();
1667  buf.insert (LC (j.getItem(), 1));
1668  for (CFListIterator i= buf; i.hasItem(); i++, k++)
1669  bufFactors[k]= i.getItem();
1670  Pi= CFArray (factors.length() - 1);
1671  CFListIterator i= buf;
1672  i++;
1673  Variable y= j.getItem().mvar();
1674  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
1675  M (1, 1)= Pi [0];
1676  k= 1;
1677  if (i.hasItem())
1678  i++;
1679  for (; i.hasItem(); i++, k++)
1680  {
1681  Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
1682  M (1, k + 1)= Pi [k];
1683  }
1684 
1685  for (int d= 1; d < l[1]; d++)
1686  henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
1687  CFList result;
1688  for (k= 1; k < factors.length(); k++)
1689  result.append (bufFactors[k]);
1690  return result;
1691 }
1692 
1693 void
1694 henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
1695  CFArray& Pi, const CFList& diophant, CFMatrix& M,
1696  const CFList& MOD)
1697 {
1698  CFArray bufFactors= CFArray (factors.length());
1699  int i= 0;
1700  CanonicalForm xToStart= power (F.mvar(), start);
1701  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1702  {
1703  if (i == 0)
1704  bufFactors[i]= mod (k.getItem(), xToStart);
1705  else
1706  bufFactors[i]= k.getItem();
1707  }
1708  for (i= start; i < end; i++)
1709  henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD);
1710 
1711  CFListIterator k= factors;
1712  for (i= 0; i < factors.length(); k++, i++)
1713  k.getItem()= bufFactors [i];
1714  factors.removeFirst();
1715  return;
1716 }
1717 
1718 CFList
1719 henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
1720  diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew)
1721 {
1722  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
1723  int k= 0;
1724  CFArray bufFactors= CFArray (factors.length());
1725  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1726  {
1727  if (k == 0)
1728  bufFactors[k]= LC (F.getLast(), 1);
1729  else
1730  bufFactors[k]= i.getItem();
1731  }
1732  CFList buf= factors;
1733  buf.removeFirst();
1734  buf.insert (LC (F.getLast(), 1));
1735  CFListIterator i= buf;
1736  i++;
1737  Variable y= F.getLast().mvar();
1738  Variable x= F.getFirst().mvar();
1739  CanonicalForm xToLOld= power (x, lOld);
1740  Pi [0]= mod (Pi[0], xToLOld);
1741  M (1, 1)= Pi [0];
1742  k= 1;
1743  if (i.hasItem())
1744  i++;
1745  for (; i.hasItem(); i++, k++)
1746  {
1747  Pi [k]= mod (Pi [k], xToLOld);
1748  M (1, k + 1)= Pi [k];
1749  }
1750 
1751  for (int d= 1; d < lNew; d++)
1752  henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
1753  CFList result;
1754  for (k= 1; k < factors.length(); k++)
1755  result.append (bufFactors[k]);
1756  return result;
1757 }
1758 
1759 CFList
1760 henselLift (const CFList& eval, const CFList& factors, int* l, int lLength,
1761  bool sort)
1762 {
1763  CFList diophant;
1764  CFList buf= factors;
1765  buf.insert (LC (eval.getFirst(), 1));
1766  if (sort)
1767  sortList (buf, Variable (1));
1768  CFArray Pi;
1769  CFMatrix M= CFMatrix (l[1], factors.length());
1770  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
1771  if (eval.length() == 2)
1772  return result;
1773  CFList MOD;
1774  for (int i= 0; i < 2; i++)
1775  MOD.append (power (Variable (i + 2), l[i]));
1776  CFListIterator j= eval;
1777  j++;
1778  CFList bufEval;
1779  bufEval.append (j.getItem());
1780  j++;
1781 
1782  for (int i= 2; i < lLength && j.hasItem(); i++, j++)
1783  {
1784  result.insert (LC (bufEval.getFirst(), 1));
1785  bufEval.append (j.getItem());
1786  M= CFMatrix (l[i], factors.length());
1787  result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
1788  MOD.append (power (Variable (i + 2), l[i]));
1789  bufEval.removeFirst();
1790  }
1791  return result;
1792 }
1793 
1794 // nonmonic
1795 
1796 void
1797 nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors,
1798  CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1799  CFArray& Pi, int j, const CFArray& /*LCs*/)
1800 {
1801  Variable x= F.mvar();
1802  CanonicalForm xToJ= power (x, j);
1803 
1804  CanonicalForm E;
1805  // compute the error
1806  if (degree (Pi [factors.length() - 2], x) > 0)
1807  E= F[j] - Pi [factors.length() - 2] [j];
1808  else
1809  E= F[j];
1810 
1811  CFArray buf= CFArray (diophant.length());
1812 
1813  int k= 0;
1814  CanonicalForm remainder;
1815  // actual lifting
1816  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1817  {
1818  if (degree (bufFactors[k], x) > 0)
1819  remainder= modNTL (E, bufFactors[k] [0]);
1820  else
1821  remainder= modNTL (E, bufFactors[k]);
1822  buf[k]= mulNTL (i.getItem(), remainder);
1823  if (degree (bufFactors[k], x) > 0)
1824  buf[k]= modNTL (buf[k], bufFactors[k] [0]);
1825  else
1826  buf[k]= modNTL (buf[k], bufFactors[k]);
1827  }
1828 
1829  for (k= 0; k < factors.length(); k++)
1830  bufFactors[k] += xToJ*buf[k];
1831 
1832  // update Pi [0]
1833  int degBuf0= degree (bufFactors[0], x);
1834  int degBuf1= degree (bufFactors[1], x);
1835  if (degBuf0 > 0 && degBuf1 > 0)
1836  {
1837  M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
1838  if (j + 2 <= M.rows())
1839  M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
1840  }
1841  else
1842  M (j + 1, 1)= 0;
1843 
1844  CanonicalForm uIZeroJ;
1845  if (degBuf0 > 0 && degBuf1 > 0)
1846  uIZeroJ= mulNTL(bufFactors[0][0], buf[1]) +
1847  mulNTL (bufFactors[1][0], buf[0]);
1848  else if (degBuf0 > 0)
1849  uIZeroJ= mulNTL (buf[0], bufFactors[1]) +
1850  mulNTL (buf[1], bufFactors[0][0]);
1851  else if (degBuf1 > 0)
1852  uIZeroJ= mulNTL (bufFactors[0], buf[1]) +
1853  mulNTL (buf[0], bufFactors[1][0]);
1854  else
1855  uIZeroJ= mulNTL (bufFactors[0], buf[1]) +
1856  mulNTL (buf[0], bufFactors[1]);
1857 
1858  Pi [0] += xToJ*uIZeroJ;
1859 
1860  CFArray tmp= CFArray (factors.length() - 1);
1861  for (k= 0; k < factors.length() - 1; k++)
1862  tmp[k]= 0;
1863  CFIterator one, two;
1864  one= bufFactors [0];
1865  two= bufFactors [1];
1866  if (degBuf0 > 0 && degBuf1 > 0)
1867  {
1868  while (one.hasTerms() && one.exp() > j) one++;
1869  while (two.hasTerms() && two.exp() > j) two++;
1870  for (k= 1; k <= (int) ceil (j/2.0); k++)
1871  {
1872  if (k != j - k + 1)
1873  {
1874  if ((one.hasTerms() && one.exp() == j - k + 1) && +
1875  (two.hasTerms() && two.exp() == j - k + 1))
1876  {
1877  tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
1878  two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
1879  one++;
1880  two++;
1881  }
1882  else if (one.hasTerms() && one.exp() == j - k + 1)
1883  {
1884  tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
1885  M (k + 1, 1);
1886  one++;
1887  }
1888  else if (two.hasTerms() && two.exp() == j - k + 1)
1889  {
1890  tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
1891  M (k + 1, 1);
1892  two++;
1893  }
1894  }
1895  else
1896  tmp[0] += M (k + 1, 1);
1897  }
1898  }
1899 
1900  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
1901  {
1902  if (j + 2 <= M.rows())
1903  tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
1904  (bufFactors [1] [j + 1] + bufFactors [1] [0]))
1905  - M(1,1) - M (j + 2,1);
1906  }
1907  else if (degBuf0 >= j + 1)
1908  {
1909  if (degBuf1 > 0)
1910  tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
1911  else
1912  tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]);
1913  }
1914  else if (degBuf1 >= j + 1)
1915  {
1916  if (degBuf0 > 0)
1917  tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
1918  else
1919  tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]);
1920  }
1921 
1922  Pi [0] += tmp[0]*xToJ*F.mvar();
1923 
1924  int degPi, degBuf;
1925  for (int l= 1; l < factors.length() - 1; l++)
1926  {
1927  degPi= degree (Pi [l - 1], x);
1928  degBuf= degree (bufFactors[l + 1], x);
1929  if (degPi > 0 && degBuf > 0)
1930  {
1931  M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
1932  if (j + 2 <= M.rows())
1933  M (j + 2, l + 1)= mulNTL (Pi [l - 1][j + 1], bufFactors[l + 1] [j + 1]);
1934  }
1935  else
1936  M (j + 1, l + 1)= 0;
1937 
1938  if (degPi > 0 && degBuf > 0)
1939  uIZeroJ= mulNTL (Pi[l - 1] [0], buf[l + 1]) +
1940  mulNTL (uIZeroJ, bufFactors[l+1] [0]);
1941  else if (degPi > 0)
1942  uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]) +
1943  mulNTL (Pi[l - 1][0], buf[l + 1]);
1944  else if (degBuf > 0)
1945  uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1][0]) +
1946  mulNTL (Pi[l - 1], buf[l + 1]);
1947  else
1948  uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]) +
1949  mulNTL (Pi[l - 1], buf[l + 1]);
1950 
1951  Pi [l] += xToJ*uIZeroJ;
1952 
1953  one= bufFactors [l + 1];
1954  two= Pi [l - 1];
1955  if (degBuf > 0 && degPi > 0)
1956  {
1957  while (one.hasTerms() && one.exp() > j) one++;
1958  while (two.hasTerms() && two.exp() > j) two++;
1959  for (k= 1; k <= (int) ceil (j/2.0); k++)
1960  {
1961  if (k != j - k + 1)
1962  {
1963  if ((one.hasTerms() && one.exp() == j - k + 1) &&
1964  (two.hasTerms() && two.exp() == j - k + 1))
1965  {
1966  tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
1967  (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1) -
1968  M (j - k + 2, l + 1);
1969  one++;
1970  two++;
1971  }
1972  else if (one.hasTerms() && one.exp() == j - k + 1)
1973  {
1974  tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
1975  Pi[l - 1] [k]) - M (k + 1, l + 1);
1976  one++;
1977  }
1978  else if (two.hasTerms() && two.exp() == j - k + 1)
1979  {
1980  tmp[l] += mulNTL (bufFactors[l + 1] [k],
1981  (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1);
1982  two++;
1983  }
1984  }
1985  else
1986  tmp[l] += M (k + 1, l + 1);
1987  }
1988  }
1989 
1990  if (degPi >= j + 1 && degBuf >= j + 1)
1991  {
1992  if (j + 2 <= M.rows())
1993  tmp [l] += mulNTL ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
1994  (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
1995  ) - M(1,l+1) - M (j + 2,l+1);
1996  }
1997  else if (degPi >= j + 1)
1998  {
1999  if (degBuf > 0)
2000  tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1] [0]);
2001  else
2002  tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1]);
2003  }
2004  else if (degBuf >= j + 1)
2005  {
2006  if (degPi > 0)
2007  tmp[l] += mulNTL (Pi [l - 1] [0], bufFactors [l + 1] [j + 1]);
2008  else
2009  tmp[l] += mulNTL (Pi [l - 1], bufFactors [l + 1] [j + 1]);
2010  }
2011 
2012  Pi[l] += tmp[l]*xToJ*F.mvar();
2013  }
2014  return;
2015 }
2016 
2017 void
2018 nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l,
2019  CFArray& Pi, CFList& diophant, CFMatrix& M,
2020  const CFArray& LCs, bool sort)
2021 {
2022  if (sort)
2023  sortList (factors, Variable (1));
2024  Pi= CFArray (factors.length() - 2);
2025  CFList bufFactors2= factors;
2026  bufFactors2.removeFirst();
2027  diophant= diophantine (F[0], bufFactors2);
2028  DEBOUTLN (cerr, "diophant= " << diophant);
2029 
2030  CFArray bufFactors= CFArray (bufFactors2.length());
2031  int i= 0;
2032  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
2033  bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
2034 
2035  Variable x= F.mvar();
2036  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
2037  {
2038  M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
2039  Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
2040  mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
2041  }
2042  else if (degree (bufFactors[0], x) > 0)
2043  {
2044  M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
2045  Pi [0]= M (1, 1) +
2046  mulNTL (bufFactors [0] [1], bufFactors[1])*x;
2047  }
2048  else if (degree (bufFactors[1], x) > 0)
2049  {
2050  M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
2051  Pi [0]= M (1, 1) +
2052  mulNTL (bufFactors [0], bufFactors[1] [1])*x;
2053  }
2054  else
2055  {
2056  M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
2057  Pi [0]= M (1, 1);
2058  }
2059 
2060  for (i= 1; i < Pi.size(); i++)
2061  {
2062  if (degree (Pi[i-1], x) > 0 && degree (bufFactors [i+1], x) > 0)
2063  {
2064  M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors[i+1] [0]);
2065  Pi [i]= M (1,i+1) + (mulNTL (Pi[i-1] [1], bufFactors[i+1] [0]) +
2066  mulNTL (Pi[i-1] [0], bufFactors [i+1] [1]))*x;
2067  }
2068  else if (degree (Pi[i-1], x) > 0)
2069  {
2070  M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors [i+1]);
2071  Pi [i]= M(1,i+1) + mulNTL (Pi[i-1] [1], bufFactors[i+1])*x;
2072  }
2073  else if (degree (bufFactors[i+1], x) > 0)
2074  {
2075  M (1,i+1)= mulNTL (Pi[i-1], bufFactors [i+1] [0]);
2076  Pi [i]= M (1,i+1) + mulNTL (Pi[i-1], bufFactors[i+1] [1])*x;
2077  }
2078  else
2079  {
2080  M (1,i+1)= mulNTL (Pi [i-1], bufFactors [i+1]);
2081  Pi [i]= M (1,i+1);
2082  }
2083  }
2084 
2085  for (i= 1; i < l; i++)
2086  nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
2087 
2088  factors= CFList();
2089  for (i= 0; i < bufFactors.size(); i++)
2090  factors.append (bufFactors[i]);
2091  return;
2092 }
2093 
2094 
2095 /// solve \f$ E=\sum_{i= 1}^r{\sigma_{i}\prod_{j=1, j\neq i}^{r}{f_{j}}} \f$
2096 /// mod M, @a products contains \f$ \prod_{j=1, j\neq i}^{r}{f_{j}} \f$
2097 CFList
2098 diophantine (const CFList& recResult, const CFList& factors,
2099  const CFList& products, const CFList& M, const CanonicalForm& E,
2100  bool& bad)
2101 {
2102  if (M.isEmpty())
2103  {
2104  CFList result;
2105  CFListIterator j= factors;
2107  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2108  {
2109  ASSERT (E.isUnivariate() || E.inCoeffDomain(),
2110  "constant or univariate poly expected");
2111  ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
2112  "constant or univariate poly expected");
2113  ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
2114  "constant or univariate poly expected");
2115  buf= mulNTL (E, i.getItem());
2116  result.append (modNTL (buf, j.getItem()));
2117  }
2118  return result;
2119  }
2120  Variable y= M.getLast().mvar();
2121  CFList bufFactors= factors;
2122  for (CFListIterator i= bufFactors; i.hasItem(); i++)
2123  i.getItem()= mod (i.getItem(), y);
2124  CFList bufProducts= products;
2125  for (CFListIterator i= bufProducts; i.hasItem(); i++)
2126  i.getItem()= mod (i.getItem(), y);
2127  CFList buf= M;
2128  buf.removeLast();
2129  CanonicalForm bufE= mod (E, y);
2130  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2131  bufE, bad);
2132 
2133  if (bad)
2134  return CFList();
2135 
2136  CanonicalForm e= E;
2137  CFListIterator j= products;
2138  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
2139  e -= j.getItem()*i.getItem();
2140 
2141  CFList result= recDiophantine;
2142  int d= degree (M.getLast());
2143  CanonicalForm coeffE;
2144  for (int i= 1; i < d; i++)
2145  {
2146  if (degree (e, y) > 0)
2147  coeffE= e[i];
2148  else
2149  coeffE= 0;
2150  if (!coeffE.isZero())
2151  {
2153  recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2154  coeffE, bad);
2155  if (bad)
2156  return CFList();
2157  CFListIterator l= products;
2158  for (j= recDiophantine; j.hasItem(); j++, k++, l++)
2159  {
2160  k.getItem() += j.getItem()*power (y, i);
2161  e -= l.getItem()*(j.getItem()*power (y, i));
2162  }
2163  }
2164  if (e.isZero())
2165  break;
2166  }
2167  if (!e.isZero())
2168  {
2169  bad= true;
2170  return CFList();
2171  }
2172  return result;
2173 }
2174 
2175 void
2176 nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
2177  CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
2178  CFArray& Pi, const CFList& products, int j,
2179  const CFList& MOD, bool& noOneToOne)
2180 {
2181  CanonicalForm E;
2182  CanonicalForm xToJ= power (F.mvar(), j);
2183  Variable x= F.mvar();
2184 
2185  // compute the error
2186 #ifdef DEBUGOUTPUT
2187  CanonicalForm test= 1;
2188  for (int i= 0; i < factors.length(); i++)
2189  {
2190  if (i == 0)
2191  test *= mod (bufFactors [i], power (x, j));
2192  else
2193  test *= bufFactors[i];
2194  }
2195  test= mod (test, power (x, j));
2196  test= mod (test, MOD);
2197  CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
2198  DEBOUTLN (cerr, "test in nonMonicHenselStep= " << test2);
2199 #endif
2200 
2201  if (degree (Pi [factors.length() - 2], x) > 0)
2202  E= F[j] - Pi [factors.length() - 2] [j];
2203  else
2204  E= F[j];
2205 
2206  CFArray buf= CFArray (diophant.length());
2207 
2208  // actual lifting
2209  TIMING_START (diotime);
2210  CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
2211  noOneToOne);
2212 
2213  if (noOneToOne)
2214  return;
2215 
2216  int k= 0;
2217  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
2218  {
2219  buf[k]= i.getItem();
2220  bufFactors[k] += xToJ*i.getItem();
2221  }
2222  TIMING_END_AND_PRINT (diotime, "time for dio: ");
2223 
2224  // update Pi [0]
2225  TIMING_START (product2);
2226  int degBuf0= degree (bufFactors[0], x);
2227  int degBuf1= degree (bufFactors[1], x);
2228  if (degBuf0 > 0 && degBuf1 > 0)
2229  {
2230  M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
2231  if (j + 2 <= M.rows())
2232  M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
2233  }
2234  else
2235  M (j + 1, 1)= 0;
2236 
2237  CanonicalForm uIZeroJ;
2238  if (degBuf0 > 0 && degBuf1 > 0)
2239  uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
2240  mulMod (bufFactors[1] [0], buf[0], MOD);
2241  else if (degBuf0 > 0)
2242  uIZeroJ= mulMod (buf[0], bufFactors[1], MOD) +
2243  mulMod (buf[1], bufFactors[0][0], MOD);
2244  else if (degBuf1 > 0)
2245  uIZeroJ= mulMod (bufFactors[0], buf[1], MOD) +
2246  mulMod (buf[0], bufFactors[1][0], MOD);
2247  else
2248  uIZeroJ= mulMod (bufFactors[0], buf[1], MOD) +
2249  mulMod (buf[0], bufFactors[1], MOD);
2250  Pi [0] += xToJ*uIZeroJ;
2251 
2252  CFArray tmp= CFArray (factors.length() - 1);
2253  for (k= 0; k < factors.length() - 1; k++)
2254  tmp[k]= 0;
2255  CFIterator one, two;
2256  one= bufFactors [0];
2257  two= bufFactors [1];
2258  if (degBuf0 > 0 && degBuf1 > 0)
2259  {
2260  while (one.hasTerms() && one.exp() > j) one++;
2261  while (two.hasTerms() && two.exp() > j) two++;
2262  for (k= 1; k <= (int) ceil (j/2.0); k++)
2263  {
2264  if (k != j - k + 1)
2265  {
2266  if ((one.hasTerms() && one.exp() == j - k + 1) &&
2267  (two.hasTerms() && two.exp() == j - k + 1))
2268  {
2269  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2270  (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
2271  M (j - k + 2, 1);
2272  one++;
2273  two++;
2274  }
2275  else if (one.hasTerms() && one.exp() == j - k + 1)
2276  {
2277  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2278  bufFactors[1] [k], MOD) - M (k + 1, 1);
2279  one++;
2280  }
2281  else if (two.hasTerms() && two.exp() == j - k + 1)
2282  {
2283  tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
2284  two.coeff()), MOD) - M (k + 1, 1);
2285  two++;
2286  }
2287  }
2288  else
2289  {
2290  tmp[0] += M (k + 1, 1);
2291  }
2292  }
2293  }
2294 
2295  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
2296  {
2297  if (j + 2 <= M.rows())
2298  tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
2299  (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
2300  - M(1,1) - M (j + 2,1);
2301  }
2302  else if (degBuf0 >= j + 1)
2303  {
2304  if (degBuf1 > 0)
2305  tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
2306  else
2307  tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
2308  }
2309  else if (degBuf1 >= j + 1)
2310  {
2311  if (degBuf0 > 0)
2312  tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
2313  else
2314  tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
2315  }
2316  Pi [0] += tmp[0]*xToJ*F.mvar();
2317 
2318  // update Pi [l]
2319  int degPi, degBuf;
2320  for (int l= 1; l < factors.length() - 1; l++)
2321  {
2322  degPi= degree (Pi [l - 1], x);
2323  degBuf= degree (bufFactors[l + 1], x);
2324  if (degPi > 0 && degBuf > 0)
2325  {
2326  M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
2327  if (j + 2 <= M.rows())
2328  M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
2329  MOD);
2330  }
2331  else
2332  M (j + 1, l + 1)= 0;
2333 
2334  if (degPi > 0 && degBuf > 0)
2335  uIZeroJ= mulMod (Pi[l - 1] [0], buf[l + 1], MOD) +
2336  mulMod (uIZeroJ, bufFactors[l + 1] [0], MOD);
2337  else if (degPi > 0)
2338  uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD) +
2339  mulMod (Pi[l - 1][0], buf[l + 1], MOD);
2340  else if (degBuf > 0)
2341  uIZeroJ= mulMod (Pi[l - 1], buf[l + 1], MOD) +
2342  mulMod (uIZeroJ, bufFactors[l + 1][0], MOD);
2343  else
2344  uIZeroJ= mulMod (Pi[l - 1], buf[l + 1], MOD) +
2345  mulMod (uIZeroJ, bufFactors[l + 1], MOD);
2346 
2347  Pi [l] += xToJ*uIZeroJ;
2348 
2349  one= bufFactors [l + 1];
2350  two= Pi [l - 1];
2351  if (degBuf > 0 && degPi > 0)
2352  {
2353  while (one.hasTerms() && one.exp() > j) one++;
2354  while (two.hasTerms() && two.exp() > j) two++;
2355  for (k= 1; k <= (int) ceil (j/2.0); k++)
2356  {
2357  if (k != j - k + 1)
2358  {
2359  if ((one.hasTerms() && one.exp() == j - k + 1) &&
2360  (two.hasTerms() && two.exp() == j - k + 1))
2361  {
2362  tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2363  (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
2364  M (j - k + 2, l + 1);
2365  one++;
2366  two++;
2367  }
2368  else if (one.hasTerms() && one.exp() == j - k + 1)
2369  {
2370  tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2371  Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
2372  one++;
2373  }
2374  else if (two.hasTerms() && two.exp() == j - k + 1)
2375  {
2376  tmp[l] += mulMod (bufFactors[l + 1] [k],
2377  (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
2378  two++;
2379  }
2380  }
2381  else
2382  tmp[l] += M (k + 1, l + 1);
2383  }
2384  }
2385 
2386  if (degPi >= j + 1 && degBuf >= j + 1)
2387  {
2388  if (j + 2 <= M.rows())
2389  tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
2390  (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
2391  , MOD) - M(1,l+1) - M (j + 2,l+1);
2392  }
2393  else if (degPi >= j + 1)
2394  {
2395  if (degBuf > 0)
2396  tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
2397  else
2398  tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
2399  }
2400  else if (degBuf >= j + 1)
2401  {
2402  if (degPi > 0)
2403  tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
2404  else
2405  tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
2406  }
2407 
2408  Pi[l] += tmp[l]*xToJ*F.mvar();
2409  }
2410  TIMING_END_AND_PRINT (product2, "time for product in hensel step: ");
2411  return;
2412 }
2413 
2414 // wrt. Variable (1)
2416 {
2417  if (degree (F, 1) <= 0)
2418  return c;
2419  else
2420  {
2421  CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
2422  result += (swapvar (c, Variable (F.level() + 1), Variable (1))
2423  - LC (result))*power (result.mvar(), degree (result));
2424  return swapvar (result, Variable (F.level() + 1), Variable (1));
2425  }
2426 }
2427 
2428 CFList
2429 nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
2430  diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
2431  const CFList& LCs2, bool& bad)
2432 {
2433  CFList buf= factors;
2434  int k= 0;
2435  int liftBoundBivar= l[k];
2436  CFList bufbuf= factors;
2437  Variable v= Variable (2);
2438 
2439  CFList MOD;
2440  MOD.append (power (Variable (2), liftBoundBivar));
2441  CFArray bufFactors= CFArray (factors.length());
2442  k= 0;
2443  CFListIterator j= eval;
2444  j++;
2445  CFListIterator iter1= LCs1;
2446  CFListIterator iter2= LCs2;
2447  iter1++;
2448  iter2++;
2449  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
2450  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
2451 
2452  CFListIterator i= buf;
2453  i++;
2454  Variable y= j.getItem().mvar();
2455  if (y.level() != 3)
2456  y= Variable (3);
2457 
2458  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
2459  M (1, 1)= Pi[0];
2460  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2461  Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2462  mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2463  else if (degree (bufFactors[0], y) > 0)
2464  Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2465  else if (degree (bufFactors[1], y) > 0)
2466  Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2467 
2468  CFList products;
2469  for (int i= 0; i < bufFactors.size(); i++)
2470  {
2471  if (degree (bufFactors[i], y) > 0)
2472  products.append (eval.getFirst()/bufFactors[i] [0]);
2473  else
2474  products.append (eval.getFirst()/bufFactors[i]);
2475  }
2476 
2477  for (int d= 1; d < l[1]; d++)
2478  {
2479  nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
2480  d, MOD, bad);
2481  if (bad)
2482  return CFList();
2483  }
2484  CFList result;
2485  for (k= 0; k < factors.length(); k++)
2486  result.append (bufFactors[k]);
2487  return result;
2488 }
2489 
2490 
2491 CFList
2492 nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
2493  CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2494  int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
2495  )
2496 {
2497  int k= 0;
2498  CFArray bufFactors= CFArray (factors.length());
2499  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
2500  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
2501  CFList buf= factors;
2502  Variable y= F.getLast().mvar();
2503  Variable x= F.getFirst().mvar();
2504  CanonicalForm xToLOld= power (x, lOld);
2505  Pi [0]= mod (Pi[0], xToLOld);
2506  M (1, 1)= Pi [0];
2507 
2508  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2509  Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2510  mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2511  else if (degree (bufFactors[0], y) > 0)
2512  Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2513  else if (degree (bufFactors[1], y) > 0)
2514  Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2515 
2516  CFList products;
2517  CanonicalForm quot;
2518  for (int i= 0; i < bufFactors.size(); i++)
2519  {
2520  if (degree (bufFactors[i], y) > 0)
2521  {
2522  if (!fdivides (bufFactors[i] [0], F.getFirst(), quot))
2523  {
2524  bad= true;
2525  return CFList();
2526  }
2527  products.append (quot);
2528  }
2529  else
2530  {
2531  if (!fdivides (bufFactors[i], F.getFirst(), quot))
2532  {
2533  bad= true;
2534  return CFList();
2535  }
2536  products.append (quot);
2537  }
2538  }
2539 
2540  for (int d= 1; d < lNew; d++)
2541  {
2542  nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
2543  d, MOD, bad);
2544  if (bad)
2545  return CFList();
2546  }
2547 
2548  CFList result;
2549  for (k= 0; k < factors.length(); k++)
2550  result.append (bufFactors[k]);
2551  return result;
2552 }
2553 
2554 CFList
2555 nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
2556  lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
2557  const CFArray& Pi, const CFList& diophant, bool& bad)
2558 {
2559  CFList bufDiophant= diophant;
2560  CFList buf= factors;
2561  if (sort)
2562  sortList (buf, Variable (1));
2563  CFArray bufPi= Pi;
2564  CFMatrix M= CFMatrix (l[1], factors.length());
2565  CFList result=
2566  nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
2567  if (bad)
2568  return CFList();
2569 
2570  if (eval.length() == 2)
2571  return result;
2572  CFList MOD;
2573  for (int i= 0; i < 2; i++)
2574  MOD.append (power (Variable (i + 2), l[i]));
2575  CFListIterator j= eval;
2576  j++;
2577  CFList bufEval;
2578  bufEval.append (j.getItem());
2579  j++;
2580  CFListIterator jj= LCs1;
2581  CFListIterator jjj= LCs2;
2582  CFList bufLCs1, bufLCs2;
2583  jj++, jjj++;
2584  bufLCs1.append (jj.getItem());
2585  bufLCs2.append (jjj.getItem());
2586  jj++, jjj++;
2587 
2588  for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
2589  {
2590  bufEval.append (j.getItem());
2591  bufLCs1.append (jj.getItem());
2592  bufLCs2.append (jjj.getItem());
2593  M= CFMatrix (l[i], factors.length());
2594  result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
2595  l[i - 1], l[i], bufLCs1, bufLCs2, bad);
2596  if (bad)
2597  return CFList();
2598  MOD.append (power (Variable (i + 2), l[i]));
2599  bufEval.removeFirst();
2600  bufLCs1.removeFirst();
2601  bufLCs2.removeFirst();
2602  }
2603  return result;
2604 }
2605 
2606 CFList
2607 nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
2608  CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
2609  int bivarLiftBound, bool& bad)
2610 {
2611  CFList bufFactors2= factors;
2612 
2613  Variable y= Variable (2);
2614  for (CFListIterator i= bufFactors2; i.hasItem(); i++)
2615  i.getItem()= mod (i.getItem(), y);
2616 
2617  CanonicalForm bufF= F;
2618  bufF= mod (bufF, y);
2619  bufF= mod (bufF, Variable (3));
2620 
2621  TIMING_START (diotime);
2622  diophant= diophantine (bufF, bufFactors2);
2623  TIMING_END_AND_PRINT (diotime, "time for dio in 23: ");
2624 
2625  CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
2626 
2627  Pi= CFArray (bufFactors2.length() - 1);
2628 
2629  CFArray bufFactors= CFArray (bufFactors2.length());
2630  CFListIterator j= LCs;
2631  int i= 0;
2632  for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
2633  bufFactors[i]= replaceLC (k.getItem(), j.getItem());
2634 
2635  //initialise Pi
2636  Variable v= Variable (3);
2637  CanonicalForm yToL= power (y, bivarLiftBound);
2638  if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
2639  {
2640  M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
2641  Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
2642  mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
2643  }
2644  else if (degree (bufFactors[0], v) > 0)
2645  {
2646  M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
2647  Pi [0]= M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
2648  }
2649  else if (degree (bufFactors[1], v) > 0)
2650  {
2651  M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
2652  Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
2653  }
2654  else
2655  {
2656  M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
2657  Pi [0]= M (1,1);
2658  }
2659 
2660  for (i= 1; i < Pi.size(); i++)
2661  {
2662  if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
2663  {
2664  M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
2665  Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
2666  mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
2667  }
2668  else if (degree (Pi[i-1], v) > 0)
2669  {
2670  M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
2671  Pi [i]= M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
2672  }
2673  else if (degree (bufFactors[i+1], v) > 0)
2674  {
2675  M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
2676  Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
2677  }
2678  else
2679  {
2680  M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
2681  Pi [i]= M (1,i+1);
2682  }
2683  }
2684 
2685  CFList products;
2686  bufF= mod (F, Variable (3));
2687  TIMING_START (product1);
2688  for (CFListIterator k= factors; k.hasItem(); k++)
2689  products.append (bufF/k.getItem());
2690  TIMING_END_AND_PRINT (product1, "time for product1 in 23: ");
2691 
2692  CFList MOD= CFList (power (v, liftBound));
2693  MOD.insert (yToL);
2694  for (int d= 1; d < liftBound; d++)
2695  {
2696  nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
2697  MOD, bad);
2698  if (bad)
2699  return CFList();
2700  }
2701 
2702  CFList result;
2703  for (i= 0; i < factors.length(); i++)
2704  result.append (bufFactors[i]);
2705  return result;
2706 }
2707 
2708 CFList
2709 nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
2710  CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2711  int& lNew, const CFList& MOD, bool& noOneToOne
2712  )
2713 {
2714 
2715  int k= 0;
2716  CFArray bufFactors= CFArray (factors.length());
2717  CFListIterator j= LCs;
2718  for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
2719  bufFactors [k]= replaceLC (i.getItem(), j.getItem());
2720 
2721  Variable y= F.getLast().mvar();
2722  Variable x= F.getFirst().mvar();
2723  CanonicalForm xToLOld= power (x, lOld);
2724 
2725  Pi [0]= mod (Pi[0], xToLOld);
2726  M (1, 1)= Pi [0];
2727 
2728  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2729  Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2730  mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2731  else if (degree (bufFactors[0], y) > 0)
2732  Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2733  else if (degree (bufFactors[1], y) > 0)
2734  Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2735 
2736  for (int i= 1; i < Pi.size(); i++)
2737  {
2738  Pi [i]= mod (Pi [i], xToLOld);
2739  M (1, i + 1)= Pi [i];
2740 
2741  if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
2742  Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
2743  mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
2744  else if (degree (Pi[i-1], y) > 0)
2745  Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
2746  else if (degree (bufFactors[i+1], y) > 0)
2747  Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
2748  }
2749 
2750  CFList products;
2751  CanonicalForm quot, bufF= F.getFirst();
2752 
2753  TIMING_START (product1);
2754  for (int i= 0; i < bufFactors.size(); i++)
2755  {
2756  if (degree (bufFactors[i], y) > 0)
2757  {
2758  if (!fdivides (bufFactors[i] [0], bufF, quot))
2759  {
2760  noOneToOne= true;
2761  return factors;
2762  }
2763  products.append (quot);
2764  }
2765  else
2766  {
2767  if (!fdivides (bufFactors[i], bufF, quot))
2768  {
2769  noOneToOne= true;
2770  return factors;
2771  }
2772  products.append (quot);
2773  }
2774  }
2775  TIMING_END_AND_PRINT (product1, "time for product1 in further hensel:" );
2776 
2777  for (int d= 1; d < lNew; d++)
2778  {
2779  nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
2780  products, d, MOD, noOneToOne);
2781  if (noOneToOne)
2782  return CFList();
2783  }
2784 
2785  CFList result;
2786  for (k= 0; k < factors.length(); k++)
2787  result.append (bufFactors[k]);
2788  return result;
2789 }
2790 
2791 CFList
2792 nonMonicHenselLift (const CFList& eval, const CFList& factors,
2793  CFList* const& LCs, CFList& diophant, CFArray& Pi,
2794  int* liftBound, int length, bool& noOneToOne
2795  )
2796 {
2797  CFList bufDiophant= diophant;
2798  CFList buf= factors;
2799  CFArray bufPi= Pi;
2800  CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
2801  int k= 0;
2802 
2803  TIMING_START (hensel23);
2804  CFList result=
2805  nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
2806  liftBound[1], liftBound[0], noOneToOne);
2807  TIMING_END_AND_PRINT (hensel23, "time for 23: ");
2808 
2809  if (noOneToOne)
2810  return CFList();
2811 
2812  if (eval.length() == 1)
2813  return result;
2814 
2815  k++;
2816  CFList MOD;
2817  for (int i= 0; i < 2; i++)
2818  MOD.append (power (Variable (i + 2), liftBound[i]));
2819 
2820  CFListIterator j= eval;
2821  CFList bufEval;
2822  bufEval.append (j.getItem());
2823  j++;
2824 
2825  for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
2826  {
2827  bufEval.append (j.getItem());
2828  M= CFMatrix (liftBound[i], factors.length() - 1);
2829  TIMING_START (hensel);
2830  result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
2831  liftBound[i-1], liftBound[i], MOD, noOneToOne);
2832  TIMING_END_AND_PRINT (hensel, "time for further hensel: ");
2833  if (noOneToOne)
2834  return result;
2835  MOD.append (power (Variable (i + 2), liftBound[i]));
2836  bufEval.removeFirst();
2837  }
2838 
2839  return result;
2840 }
2841 
2842 #endif
2843 /* HAVE_NTL */
2844 
int status int void size_t count
Definition: si_signals.h:59
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
List< CanonicalForm > CFList
int cf_getNumBigPrimes()
Definition: cf_primes.cc:45
nmod_poly_init(FLINTmipo, getCharacteristic())
const CanonicalForm int s
Definition: facAbsFact.cc:55
int level() const
Definition: variable.h:49
ZZ convertFacCF2NTLZZ(const CanonicalForm &f)
NAME: convertFacCF2NTLZZX.
Definition: NTLconvert.cc:666
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
int j
Definition: facHensel.cc:105
GCD over Q(a)
Conversion to and from NTL.
This file defines functions for Hensel lifting.
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
const poly a
Definition: syzextra.cc:212
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
static CFList replacevar(const CFList &L, const Variable &a, const Variable &b)
Definition: facHensel.cc:281
void Off(int sw)
switches
bool inCoeffDomain() const
Matrix< CanonicalForm > CFMatrix
CanonicalForm buf1
Definition: facFqBivar.cc:71
functions to print debug output
CanonicalForm extgcd(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &a, CanonicalForm &b)
CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a...
Definition: cfUnivarGcd.cc:173
void setReduce(const Variable &alpha, bool reduce)
Definition: variable.cc:238
#define TIMING_END_AND_PRINT(t, msg)
Definition: timing.h:89
int getp() const
Definition: fac_util.h:35
CanonicalForm reduce(const CanonicalForm &f, const CanonicalForm &M)
polynomials in M.mvar() are considered coefficients M univariate monic polynomial the coefficients of...
Definition: cf_ops.cc:646
return P p
Definition: myNF.cc:203
f
Definition: cfModGcd.cc:4022
fq_nmod_t buf
Definition: facHensel.cc:96
factory&#39;s class for variables
Definition: variable.h:32
static void chineseRemainder(const CFList &x1, const CanonicalForm &q1, const CFList &x2, const CanonicalForm &q2, CFList &xnew, CanonicalForm &qnew)
Definition: facHensel.cc:256
int size() const
Definition: ftmpl_array.cc:92
bool bad
Definition: facFactorize.cc:65
void nonMonicHenselStep(const CanonicalForm &F, const CFList &factors, CFArray &bufFactors, const CFList &diophant, CFMatrix &M, CFArray &Pi, const CFList &products, int j, const CFList &MOD, bool &noOneToOne)
Definition: facHensel.cc:2176
void henselStep12(const CanonicalForm &F, const CFList &factors, CFArray &bufFactors, const CFList &diophant, CFMatrix &M, CFArray &Pi, int j, const modpk &b)
Definition: facHensel.cc:945
CanonicalForm replaceLc(const CanonicalForm &f, const CanonicalForm &c)
Definition: fac_util.cc:89
void convertFacCF2Fq_nmod_poly_t(fq_nmod_poly_t result, const CanonicalForm &f, const fq_nmod_ctx_t ctx)
conversion of a factory univariate poly over F_q to a FLINT fq_nmod_poly_t
ZZ_pEX convertFacCF2NTLZZ_pEX(const CanonicalForm &f, const ZZ_pX &mipo)
CanonicalForm in Z_p(a)[X] to NTL ZZ_pEX.
Definition: NTLconvert.cc:1038
factory&#39;s main class
Definition: canonicalform.h:75
nmod_poly_clear(FLINTmipo)
assertions for Factory
fq_nmod_ctx_t fq_con
Definition: facHensel.cc:94
CFList diophantineHensel(const CanonicalForm &F, const CFList &factors, const modpk &b)
Definition: facHensel.cc:472
Array< CanonicalForm > CFArray
CFList biDiophantine(const CanonicalForm &F, const CFList &factors, int d)
Definition: facHensel.cc:1239
g
Definition: cfModGcd.cc:4031
ZZX convertFacCF2NTLZZX(const CanonicalForm &f)
Definition: NTLconvert.cc:687
int k
Definition: cfEzgcd.cc:93
void nonMonicHenselStep12(const CanonicalForm &F, const CFList &factors, CFArray &bufFactors, const CFList &diophant, CFMatrix &M, CFArray &Pi, int j, const CFArray &)
Definition: facHensel.cc:1797
static void henselStep(const CanonicalForm &F, const CFList &factors, CFArray &bufFactors, const CFList &diophant, CFMatrix &M, CFArray &Pi, int j, const CFList &MOD)
Definition: facHensel.cc:1428
Variable alpha
Definition: facAbsBiFact.cc:52
static int mod(const CFList &L, const CanonicalForm &p)
Definition: facHensel.cc:244
CFList nonMonicHenselLift(const CFList &F, const CFList &factors, const CFList &LCs, CFList &diophant, CFArray &Pi, CFMatrix &M, int lOld, int &lNew, const CFList &MOD, bool &noOneToOne)
Definition: facHensel.cc:2709
CF_NO_INLINE bool isZero() const
Definition: cf_inline.cc:372
void insert(const T &)
Definition: ftmpl_list.cc:193
CanonicalForm convertFq_nmod_poly_t2FacCF(const fq_nmod_poly_t p, const Variable &x, const Variable &alpha, const fq_nmod_ctx_t ctx)
conversion of a FLINT poly over Fq to a CanonicalForm with alg. variable alpha and polynomial variabl...
static TreeM * G
Definition: janet.cc:38
CanonicalForm Lc(const CanonicalForm &f)
void tryInvert(const CanonicalForm &F, const CanonicalForm &M, CanonicalForm &inv, bool &fail)
Definition: cfGcdAlgExt.cc:220
void setCharacteristic(int c)
Definition: cf_char.cc:23
int getCharacteristic()
Definition: cf_char.cc:51
CanonicalForm inverse(const CanonicalForm &f, bool symmetric=true) const
Definition: fac_util.cc:58
void removeFirst()
Definition: ftmpl_list.cc:287
fq_nmod_ctx_clear(fq_con)
void sortList(CFList &list, const Variable &x)
sort a list of polynomials by their degree in x.
Definition: facHensel.cc:441
CF_NO_INLINE bool isOne() const
CF_INLINE bool CanonicalForm::isOne, isZero () const.
Definition: cf_inline.cc:354
int level() const
level() returns the level of CO.
TIMING_DEFINE_PRINT(diotime) TIMING_DEFINE_PRINT(product1) TIMING_DEFINE_PRINT(product2) TIMING_DEFINE_PRINT(hensel23) TIMING_DEFINE_PRINT(hensel) static CFList productsFLINT(const CFList &factors
CF_NO_INLINE int hasTerms() const
check if iterator has reached < the end of CanonicalForm
int isEmpty() const
Definition: ftmpl_list.cc:267
CanonicalForm lc(const CanonicalForm &f)
static CFList Farey(const CFList &L, const CanonicalForm &q)
Definition: facHensel.cc:272
CFList nonMonicHenselLift2(const CFList &F, const CFList &factors, const CFList &MOD, CFList &diophant, CFArray &Pi, CFMatrix &M, int lOld, int &lNew, const CFList &LCs1, const CFList &LCs2, bool &bad)
Definition: facHensel.cc:2492
return modpk(p, k)
void tryNTLXGCD(zz_pEX &d, zz_pEX &s, zz_pEX &t, const zz_pEX &a, const zz_pEX &b, bool &fail)
compute the extended GCD d=s*a+t*b, fail is set to true if a zero divisor is encountered ...
CF_NO_INLINE CanonicalForm coeff() const
get the current coefficient
bool equal
Definition: cfModGcd.cc:4067
CanonicalForm swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
fq_nmod_poly_t * vec
Definition: facHensel.cc:103
void removeLast()
Definition: ftmpl_list.cc:317
CanonicalForm mulNTL(const CanonicalForm &F, const CanonicalForm &G, const modpk &b)
multiplication of univariate polys using FLINT/NTL over F_p, F_q, Z/p^k, Z/p^k[t]/(f), Z, Q, Q(a), if we are in GF factory&#39;s default multiplication is used. If b!= 0 and getCharacteristic() == 0 the input will be considered as elements over Z/p^k or Z/p^k[t]/(f).
Definition: facMul.cc:407
CFList multiRecDiophantine(const CanonicalForm &F, const CFList &factors, const CFList &recResult, const CFList &M, int d)
Definition: facHensel.cc:1339
univariate Gcd over finite fields and Z, extended GCD over finite fields and Q
CanonicalForm convertNTLZZ_pEX2CF(const ZZ_pEX &f, const Variable &x, const Variable &alpha)
Definition: NTLconvert.cc:1116
CanonicalForm divNTL(const CanonicalForm &F, const CanonicalForm &G, const modpk &b)
division of univariate polys using FLINT/NTL over F_p, F_q, Z/p^k, Z/p^k[t]/(f), Z, Q, Q(a), if we are in GF factory&#39;s default multiplication is used. If b!= 0 and getCharacteristic() == 0 the input will be considered as elements over Z/p^k or Z/p^k[t]/(f); in this case invertiblity of Lc(G) is not checked
Definition: facMul.cc:850
convertFacCF2nmod_poly_t(FLINTmipo, M)
This file defines functions for fast multiplication and division with remainder.
CanonicalForm modNTL(const CanonicalForm &F, const CanonicalForm &G, const modpk &b)
mod of univariate polys using FLINT/NTL over F_p, F_q, Z/p^k, Z/p^k[t]/(f), Z, Q, Q(a)...
Definition: facMul.cc:676
zz_pEX convertFacCF2NTLzz_pEX(const CanonicalForm &f, const zz_pX &mipo)
Definition: NTLconvert.cc:1065
const CanonicalForm & M
Definition: facHensel.cc:92
#define A
Definition: sirandom.c:23
CFList henselLift23(const CFList &eval, const CFList &factors, int *l, CFList &diophant, CFArray &Pi, CFMatrix &M)
Hensel lifting from bivariate to trivariate.
Definition: facHensel.cc:1653
int getk() const
Definition: fac_util.h:36
CFList & eval
Definition: facFactorize.cc:48
CFList result
Definition: facHensel.cc:121
fq_nmod_ctx_init_modulus(fq_con, FLINTmipo,"Z")
void convertFacCF2Fq_nmod_t(fq_nmod_t result, const CanonicalForm &f, const fq_nmod_ctx_t ctx)
conversion of a factory element of F_q to a FLINT fq_nmod_t, does not do any memory allocation for po...
CFList nonMonicHenselLift232(const CFList &eval, const CFList &factors, int *l, CFList &diophant, CFArray &Pi, CFMatrix &M, const CFList &LCs1, const CFList &LCs2, bool &bad)
Definition: facHensel.cc:2429
bool hasFirstAlgVar(const CanonicalForm &f, Variable &a)
check if poly f contains an algebraic variable a
Definition: cf_ops.cc:665
CanonicalForm mulMod2(const CanonicalForm &A, const CanonicalForm &B, const CanonicalForm &M)
Karatsuba style modular multiplication for bivariate polynomials.
Definition: facMul.cc:2853
static const int SW_RATIONAL
set to 1 for computations over Q
Definition: cf_defs.h:28
int m
Definition: cfEzgcd.cc:119
CanonicalForm replaceLC(const CanonicalForm &F, const CanonicalForm &c)
Definition: facHensel.cc:2415
bool isOn(int sw)
switches
void On(int sw)
switches
T getFirst() const
Definition: ftmpl_list.cc:279
CanonicalForm mapinto() const
int length() const
Definition: ftmpl_list.cc:273
int i
Definition: cfEzgcd.cc:123
CanonicalForm getpk() const
Definition: fac_util.h:38
CFList tmp2
Definition: facFqBivar.cc:70
CFList diophantine(const CanonicalForm &F, const CFList &factors)
Definition: facHensel.cc:938
declarations of higher level algorithms.
This file defines functions for univariate GCD and extended GCD over Z/p[t]/(f)[x] for reducible f...
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
Variable x
Definition: facHensel.cc:122
fq_nmod_poly_clear(prod, fq_con)
CanonicalForm convertNTLzz_pEX2CF(const zz_pEX &f, const Variable &x, const Variable &alpha)
Definition: NTLconvert.cc:1093
int rows() const
Definition: ftmpl_matrix.h:45
bool isUnivariate() const
static CFList mapinto(const CFList &L)
Definition: facHensel.cc:235
CanonicalForm maxNorm(const CanonicalForm &f)
CanonicalForm maxNorm ( const CanonicalForm & f )
CFList diophantineQa(const CanonicalForm &F, const CanonicalForm &G, const CFList &factors, modpk &b, const Variable &alpha)
solve mod over by first computing mod and if no zero divisor occurred compute it mod ...
Definition: facHensel.cc:776
class to iterate through CanonicalForm&#39;s
Definition: cf_iter.h:44
CanonicalForm test
Definition: cfModGcd.cc:4037
CanonicalForm buf2
Definition: facFqBivar.cc:71
void insert(const T &)
Definition: ftmpl_list.cc:492
Variable beta
Definition: facAbsFact.cc:99
void henselLift12(const CanonicalForm &F, CFList &factors, int l, CFArray &Pi, CFList &diophant, CFMatrix &M, modpk &b, bool sort)
Hensel lift from univariate to bivariate.
Definition: facHensel.cc:1148
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
CF_NO_INLINE CanonicalForm div(const CanonicalForm &, const CanonicalForm &)
CF_INLINE CanonicalForm div, mod ( const CanonicalForm & lhs, const CanonicalForm & rhs ) ...
Definition: cf_inline.cc:553
CanonicalForm mulMod(const CanonicalForm &A, const CanonicalForm &B, const CFList &MOD)
Karatsuba style modular multiplication for multivariate polynomials.
Definition: facMul.cc:2947
access to prime tables
REvaluation E(1, terms.length(), IntRandom(25))
fq_nmod_poly_init(prod, fq_con)
CF_NO_INLINE int exp() const
get the current exponent
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:103
CanonicalForm den(const CanonicalForm &f)
T & getItem() const
Definition: ftmpl_list.cc:431
fq_nmod_poly_t prod
Definition: facHensel.cc:95
#define TIMING_START(t)
Definition: timing.h:87
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
CFList tmp1
Definition: facFqBivar.cc:70
const CanonicalForm & w
Definition: facAbsFact.cc:55
T getLast() const
Definition: ftmpl_list.cc:309
void henselLiftResume(const CanonicalForm &F, CFList &factors, int start, int end, CFArray &Pi, const CFList &diophant, CFMatrix &M, const CFList &MOD)
resume Hensel lifting.
Definition: facHensel.cc:1694
void findGoodPrime(const CanonicalForm &f, int &start)
find a big prime p from our tables such that no term of f vanishes mod p
Definition: facBivar.cc:60
void henselLiftResume12(const CanonicalForm &F, CFList &factors, int start, int end, CFArray &Pi, const CFList &diophant, CFMatrix &M, const modpk &b)
resume Hensel lift from univariate to bivariate. Assumes factors are lifted to precision Variable (2)...
Definition: facHensel.cc:1214
void sort(CFArray &A, int l=0)
quick sort A
CFList diophantineHenselQa(const CanonicalForm &F, const CanonicalForm &G, const CFList &factors, modpk &b, const Variable &alpha)
solve mod over by p-adic lifting
Definition: facHensel.cc:564
#define ASSERT(expression, message)
Definition: cf_assert.h:99
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
modpk coeffBound(const CanonicalForm &f, int p, const CanonicalForm &mipo)
compute p^k larger than the bound on the coefficients of a factor of f over Q (mipo) ...
Definition: facBivar.cc:96
operations mod p^k and some other useful functions for factorization
bivariate factorization over Q(a)
static void tryDiophantine(CFList &result, const CanonicalForm &F, const CFList &factors, const CanonicalForm &M, bool &fail)
Definition: facHensel.cc:148
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
void append(const T &)
Definition: ftmpl_list.cc:256
long fac_NTL_char
Definition: NTLconvert.cc:44
int degree(const CanonicalForm &f)
static jList * T
Definition: janet.cc:37
CanonicalForm LC(const CanonicalForm &f)
int hasAlgVar(const CanonicalForm &f, const Variable &v)
CFList modularDiophant(const CanonicalForm &f, const CFList &factors, const CanonicalForm &M)
Definition: facHensel.cc:290
CFList nonMonicHenselLift23(const CanonicalForm &F, const CFList &factors, const CFList &LCs, CFList &diophant, CFArray &Pi, int liftBound, int bivarLiftBound, bool &bad)
Definition: facHensel.cc:2607
const poly b
Definition: syzextra.cc:213
class to do operations mod p^k for int&#39;s p and k
Definition: fac_util.h:22
void nonMonicHenselLift12(const CanonicalForm &F, CFList &factors, int l, CFArray &Pi, CFList &diophant, CFMatrix &M, const CFArray &LCs, bool sort)
Hensel lifting from univariate to bivariate, factors need not to be monic.
Definition: facHensel.cc:2018
int l
Definition: cfEzgcd.cc:94
Variable rootOf(const CanonicalForm &mipo, char name)
returns a symbolic root of polynomial with name name Use it to define algebraic variables ...
Definition: variable.cc:162
void divrem(const CanonicalForm &f, const CanonicalForm &g, CanonicalForm &q, CanonicalForm &r)
CFList henselLift(const CFList &F, const CFList &factors, const CFList &MOD, CFList &diophant, CFArray &Pi, CFMatrix &M, int lOld, int lNew)
Hensel lifting.
Definition: facHensel.cc:1719