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