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 #ifndef HAVE_NTL
825 next_prime:
826 #endif
827  int i= 0;
828  while (cf_getBigPrime (i) <= p)
829  i++;
830  findGoodPrime (F, i);
831  findGoodPrime (G, i);
832  p=cf_getBigPrime(i);
833  b = coeffBound( G, p, mipo );
834  modpk bb= coeffBound (F, p, mipo );
835  if (bb.getk() > b.getk() ) b=bb;
836  fail= false;
837  }
838  else
839  break;
840  }
841  setReduce (alpha, true);
842  setCharacteristic (0);
843 
844  Variable gamma= alpha;
846  if (mipoHasDen)
847  {
848  On (SW_RATIONAL);
849  modMipo= getMipo (alpha);
850  den= bCommonDen (modMipo);
851  modMipo *= den;
852  Off (SW_RATIONAL);
853  setReduce (alpha, false);
854  gamma= rootOf (b (modMipo*b.inverse (den)));
855  setReduce (alpha, true);
856  }
857 
858  Variable x= Variable (1);
859  CanonicalForm buf1, buf2, buf3, S;
860  CFList bufFactors= factors;
861  CFListIterator i= bufFactors;
862  if (mipoHasDen)
863  {
864  for (; i.hasItem(); i++)
865  i.getItem()= replacevar (i.getItem(), alpha, gamma);
866  }
867  i= bufFactors;
868  CFList result;
869  if (i.hasItem())
870  i++;
871  buf1= 0;
872  CanonicalForm Freplaced;
873  if (mipoHasDen)
874  {
875  Freplaced= replacevar (F, alpha, gamma);
876  buf2= divNTL (Freplaced, replacevar (i.getItem(), alpha, gamma), b);
877  }
878  else
879  buf2= divNTL (F, i.getItem(), b);
880 
881 #ifdef HAVE_NTL
882  ZZ_p::init (convertFacCF2NTLZZ (b.getpk()));
883  ZZ_pX NTLmipo= to_ZZ_pX (convertFacCF2NTLZZX (getMipo (gamma)));
884  ZZ_pE::init (NTLmipo);
885  ZZ_pEX NTLS, NTLT, NTLbuf3;
886  ZZ_pEX NTLbuf1= convertFacCF2NTLZZ_pEX (buf1, NTLmipo);
887  ZZ_pEX NTLbuf2= convertFacCF2NTLZZ_pEX (buf2, NTLmipo);
888  XGCD (NTLbuf3, NTLS, NTLT, NTLbuf1, NTLbuf2);
889  result.append (b (convertNTLZZ_pEX2CF (NTLS, x, gamma)));
890  result.append (b (convertNTLZZ_pEX2CF (NTLT, x, gamma)));
891 #else // HAVE_FLINT
892 
893  fmpz_init(bigpk); // does convert expect an initalized object?
894  convertCF2initFmpz(bigpk, b.getpk());
895  fmpz_mod_poly_t FLINTmipo;
896  convertFacCF2Fmpz_mod_poly_t(FLINTmipo, getMipo(gamma), bigpk);
897 #if __FLINT_RELEASE >= 20700
898  fmpz_mod_ctx_t bigpk_ctx;
899  fmpz_mod_ctx_init(bigpk_ctx, bigpk);
900  fq_ctx_init_modulus(fqctx, FLINTmipo, bigpk_ctx, "Z");
901  fmpz_mod_ctx_clear(bigpk_ctx);
902  fmpz_mod_poly_clear(FLINTmipo, bigpk_ctx);
903 #else
904  fq_ctx_init_modulus(fqctx, FLINTmipo, "Z");
905  fmpz_mod_poly_clear(FLINTmipo);
906 #endif
907 
908  fq_init(fcheck, fqctx);
909  fq_poly_init(FLINTS, fqctx);
910  fq_poly_init(FLINTT, fqctx);
911  fq_poly_init(FLINTbuf3, fqctx);
912  //fq_poly_init(FLINTbuf1, fqctx); //convert expects uninitialized!
913  //fq_poly_init(FLINTbuf2, fqctx); //convert expects uninitialized!
914  convertFacCF2Fq_poly_t(FLINTbuf1, buf1, fqctx);
915  convertFacCF2Fq_poly_t(FLINTbuf2, buf2, fqctx);
916 
917  fq_poly_xgcd_euclidean_f(fcheck, FLINTbuf3, FLINTS, FLINTT,
918  FLINTbuf1, FLINTbuf2, fqctx);
919  if (!fq_is_one(fcheck, fqctx))
920  {
921  fmpz_clear(bigpk);
922  fq_clear(fcheck, fqctx);
923  fq_poly_clear(FLINTS, fqctx);
924  fq_poly_clear(FLINTT, fqctx);
925  fq_poly_clear(FLINTbuf3, fqctx);
926  fq_poly_clear(FLINTbuf1, fqctx);
927  fq_poly_clear(FLINTbuf2, fqctx);
928  fq_ctx_clear(fqctx);
929  setReduce (alpha, false);
930  fail = true;
931  goto next_prime;
932  }
933 
934  result.append(b(convertFq_poly_t2FacCF(FLINTS, x, alpha, fqctx)));
935  result.append(b(convertFq_poly_t2FacCF(FLINTT, x, alpha, fqctx)));
936 #endif
937 
938  if (i.hasItem())
939  i++;
940  for (; i.hasItem(); i++)
941  {
942  if (mipoHasDen)
943  buf1= divNTL (Freplaced, i.getItem(), b);
944  else
945  buf1= divNTL (F, i.getItem(), b);
946 
947 #ifdef HAVE_NTL
948  XGCD (NTLbuf3, NTLS, NTLT, NTLbuf3, convertFacCF2NTLZZ_pEX (buf1, NTLmipo));
949  S= convertNTLZZ_pEX2CF (NTLS, x, gamma);
950 #else
951  fq_poly_clear(FLINTbuf1, fqctx); //convert expects uninitialized!
952  convertFacCF2Fq_poly_t(FLINTbuf1, buf1, fqctx);
953 
954  // xgcd aliasing bug in <= 2.7.1
955  fq_poly_xgcd_euclidean_f(fcheck, FLINTbuf2, FLINTS, FLINTT,
956  FLINTbuf3, FLINTbuf1, fqctx);
957  fq_poly_swap(FLINTbuf3, FLINTbuf2, fqctx);
958 
959  if (!fq_is_one(fcheck, fqctx))
960  {
961  fmpz_clear(bigpk);
962  fq_clear(fcheck, fqctx);
963  fq_poly_clear(FLINTS, fqctx);
964  fq_poly_clear(FLINTT, fqctx);
965  fq_poly_clear(FLINTbuf3, fqctx);
966  fq_poly_clear(FLINTbuf1, fqctx);
967  fq_poly_clear(FLINTbuf2, fqctx);
968  fq_ctx_clear(fqctx);
969  setReduce (alpha, false);
970  fail = true;
971  goto next_prime;
972  }
973 
974  S= convertFq_poly_t2FacCF(FLINTS, x, alpha, fqctx);
975 #endif
976 
977  CFListIterator k= bufFactors;
978  for (CFListIterator j= result; j.hasItem(); j++, k++)
979  {
980  j.getItem()= mulNTL (j.getItem(), S, b);
981  j.getItem()= modNTL (j.getItem(), k.getItem(), b);
982  }
983 #if HAVE_NTL
984  result.append (b (convertNTLZZ_pEX2CF (NTLT, x, gamma)));
985 #else
986  result.append (b (convertFq_poly_t2FacCF(FLINTT, x, alpha, fqctx)));
987 #endif
988  }
989 
990 #if HAVE_NTL
991  // no cleanup for ntl
992 #else
993  fmpz_clear(bigpk);
994  fq_clear(fcheck, fqctx);
995  fq_poly_clear(FLINTS, fqctx);
996  fq_poly_clear(FLINTT, fqctx);
997  fq_poly_clear(FLINTbuf3, fqctx);
998  fq_poly_clear(FLINTbuf1, fqctx);
999  fq_poly_clear(FLINTbuf2, fqctx);
1000  fq_ctx_clear(fqctx);
1001 #endif
1002 
1003  return result;
1004 }
1005 #endif
1006 
1007 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // diophantineQa
1008 CFList
1010  const CFList& factors, modpk& b)
1011 {
1012  if (getCharacteristic() == 0)
1013  {
1014  Variable v;
1015  bool hasAlgVar= hasFirstAlgVar (F, v);
1016  for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++)
1017  hasAlgVar= hasFirstAlgVar (i.getItem(), v);
1018  if (hasAlgVar)
1019  {
1020  if (b.getp() != 0)
1021  {
1022  CFList result= diophantineQa (F, G, factors, b, v);
1023  return result;
1024  }
1025  CFList result= modularDiophant (F, factors, getMipo (v));
1026  return result;
1027  }
1028  if (b.getp() != 0)
1029  return diophantineHensel (F, factors, b);
1030  }
1031 
1032  CanonicalForm buf1, buf2, buf3, S, T;
1033  CFListIterator i= factors;
1034  CFList result;
1035  if (i.hasItem())
1036  i++;
1037  buf1= F/factors.getFirst();
1038  buf2= divNTL (F, i.getItem());
1039  buf3= extgcd (buf1, buf2, S, T);
1040  result.append (S);
1041  result.append (T);
1042  if (i.hasItem())
1043  i++;
1044  for (; i.hasItem(); i++)
1045  {
1046  buf1= divNTL (F, i.getItem());
1047  buf3= extgcd (buf3, buf1, S, T);
1048  CFListIterator k= factors;
1049  for (CFListIterator j= result; j.hasItem(); j++, k++)
1050  {
1051  j.getItem()= mulNTL (j.getItem(), S);
1052  j.getItem()= modNTL (j.getItem(), k.getItem());
1053  }
1054  result.append (T);
1055  }
1056  return result;
1057 }
1058 #endif
1059 
1060 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // diophantineQa
1061 CFList
1062 diophantine (const CanonicalForm& F, const CFList& factors)
1063 {
1064  modpk b= modpk();
1065  return diophantine (F, 1, factors, b);
1066 }
1067 #endif
1068 
1069 void
1070 henselStep12 (const CanonicalForm& F, const CFList& factors,
1071  CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1072  CFArray& Pi, int j, const modpk& b)
1073 {
1074  CanonicalForm E;
1075  CanonicalForm xToJ= power (F.mvar(), j);
1076  Variable x= F.mvar();
1077  // compute the error
1078  if (j == 1)
1079  E= F[j];
1080  else
1081  {
1082  if (degree (Pi [factors.length() - 2], x) > 0)
1083  E= F[j] - Pi [factors.length() - 2] [j];
1084  else
1085  E= F[j];
1086  }
1087 
1088  if (b.getp() != 0)
1089  E= b(E);
1090  CFArray buf= CFArray (diophant.length());
1091  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1092  int k= 0;
1094  // actual lifting
1095  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1096  {
1097  if (degree (bufFactors[k], x) > 0)
1098  {
1099  if (k > 0)
1100  remainder= modNTL (E, bufFactors[k] [0], b); //TODO precompute inverses of bufFactors[k][0]
1101  else
1102  remainder= E;
1103  }
1104  else
1105  remainder= modNTL (E, bufFactors[k], b);
1106 
1107  buf[k]= mulNTL (i.getItem(), remainder, b);
1108  if (degree (bufFactors[k], x) > 0)
1109  buf[k]= modNTL (buf[k], bufFactors[k] [0], b);
1110  else
1111  buf[k]= modNTL (buf[k], bufFactors[k], b);
1112  }
1113  for (k= 1; k < factors.length(); k++)
1114  {
1115  bufFactors[k] += xToJ*buf[k];
1116  if (b.getp() != 0)
1117  bufFactors[k]= b(bufFactors[k]);
1118  }
1119 
1120  // update Pi [0]
1121  int degBuf0= degree (bufFactors[0], x);
1122  int degBuf1= degree (bufFactors[1], x);
1123  if (degBuf0 > 0 && degBuf1 > 0)
1124  M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j], b);
1125  CanonicalForm uIZeroJ;
1126 
1127  if (degBuf0 > 0 && degBuf1 > 0)
1128  uIZeroJ= mulNTL ((bufFactors[0] [0] + bufFactors[0] [j]),
1129  (bufFactors[1] [0] + buf[1]), b) - M(1, 1) - M(j + 1, 1);
1130  else if (degBuf0 > 0)
1131  uIZeroJ= mulNTL (bufFactors[0] [j], bufFactors[1], b);
1132  else if (degBuf1 > 0)
1133  uIZeroJ= mulNTL (bufFactors[0], buf[1], b);
1134  else
1135  uIZeroJ= 0;
1136  if (b.getp() != 0)
1137  uIZeroJ= b (uIZeroJ);
1138  Pi [0] += xToJ*uIZeroJ;
1139  if (b.getp() != 0)
1140  Pi [0]= b (Pi[0]);
1141 
1142  CFArray tmp= CFArray (factors.length() - 1);
1143  for (k= 0; k < factors.length() - 1; k++)
1144  tmp[k]= 0;
1145  CFIterator one, two;
1146  one= bufFactors [0];
1147  two= bufFactors [1];
1148  if (degBuf0 > 0 && degBuf1 > 0)
1149  {
1150  for (k= 1; k <= (j+1)/2; k++)
1151  {
1152  if (k != j - k + 1)
1153  {
1154  if ((one.hasTerms() && one.exp() == j - k + 1)
1155  && (two.hasTerms() && two.exp() == j - k + 1))
1156  {
1157  tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), (bufFactors[1][k]+
1158  two.coeff()), b) - M (k + 1, 1) - M (j - k + 2, 1);
1159  one++;
1160  two++;
1161  }
1162  else if (one.hasTerms() && one.exp() == j - k + 1)
1163  {
1164  tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1][k], b)
1165  - M (k + 1, 1);
1166  one++;
1167  }
1168  else if (two.hasTerms() && two.exp() == j - k + 1)
1169  {
1170  tmp[0] += mulNTL (bufFactors[0][k], (bufFactors[1][k]+two.coeff()), b)
1171  - M (k + 1, 1);
1172  two++;
1173  }
1174  }
1175  else
1176  {
1177  tmp[0] += M (k + 1, 1);
1178  }
1179  }
1180  }
1181  if (b.getp() != 0)
1182  tmp[0]= b (tmp[0]);
1183  Pi [0] += tmp[0]*xToJ*F.mvar();
1184 
1185  // update Pi [l]
1186  int degPi, degBuf;
1187  for (int l= 1; l < factors.length() - 1; l++)
1188  {
1189  degPi= degree (Pi [l - 1], x);
1190  degBuf= degree (bufFactors[l + 1], x);
1191  if (degPi > 0 && degBuf > 0)
1192  M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j], b);
1193  if (j == 1)
1194  {
1195  if (degPi > 0 && degBuf > 0)
1196  Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
1197  bufFactors[l + 1] [0] + buf[l + 1], b) - M (j + 1, l +1) -
1198  M (1, l + 1));
1199  else if (degPi > 0)
1200  Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1], b));
1201  else if (degBuf > 0)
1202  Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1], b));
1203  }
1204  else
1205  {
1206  if (degPi > 0 && degBuf > 0)
1207  {
1208  uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0], b);
1209  uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1], b);
1210  }
1211  else if (degPi > 0)
1212  uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1], b);
1213  else if (degBuf > 0)
1214  {
1215  uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0], b);
1216  uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1], b);
1217  }
1218  Pi[l] += xToJ*uIZeroJ;
1219  }
1220  one= bufFactors [l + 1];
1221  two= Pi [l - 1];
1222  if (two.hasTerms() && two.exp() == j + 1)
1223  {
1224  if (degBuf > 0 && degPi > 0)
1225  {
1226  tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0], b);
1227  two++;
1228  }
1229  else if (degPi > 0)
1230  {
1231  tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1], b);
1232  two++;
1233  }
1234  }
1235  if (degBuf > 0 && degPi > 0)
1236  {
1237  for (k= 1; k <= (j+1)/2; k++)
1238  {
1239  if (k != j - k + 1)
1240  {
1241  if ((one.hasTerms() && one.exp() == j - k + 1) &&
1242  (two.hasTerms() && two.exp() == j - k + 1))
1243  {
1244  tmp[l] += mulNTL ((bufFactors[l+1][k] + one.coeff()), (Pi[l-1][k] +
1245  two.coeff()),b) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
1246  one++;
1247  two++;
1248  }
1249  else if (one.hasTerms() && one.exp() == j - k + 1)
1250  {
1251  tmp[l] += mulNTL ((bufFactors[l+1][k]+one.coeff()), Pi[l-1][k], b) -
1252  M (k + 1, l + 1);
1253  one++;
1254  }
1255  else if (two.hasTerms() && two.exp() == j - k + 1)
1256  {
1257  tmp[l] += mulNTL (bufFactors[l+1][k], (Pi[l-1][k] + two.coeff()), b)
1258  - M (k + 1, l + 1);
1259  two++;
1260  }
1261  }
1262  else
1263  tmp[l] += M (k + 1, l + 1);
1264  }
1265  }
1266  if (b.getp() != 0)
1267  tmp[l]= b (tmp[l]);
1268  Pi[l] += tmp[l]*xToJ*F.mvar();
1269  }
1270 }
1271 
1272 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // diopantineQa
1273 void
1274 henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1275  CFList& diophant, CFMatrix& M, modpk& b, bool sort)
1276 {
1277  if (sort)
1278  sortList (factors, Variable (1));
1279  Pi= CFArray (factors.length() - 1);
1280  CFListIterator j= factors;
1281  diophant= diophantine (F[0], F, factors, b);
1282  CanonicalForm bufF= F;
1283  if (getCharacteristic() == 0 && b.getp() != 0)
1284  {
1285  Variable v;
1286  bool hasAlgVar= hasFirstAlgVar (F, v);
1287  for (CFListIterator i= factors; i.hasItem() && !hasAlgVar; i++)
1288  hasAlgVar= hasFirstAlgVar (i.getItem(), v);
1289  Variable w;
1290  bool hasAlgVar2= false;
1291  for (CFListIterator i= diophant; i.hasItem() && !hasAlgVar2; i++)
1292  hasAlgVar2= hasFirstAlgVar (i.getItem(), w);
1293  if (hasAlgVar && hasAlgVar2 && v!=w)
1294  {
1295  bufF= replacevar (bufF, v, w);
1296  for (CFListIterator i= factors; i.hasItem(); i++)
1297  i.getItem()= replacevar (i.getItem(), v, w);
1298  }
1299  }
1300 
1301  DEBOUTLN (cerr, "diophant= " << diophant);
1302  j++;
1303  Pi [0]= mulNTL (j.getItem(), mod (factors.getFirst(), F.mvar()), b);
1304  M (1, 1)= Pi [0];
1305  int i= 1;
1306  if (j.hasItem())
1307  j++;
1308  for (; j.hasItem(); j++, i++)
1309  {
1310  Pi [i]= mulNTL (Pi [i - 1], j.getItem(), b);
1311  M (1, i + 1)= Pi [i];
1312  }
1313  CFArray bufFactors= CFArray (factors.length());
1314  i= 0;
1315  for (CFListIterator k= factors; k.hasItem(); i++, k++)
1316  {
1317  if (i == 0)
1318  bufFactors[i]= mod (k.getItem(), F.mvar());
1319  else
1320  bufFactors[i]= k.getItem();
1321  }
1322  for (i= 1; i < l; i++)
1323  henselStep12 (bufF, factors, bufFactors, diophant, M, Pi, i, b);
1324 
1325  CFListIterator k= factors;
1326  for (i= 0; i < factors.length (); i++, k++)
1327  k.getItem()= bufFactors[i];
1328  factors.removeFirst();
1329 }
1330 #endif
1331 
1332 #if defined(HAVE_NTL) || defined(HAVE_FLINT) //henselLift12
1333 void
1334 henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
1335  CFList& diophant, CFMatrix& M, bool sort)
1336 {
1337  modpk dummy= modpk();
1338  henselLift12 (F, factors, l, Pi, diophant, M, dummy, sort);
1339 }
1340 #endif
1341 
1342 void
1343 henselLiftResume12 (const CanonicalForm& F, CFList& factors, int start, int
1344  end, CFArray& Pi, const CFList& diophant, CFMatrix& M,
1345  const modpk& b)
1346 {
1347  CFArray bufFactors= CFArray (factors.length());
1348  int i= 0;
1349  CanonicalForm xToStart= power (F.mvar(), start);
1350  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1351  {
1352  if (i == 0)
1353  bufFactors[i]= mod (k.getItem(), xToStart);
1354  else
1355  bufFactors[i]= k.getItem();
1356  }
1357  for (i= start; i < end; i++)
1358  henselStep12 (F, factors, bufFactors, diophant, M, Pi, i, b);
1359 
1360  CFListIterator k= factors;
1361  for (i= 0; i < factors.length(); k++, i++)
1362  k.getItem()= bufFactors [i];
1363  factors.removeFirst();
1364  return;
1365 }
1366 
1367 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // diophantine
1368 CFList
1369 biDiophantine (const CanonicalForm& F, const CFList& factors, int d)
1370 {
1371  Variable y= F.mvar();
1372  CFList result;
1373  if (y.level() == 1)
1374  {
1375  result= diophantine (F, factors);
1376  return result;
1377  }
1378  else
1379  {
1380  CFList buf= factors;
1381  for (CFListIterator i= buf; i.hasItem(); i++)
1382  i.getItem()= mod (i.getItem(), y);
1383  CanonicalForm A= mod (F, y);
1384  int bufD= 1;
1385  CFList recResult= biDiophantine (A, buf, bufD);
1386  CanonicalForm e= 1;
1387  CFList p;
1388  CFArray bufFactors= CFArray (factors.length());
1389  CanonicalForm yToD= power (y, d);
1390  int k= 0;
1391  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1392  {
1393  bufFactors [k]= i.getItem();
1394  }
1395  CanonicalForm b, quot;
1396  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1397  {
1398  b= 1;
1399  if (fdivides (bufFactors[k], F, quot))
1400  b= quot;
1401  else
1402  {
1403  for (int l= 0; l < factors.length(); l++)
1404  {
1405  if (l == k)
1406  continue;
1407  else
1408  {
1409  b= mulMod2 (b, bufFactors[l], yToD);
1410  }
1411  }
1412  }
1413  p.append (b);
1414  }
1415 
1416  CFListIterator j= p;
1417  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1418  e -= i.getItem()*j.getItem();
1419 
1420  if (e.isZero())
1421  return recResult;
1422  CanonicalForm coeffE;
1423  CFList s;
1424  result= recResult;
1425  CanonicalForm g;
1426  for (int i= 1; i < d; i++)
1427  {
1428  if (degree (e, y) > 0)
1429  coeffE= e[i];
1430  else
1431  coeffE= 0;
1432  if (!coeffE.isZero())
1433  {
1435  CFListIterator l= p;
1436  int ii= 0;
1437  j= recResult;
1438  for (; j.hasItem(); j++, k++, l++, ii++)
1439  {
1440  g= coeffE*j.getItem();
1441  if (degree (bufFactors[ii], y) <= 0)
1442  g= mod (g, bufFactors[ii]);
1443  else
1444  g= mod (g, bufFactors[ii][0]);
1445  k.getItem() += g*power (y, i);
1446  e -= mulMod2 (g*power(y, i), l.getItem(), yToD);
1447  DEBOUTLN (cerr, "mod (e, power (y, i + 1))= " <<
1448  mod (e, power (y, i + 1)));
1449  }
1450  }
1451  if (e.isZero())
1452  break;
1453  }
1454 
1455  DEBOUTLN (cerr, "mod (e, y)= " << mod (e, y));
1456 
1457 #ifdef DEBUGOUTPUT
1458  CanonicalForm test= 0;
1459  j= p;
1460  for (CFListIterator i= result; i.hasItem(); i++, j++)
1461  test += mod (i.getItem()*j.getItem(), power (y, d));
1462  DEBOUTLN (cerr, "test= " << test);
1463 #endif
1464  return result;
1465  }
1466 }
1467 #endif
1468 
1469 CFList
1470 multiRecDiophantine (const CanonicalForm& F, const CFList& factors,
1471  const CFList& recResult, const CFList& M, int d)
1472 {
1473  Variable y= F.mvar();
1474  CFList result;
1475  CFListIterator i;
1476  CanonicalForm e= 1;
1477  CFListIterator j= factors;
1478  CFList p;
1479  CFArray bufFactors= CFArray (factors.length());
1480  CanonicalForm yToD= power (y, d);
1481  int k= 0;
1482  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1483  bufFactors [k]= i.getItem();
1484  CanonicalForm b, quot;
1485  CFList buf= M;
1486  buf.removeLast();
1487  buf.append (yToD);
1488  for (k= 0; k < factors.length(); k++) //TODO compute b's faster
1489  {
1490  b= 1;
1491  if (fdivides (bufFactors[k], F, quot))
1492  b= quot;
1493  else
1494  {
1495  for (int l= 0; l < factors.length(); l++)
1496  {
1497  if (l == k)
1498  continue;
1499  else
1500  {
1501  b= mulMod (b, bufFactors[l], buf);
1502  }
1503  }
1504  }
1505  p.append (b);
1506  }
1507  j= p;
1508  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
1509  e -= mulMod (i.getItem(), j.getItem(), M);
1510 
1511  if (e.isZero())
1512  return recResult;
1513  CanonicalForm coeffE;
1514  CFList s;
1515  result= recResult;
1516  CanonicalForm g;
1517  for (int i= 1; i < d; i++)
1518  {
1519  if (degree (e, y) > 0)
1520  coeffE= e[i];
1521  else
1522  coeffE= 0;
1523  if (!coeffE.isZero())
1524  {
1526  CFListIterator l= p;
1527  j= recResult;
1528  int ii= 0;
1529  CanonicalForm dummy;
1530  for (; j.hasItem(); j++, k++, l++, ii++)
1531  {
1532  g= mulMod (coeffE, j.getItem(), M);
1533  if (degree (bufFactors[ii], y) <= 0)
1534  divrem (g, mod (bufFactors[ii], Variable (y.level() - 1)), dummy,
1535  g, M);
1536  else
1537  divrem (g, bufFactors[ii][0], dummy, g, M);
1538  k.getItem() += g*power (y, i);
1539  e -= mulMod (g*power (y, i), l.getItem(), M);
1540  }
1541  }
1542 
1543  if (e.isZero())
1544  break;
1545  }
1546 
1547 #ifdef DEBUGOUTPUT
1548  CanonicalForm test= 0;
1549  j= p;
1550  for (CFListIterator i= result; i.hasItem(); i++, j++)
1551  test += mod (i.getItem()*j.getItem(), power (y, d));
1552  DEBOUTLN (cerr, "test in multiRecDiophantine= " << test);
1553 #endif
1554  return result;
1555 }
1556 
1557 static
1558 void
1559 henselStep (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
1560  const CFList& diophant, CFMatrix& M, CFArray& Pi, int j,
1561  const CFList& MOD)
1562 {
1563  CanonicalForm E;
1564  CanonicalForm xToJ= power (F.mvar(), j);
1565  Variable x= F.mvar();
1566  // compute the error
1567  if (j == 1)
1568  {
1569  E= F[j];
1570 #ifdef DEBUGOUTPUT
1571  CanonicalForm test= 1;
1572  for (int i= 0; i < factors.length(); i++)
1573  {
1574  if (i == 0)
1575  test= mulMod (test, mod (bufFactors [i], xToJ), MOD);
1576  else
1577  test= mulMod (test, bufFactors[i], MOD);
1578  }
1579  CanonicalForm test2= mod (F-test, xToJ);
1580 
1581  test2= mod (test2, MOD);
1582  DEBOUTLN (cerr, "test in henselStep= " << test2);
1583 #endif
1584  }
1585  else
1586  {
1587 #ifdef DEBUGOUTPUT
1588  CanonicalForm test= 1;
1589  for (int i= 0; i < factors.length(); i++)
1590  {
1591  if (i == 0)
1592  test *= mod (bufFactors [i], power (x, j));
1593  else
1594  test *= bufFactors[i];
1595  }
1596  test= mod (test, power (x, j));
1597  test= mod (test, MOD);
1598  CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
1599  DEBOUTLN (cerr, "test in henselStep= " << test2);
1600 #endif
1601 
1602  if (degree (Pi [factors.length() - 2], x) > 0)
1603  E= F[j] - Pi [factors.length() - 2] [j];
1604  else
1605  E= F[j];
1606  }
1607 
1608  CFArray buf= CFArray (diophant.length());
1609  bufFactors[0]= mod (factors.getFirst(), power (F.mvar(), j + 1));
1610  int k= 0;
1611  // actual lifting
1612  CanonicalForm dummy, rest1;
1613  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1614  {
1615  if (degree (bufFactors[k], x) > 0)
1616  {
1617  if (k > 0)
1618  divrem (E, bufFactors[k] [0], dummy, rest1, MOD);
1619  else
1620  rest1= E;
1621  }
1622  else
1623  divrem (E, bufFactors[k], dummy, rest1, MOD);
1624 
1625  buf[k]= mulMod (i.getItem(), rest1, MOD);
1626 
1627  if (degree (bufFactors[k], x) > 0)
1628  divrem (buf[k], bufFactors[k] [0], dummy, buf[k], MOD);
1629  else
1630  divrem (buf[k], bufFactors[k], dummy, buf[k], MOD);
1631  }
1632  for (k= 1; k < factors.length(); k++)
1633  bufFactors[k] += xToJ*buf[k];
1634 
1635  // update Pi [0]
1636  int degBuf0= degree (bufFactors[0], x);
1637  int degBuf1= degree (bufFactors[1], x);
1638  if (degBuf0 > 0 && degBuf1 > 0)
1639  M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
1640  CanonicalForm uIZeroJ;
1641 
1642  if (degBuf0 > 0 && degBuf1 > 0)
1643  uIZeroJ= mulMod ((bufFactors[0] [0] + bufFactors[0] [j]),
1644  (bufFactors[1] [0] + buf[1]), MOD) - M(1, 1) - M(j + 1, 1);
1645  else if (degBuf0 > 0)
1646  uIZeroJ= mulMod (bufFactors[0] [j], bufFactors[1], MOD);
1647  else if (degBuf1 > 0)
1648  uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
1649  else
1650  uIZeroJ= 0;
1651  Pi [0] += xToJ*uIZeroJ;
1652 
1653  CFArray tmp= CFArray (factors.length() - 1);
1654  for (k= 0; k < factors.length() - 1; k++)
1655  tmp[k]= 0;
1656  CFIterator one, two;
1657  one= bufFactors [0];
1658  two= bufFactors [1];
1659  if (degBuf0 > 0 && degBuf1 > 0)
1660  {
1661  for (k= 1; k <= (j+1)/2; k++)
1662  {
1663  if (k != j - k + 1)
1664  {
1665  if ((one.hasTerms() && one.exp() == j - k + 1) &&
1666  (two.hasTerms() && two.exp() == j - k + 1))
1667  {
1668  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
1669  (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
1670  M (j - k + 2, 1);
1671  one++;
1672  two++;
1673  }
1674  else if (one.hasTerms() && one.exp() == j - k + 1)
1675  {
1676  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
1677  bufFactors[1] [k], MOD) - M (k + 1, 1);
1678  one++;
1679  }
1680  else if (two.hasTerms() && two.exp() == j - k + 1)
1681  {
1682  tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
1683  two.coeff()), MOD) - M (k + 1, 1);
1684  two++;
1685  }
1686  }
1687  else
1688  {
1689  tmp[0] += M (k + 1, 1);
1690  }
1691  }
1692  }
1693  Pi [0] += tmp[0]*xToJ*F.mvar();
1694 
1695  // update Pi [l]
1696  int degPi, degBuf;
1697  for (int l= 1; l < factors.length() - 1; l++)
1698  {
1699  degPi= degree (Pi [l - 1], x);
1700  degBuf= degree (bufFactors[l + 1], x);
1701  if (degPi > 0 && degBuf > 0)
1702  M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
1703  if (j == 1)
1704  {
1705  if (degPi > 0 && degBuf > 0)
1706  Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
1707  (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
1708  M (1, l + 1));
1709  else if (degPi > 0)
1710  Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
1711  else if (degBuf > 0)
1712  Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
1713  }
1714  else
1715  {
1716  if (degPi > 0 && degBuf > 0)
1717  {
1718  uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
1719  uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
1720  }
1721  else if (degPi > 0)
1722  uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
1723  else if (degBuf > 0)
1724  {
1725  uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
1726  uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
1727  }
1728  Pi[l] += xToJ*uIZeroJ;
1729  }
1730  one= bufFactors [l + 1];
1731  two= Pi [l - 1];
1732  if (two.hasTerms() && two.exp() == j + 1)
1733  {
1734  if (degBuf > 0 && degPi > 0)
1735  {
1736  tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
1737  two++;
1738  }
1739  else if (degPi > 0)
1740  {
1741  tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
1742  two++;
1743  }
1744  }
1745  if (degBuf > 0 && degPi > 0)
1746  {
1747  for (k= 1; k <= (j+1)/2; k++)
1748  {
1749  if (k != j - k + 1)
1750  {
1751  if ((one.hasTerms() && one.exp() == j - k + 1) &&
1752  (two.hasTerms() && two.exp() == j - k + 1))
1753  {
1754  tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1755  (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
1756  M (j - k + 2, l + 1);
1757  one++;
1758  two++;
1759  }
1760  else if (one.hasTerms() && one.exp() == j - k + 1)
1761  {
1762  tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
1763  Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
1764  one++;
1765  }
1766  else if (two.hasTerms() && two.exp() == j - k + 1)
1767  {
1768  tmp[l] += mulMod (bufFactors[l + 1] [k],
1769  (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
1770  two++;
1771  }
1772  }
1773  else
1774  tmp[l] += M (k + 1, l + 1);
1775  }
1776  }
1777  Pi[l] += tmp[l]*xToJ*F.mvar();
1778  }
1779 
1780  return;
1781 }
1782 
1783 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // biDiophantine
1784 CFList
1785 henselLift23 (const CFList& eval, const CFList& factors, int* l, CFList&
1786  diophant, CFArray& Pi, CFMatrix& M)
1787 {
1788  CFList buf= factors;
1789  int k= 0;
1790  int liftBoundBivar= l[k];
1791  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
1792  CFList MOD;
1793  MOD.append (power (Variable (2), liftBoundBivar));
1794  CFArray bufFactors= CFArray (factors.length());
1795  k= 0;
1797  j++;
1798  buf.removeFirst();
1799  buf.insert (LC (j.getItem(), 1));
1800  for (CFListIterator i= buf; i.hasItem(); i++, k++)
1801  bufFactors[k]= i.getItem();
1802  Pi= CFArray (factors.length() - 1);
1803  CFListIterator i= buf;
1804  i++;
1805  Variable y= j.getItem().mvar();
1806  Pi [0]= mulMod (i.getItem(), mod (buf.getFirst(), y), MOD);
1807  M (1, 1)= Pi [0];
1808  k= 1;
1809  if (i.hasItem())
1810  i++;
1811  for (; i.hasItem(); i++, k++)
1812  {
1813  Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
1814  M (1, k + 1)= Pi [k];
1815  }
1816 
1817  for (int d= 1; d < l[1]; d++)
1818  henselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, d, MOD);
1819  CFList result;
1820  for (k= 1; k < factors.length(); k++)
1821  result.append (bufFactors[k]);
1822  return result;
1823 }
1824 #endif
1825 
1826 void
1827 henselLiftResume (const CanonicalForm& F, CFList& factors, int start, int end,
1828  CFArray& Pi, const CFList& diophant, CFMatrix& M,
1829  const CFList& MOD)
1830 {
1831  CFArray bufFactors= CFArray (factors.length());
1832  int i= 0;
1833  CanonicalForm xToStart= power (F.mvar(), start);
1834  for (CFListIterator k= factors; k.hasItem(); k++, i++)
1835  {
1836  if (i == 0)
1837  bufFactors[i]= mod (k.getItem(), xToStart);
1838  else
1839  bufFactors[i]= k.getItem();
1840  }
1841  for (i= start; i < end; i++)
1842  henselStep (F, factors, bufFactors, diophant, M, Pi, i, MOD);
1843 
1844  CFListIterator k= factors;
1845  for (i= 0; i < factors.length(); k++, i++)
1846  k.getItem()= bufFactors [i];
1847  factors.removeFirst();
1848  return;
1849 }
1850 
1851 CFList
1852 henselLift (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
1853  diophant, CFArray& Pi, CFMatrix& M, int lOld, int lNew)
1854 {
1855  diophant= multiRecDiophantine (F.getFirst(), factors, diophant, MOD, lOld);
1856  int k= 0;
1857  CFArray bufFactors= CFArray (factors.length());
1858  for (CFListIterator i= factors; i.hasItem(); i++, k++)
1859  {
1860  if (k == 0)
1861  bufFactors[k]= LC (F.getLast(), 1);
1862  else
1863  bufFactors[k]= i.getItem();
1864  }
1865  CFList buf= factors;
1866  buf.removeFirst();
1867  buf.insert (LC (F.getLast(), 1));
1868  CFListIterator i= buf;
1869  i++;
1870  Variable y= F.getLast().mvar();
1871  Variable x= F.getFirst().mvar();
1872  CanonicalForm xToLOld= power (x, lOld);
1873  Pi [0]= mod (Pi[0], xToLOld);
1874  M (1, 1)= Pi [0];
1875  k= 1;
1876  if (i.hasItem())
1877  i++;
1878  for (; i.hasItem(); i++, k++)
1879  {
1880  Pi [k]= mod (Pi [k], xToLOld);
1881  M (1, k + 1)= Pi [k];
1882  }
1883 
1884  for (int d= 1; d < lNew; d++)
1885  henselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, d, MOD);
1886  CFList result;
1887  for (k= 1; k < factors.length(); k++)
1888  result.append (bufFactors[k]);
1889  return result;
1890 }
1891 
1892 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // henselLift23
1893 CFList
1894 henselLift (const CFList& eval, const CFList& factors, int* l, int lLength,
1895  bool sort)
1896 {
1897  CFList diophant;
1898  CFList buf= factors;
1899  buf.insert (LC (eval.getFirst(), 1));
1900  if (sort)
1901  sortList (buf, Variable (1));
1902  CFArray Pi;
1903  CFMatrix M= CFMatrix (l[1], factors.length());
1904  CFList result= henselLift23 (eval, buf, l, diophant, Pi, M);
1905  if (eval.length() == 2)
1906  return result;
1907  CFList MOD;
1908  for (int i= 0; i < 2; i++)
1909  MOD.append (power (Variable (i + 2), l[i]));
1911  j++;
1912  CFList bufEval;
1913  bufEval.append (j.getItem());
1914  j++;
1915 
1916  for (int i= 2; i < lLength && j.hasItem(); i++, j++)
1917  {
1918  result.insert (LC (bufEval.getFirst(), 1));
1919  bufEval.append (j.getItem());
1920  M= CFMatrix (l[i], factors.length());
1921  result= henselLift (bufEval, result, MOD, diophant, Pi, M, l[i - 1], l[i]);
1922  MOD.append (power (Variable (i + 2), l[i]));
1923  bufEval.removeFirst();
1924  }
1925  return result;
1926 }
1927 #endif
1928 
1929 // nonmonic
1930 
1931 void
1932 nonMonicHenselStep12 (const CanonicalForm& F, const CFList& factors,
1933  CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
1934  CFArray& Pi, int j, const CFArray& /*LCs*/)
1935 {
1936  Variable x= F.mvar();
1937  CanonicalForm xToJ= power (x, j);
1938 
1939  CanonicalForm E;
1940  // compute the error
1941  if (degree (Pi [factors.length() - 2], x) > 0)
1942  E= F[j] - Pi [factors.length() - 2] [j];
1943  else
1944  E= F[j];
1945 
1946  CFArray buf= CFArray (diophant.length());
1947 
1948  int k= 0;
1950  // actual lifting
1951  for (CFListIterator i= diophant; i.hasItem(); i++, k++)
1952  {
1953  if (degree (bufFactors[k], x) > 0)
1954  remainder= modNTL (E, bufFactors[k] [0]);
1955  else
1956  remainder= modNTL (E, bufFactors[k]);
1957  buf[k]= mulNTL (i.getItem(), remainder);
1958  if (degree (bufFactors[k], x) > 0)
1959  buf[k]= modNTL (buf[k], bufFactors[k] [0]);
1960  else
1961  buf[k]= modNTL (buf[k], bufFactors[k]);
1962  }
1963 
1964  for (k= 0; k < factors.length(); k++)
1965  bufFactors[k] += xToJ*buf[k];
1966 
1967  // update Pi [0]
1968  int degBuf0= degree (bufFactors[0], x);
1969  int degBuf1= degree (bufFactors[1], x);
1970  if (degBuf0 > 0 && degBuf1 > 0)
1971  {
1972  M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
1973  if (j + 2 <= M.rows())
1974  M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
1975  }
1976  else
1977  M (j + 1, 1)= 0;
1978 
1979  CanonicalForm uIZeroJ;
1980  if (degBuf0 > 0 && degBuf1 > 0)
1981  uIZeroJ= mulNTL(bufFactors[0][0], buf[1]) +
1982  mulNTL (bufFactors[1][0], buf[0]);
1983  else if (degBuf0 > 0)
1984  uIZeroJ= mulNTL (buf[0], bufFactors[1]) +
1985  mulNTL (buf[1], bufFactors[0][0]);
1986  else if (degBuf1 > 0)
1987  uIZeroJ= mulNTL (bufFactors[0], buf[1]) +
1988  mulNTL (buf[0], bufFactors[1][0]);
1989  else
1990  uIZeroJ= mulNTL (bufFactors[0], buf[1]) +
1991  mulNTL (buf[0], bufFactors[1]);
1992 
1993  Pi [0] += xToJ*uIZeroJ;
1994 
1995  CFArray tmp= CFArray (factors.length() - 1);
1996  for (k= 0; k < factors.length() - 1; k++)
1997  tmp[k]= 0;
1998  CFIterator one, two;
1999  one= bufFactors [0];
2000  two= bufFactors [1];
2001  if (degBuf0 > 0 && degBuf1 > 0)
2002  {
2003  while (one.hasTerms() && one.exp() > j) one++;
2004  while (two.hasTerms() && two.exp() > j) two++;
2005  for (k= 1; k <= (j+1)/2; k++)
2006  {
2007  if (k != j - k + 1)
2008  {
2009  if ((one.hasTerms() && one.exp() == j - k + 1) && +
2010  (two.hasTerms() && two.exp() == j - k + 1))
2011  {
2012  tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()),(bufFactors[1][k] +
2013  two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
2014  one++;
2015  two++;
2016  }
2017  else if (one.hasTerms() && one.exp() == j - k + 1)
2018  {
2019  tmp[0] += mulNTL ((bufFactors[0][k]+one.coeff()), bufFactors[1] [k]) -
2020  M (k + 1, 1);
2021  one++;
2022  }
2023  else if (two.hasTerms() && two.exp() == j - k + 1)
2024  {
2025  tmp[0] += mulNTL (bufFactors[0][k],(bufFactors[1][k] + two.coeff())) -
2026  M (k + 1, 1);
2027  two++;
2028  }
2029  }
2030  else
2031  tmp[0] += M (k + 1, 1);
2032  }
2033  }
2034 
2035  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
2036  {
2037  if (j + 2 <= M.rows())
2038  tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
2039  (bufFactors [1] [j + 1] + bufFactors [1] [0]))
2040  - M(1,1) - M (j + 2,1);
2041  }
2042  else if (degBuf0 >= j + 1)
2043  {
2044  if (degBuf1 > 0)
2045  tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
2046  else
2047  tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1]);
2048  }
2049  else if (degBuf1 >= j + 1)
2050  {
2051  if (degBuf0 > 0)
2052  tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
2053  else
2054  tmp[0] += mulNTL (bufFactors [0], bufFactors [1] [j + 1]);
2055  }
2056 
2057  Pi [0] += tmp[0]*xToJ*F.mvar();
2058 
2059  int degPi, degBuf;
2060  for (int l= 1; l < factors.length() - 1; l++)
2061  {
2062  degPi= degree (Pi [l - 1], x);
2063  degBuf= degree (bufFactors[l + 1], x);
2064  if (degPi > 0 && degBuf > 0)
2065  {
2066  M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]);
2067  if (j + 2 <= M.rows())
2068  M (j + 2, l + 1)= mulNTL (Pi [l - 1][j + 1], bufFactors[l + 1] [j + 1]);
2069  }
2070  else
2071  M (j + 1, l + 1)= 0;
2072 
2073  if (degPi > 0 && degBuf > 0)
2074  uIZeroJ= mulNTL (Pi[l - 1] [0], buf[l + 1]) +
2075  mulNTL (uIZeroJ, bufFactors[l+1] [0]);
2076  else if (degPi > 0)
2077  uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]) +
2078  mulNTL (Pi[l - 1][0], buf[l + 1]);
2079  else if (degBuf > 0)
2080  uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1][0]) +
2081  mulNTL (Pi[l - 1], buf[l + 1]);
2082  else
2083  uIZeroJ= mulNTL (uIZeroJ, bufFactors[l + 1]) +
2084  mulNTL (Pi[l - 1], buf[l + 1]);
2085 
2086  Pi [l] += xToJ*uIZeroJ;
2087 
2088  one= bufFactors [l + 1];
2089  two= Pi [l - 1];
2090  if (degBuf > 0 && degPi > 0)
2091  {
2092  while (one.hasTerms() && one.exp() > j) one++;
2093  while (two.hasTerms() && two.exp() > j) two++;
2094  for (k= 1; k <= (j+1)/2; k++)
2095  {
2096  if (k != j - k + 1)
2097  {
2098  if ((one.hasTerms() && one.exp() == j - k + 1) &&
2099  (two.hasTerms() && two.exp() == j - k + 1))
2100  {
2101  tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
2102  (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1) -
2103  M (j - k + 2, l + 1);
2104  one++;
2105  two++;
2106  }
2107  else if (one.hasTerms() && one.exp() == j - k + 1)
2108  {
2109  tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()),
2110  Pi[l - 1] [k]) - M (k + 1, l + 1);
2111  one++;
2112  }
2113  else if (two.hasTerms() && two.exp() == j - k + 1)
2114  {
2115  tmp[l] += mulNTL (bufFactors[l + 1] [k],
2116  (Pi[l - 1] [k] + two.coeff())) - M (k + 1, l + 1);
2117  two++;
2118  }
2119  }
2120  else
2121  tmp[l] += M (k + 1, l + 1);
2122  }
2123  }
2124 
2125  if (degPi >= j + 1 && degBuf >= j + 1)
2126  {
2127  if (j + 2 <= M.rows())
2128  tmp [l] += mulNTL ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
2129  (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
2130  ) - M(1,l+1) - M (j + 2,l+1);
2131  }
2132  else if (degPi >= j + 1)
2133  {
2134  if (degBuf > 0)
2135  tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1] [0]);
2136  else
2137  tmp[l] += mulNTL (Pi [l - 1] [j+1], bufFactors [l + 1]);
2138  }
2139  else if (degBuf >= j + 1)
2140  {
2141  if (degPi > 0)
2142  tmp[l] += mulNTL (Pi [l - 1] [0], bufFactors [l + 1] [j + 1]);
2143  else
2144  tmp[l] += mulNTL (Pi [l - 1], bufFactors [l + 1] [j + 1]);
2145  }
2146 
2147  Pi[l] += tmp[l]*xToJ*F.mvar();
2148  }
2149  return;
2150 }
2151 
2152 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // diophantine
2153 void
2154 nonMonicHenselLift12 (const CanonicalForm& F, CFList& factors, int l,
2155  CFArray& Pi, CFList& diophant, CFMatrix& M,
2156  const CFArray& LCs, bool sort)
2157 {
2158  if (sort)
2159  sortList (factors, Variable (1));
2160  Pi= CFArray (factors.length() - 2);
2161  CFList bufFactors2= factors;
2162  bufFactors2.removeFirst();
2163  diophant= diophantine (F[0], bufFactors2);
2164  DEBOUTLN (cerr, "diophant= " << diophant);
2165 
2166  CFArray bufFactors= CFArray (bufFactors2.length());
2167  int i= 0;
2168  for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
2169  bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
2170 
2171  Variable x= F.mvar();
2172  if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
2173  {
2174  M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
2175  Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
2176  mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
2177  }
2178  else if (degree (bufFactors[0], x) > 0)
2179  {
2180  M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
2181  Pi [0]= M (1, 1) +
2182  mulNTL (bufFactors [0] [1], bufFactors[1])*x;
2183  }
2184  else if (degree (bufFactors[1], x) > 0)
2185  {
2186  M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
2187  Pi [0]= M (1, 1) +
2188  mulNTL (bufFactors [0], bufFactors[1] [1])*x;
2189  }
2190  else
2191  {
2192  M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
2193  Pi [0]= M (1, 1);
2194  }
2195 
2196  for (i= 1; i < Pi.size(); i++)
2197  {
2198  if (degree (Pi[i-1], x) > 0 && degree (bufFactors [i+1], x) > 0)
2199  {
2200  M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors[i+1] [0]);
2201  Pi [i]= M (1,i+1) + (mulNTL (Pi[i-1] [1], bufFactors[i+1] [0]) +
2202  mulNTL (Pi[i-1] [0], bufFactors [i+1] [1]))*x;
2203  }
2204  else if (degree (Pi[i-1], x) > 0)
2205  {
2206  M (1,i+1)= mulNTL (Pi[i-1] [0], bufFactors [i+1]);
2207  Pi [i]= M(1,i+1) + mulNTL (Pi[i-1] [1], bufFactors[i+1])*x;
2208  }
2209  else if (degree (bufFactors[i+1], x) > 0)
2210  {
2211  M (1,i+1)= mulNTL (Pi[i-1], bufFactors [i+1] [0]);
2212  Pi [i]= M (1,i+1) + mulNTL (Pi[i-1], bufFactors[i+1] [1])*x;
2213  }
2214  else
2215  {
2216  M (1,i+1)= mulNTL (Pi [i-1], bufFactors [i+1]);
2217  Pi [i]= M (1,i+1);
2218  }
2219  }
2220 
2221  for (i= 1; i < l; i++)
2222  nonMonicHenselStep12 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
2223 
2224  factors= CFList();
2225  for (i= 0; i < bufFactors.size(); i++)
2226  factors.append (bufFactors[i]);
2227  return;
2228 }
2229 #endif
2230 
2231 /// solve \f$ E=\sum_{i= 1}^r{\sigma_{i}\prod_{j=1, j\neq i}^{r}{f_{j}}} \f$
2232 /// mod M, @a products contains \f$ \prod_{j=1, j\neq i}^{r}{f_{j}} \f$
2233 CFList
2234 diophantine (const CFList& recResult, const CFList& factors,
2235  const CFList& products, const CFList& M, const CanonicalForm& E,
2236  bool& bad)
2237 {
2238  if (M.isEmpty())
2239  {
2240  CFList result;
2241  CFListIterator j= factors;
2243  for (CFListIterator i= recResult; i.hasItem(); i++, j++)
2244  {
2245  ASSERT (E.isUnivariate() || E.inCoeffDomain(),
2246  "constant or univariate poly expected");
2247  ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
2248  "constant or univariate poly expected");
2249  ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
2250  "constant or univariate poly expected");
2251  buf= mulNTL (E, i.getItem());
2252  result.append (modNTL (buf, j.getItem()));
2253  }
2254  return result;
2255  }
2256  Variable y= M.getLast().mvar();
2257  CFList bufFactors= factors;
2258  for (CFListIterator i= bufFactors; i.hasItem(); i++)
2259  i.getItem()= mod (i.getItem(), y);
2260  CFList bufProducts= products;
2261  for (CFListIterator i= bufProducts; i.hasItem(); i++)
2262  i.getItem()= mod (i.getItem(), y);
2263  CFList buf= M;
2264  buf.removeLast();
2265  CanonicalForm bufE= mod (E, y);
2266  CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2267  bufE, bad);
2268 
2269  if (bad)
2270  return CFList();
2271 
2272  CanonicalForm e= E;
2273  CFListIterator j= products;
2274  for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
2275  e -= j.getItem()*i.getItem();
2276 
2277  CFList result= recDiophantine;
2278  int d= degree (M.getLast());
2279  CanonicalForm coeffE;
2280  for (int i= 1; i < d; i++)
2281  {
2282  if (degree (e, y) > 0)
2283  coeffE= e[i];
2284  else
2285  coeffE= 0;
2286  if (!coeffE.isZero())
2287  {
2289  recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
2290  coeffE, bad);
2291  if (bad)
2292  return CFList();
2293  CFListIterator l= products;
2294  for (j= recDiophantine; j.hasItem(); j++, k++, l++)
2295  {
2296  k.getItem() += j.getItem()*power (y, i);
2297  e -= l.getItem()*(j.getItem()*power (y, i));
2298  }
2299  }
2300  if (e.isZero())
2301  break;
2302  }
2303  if (!e.isZero())
2304  {
2305  bad= true;
2306  return CFList();
2307  }
2308  return result;
2309 }
2310 
2311 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // diophantine
2312 void
2313 nonMonicHenselStep (const CanonicalForm& F, const CFList& factors,
2314  CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
2315  CFArray& Pi, const CFList& products, int j,
2316  const CFList& MOD, bool& noOneToOne)
2317 {
2318  CanonicalForm E;
2319  CanonicalForm xToJ= power (F.mvar(), j);
2320  Variable x= F.mvar();
2321 
2322  // compute the error
2323 #ifdef DEBUGOUTPUT
2324  CanonicalForm test= 1;
2325  for (int i= 0; i < factors.length(); i++)
2326  {
2327  if (i == 0)
2328  test *= mod (bufFactors [i], power (x, j));
2329  else
2330  test *= bufFactors[i];
2331  }
2332  test= mod (test, power (x, j));
2333  test= mod (test, MOD);
2334  CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
2335  DEBOUTLN (cerr, "test in nonMonicHenselStep= " << test2);
2336 #endif
2337 
2338  if (degree (Pi [factors.length() - 2], x) > 0)
2339  E= F[j] - Pi [factors.length() - 2] [j];
2340  else
2341  E= F[j];
2342 
2343  CFArray buf= CFArray (diophant.length());
2344 
2345  // actual lifting
2346  TIMING_START (diotime);
2347  CFList diophantine2= diophantine (diophant, factors, products, MOD, E,
2348  noOneToOne);
2349 
2350  if (noOneToOne)
2351  return;
2352 
2353  int k= 0;
2354  for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
2355  {
2356  buf[k]= i.getItem();
2357  bufFactors[k] += xToJ*i.getItem();
2358  }
2359  TIMING_END_AND_PRINT (diotime, "time for dio: ");
2360 
2361  // update Pi [0]
2362  TIMING_START (product2);
2363  int degBuf0= degree (bufFactors[0], x);
2364  int degBuf1= degree (bufFactors[1], x);
2365  if (degBuf0 > 0 && degBuf1 > 0)
2366  {
2367  M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
2368  if (j + 2 <= M.rows())
2369  M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
2370  }
2371  else
2372  M (j + 1, 1)= 0;
2373 
2374  CanonicalForm uIZeroJ;
2375  if (degBuf0 > 0 && degBuf1 > 0)
2376  uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
2377  mulMod (bufFactors[1] [0], buf[0], MOD);
2378  else if (degBuf0 > 0)
2379  uIZeroJ= mulMod (buf[0], bufFactors[1], MOD) +
2380  mulMod (buf[1], bufFactors[0][0], MOD);
2381  else if (degBuf1 > 0)
2382  uIZeroJ= mulMod (bufFactors[0], buf[1], MOD) +
2383  mulMod (buf[0], bufFactors[1][0], MOD);
2384  else
2385  uIZeroJ= mulMod (bufFactors[0], buf[1], MOD) +
2386  mulMod (buf[0], bufFactors[1], MOD);
2387  Pi [0] += xToJ*uIZeroJ;
2388 
2389  CFArray tmp= CFArray (factors.length() - 1);
2390  for (k= 0; k < factors.length() - 1; k++)
2391  tmp[k]= 0;
2392  CFIterator one, two;
2393  one= bufFactors [0];
2394  two= bufFactors [1];
2395  if (degBuf0 > 0 && degBuf1 > 0)
2396  {
2397  while (one.hasTerms() && one.exp() > j) one++;
2398  while (two.hasTerms() && two.exp() > j) two++;
2399  for (k= 1; k <= (j+1)/2; k++)
2400  {
2401  if (k != j - k + 1)
2402  {
2403  if ((one.hasTerms() && one.exp() == j - k + 1) &&
2404  (two.hasTerms() && two.exp() == j - k + 1))
2405  {
2406  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2407  (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
2408  M (j - k + 2, 1);
2409  one++;
2410  two++;
2411  }
2412  else if (one.hasTerms() && one.exp() == j - k + 1)
2413  {
2414  tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
2415  bufFactors[1] [k], MOD) - M (k + 1, 1);
2416  one++;
2417  }
2418  else if (two.hasTerms() && two.exp() == j - k + 1)
2419  {
2420  tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
2421  two.coeff()), MOD) - M (k + 1, 1);
2422  two++;
2423  }
2424  }
2425  else
2426  {
2427  tmp[0] += M (k + 1, 1);
2428  }
2429  }
2430  }
2431 
2432  if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
2433  {
2434  if (j + 2 <= M.rows())
2435  tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
2436  (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
2437  - M(1,1) - M (j + 2,1);
2438  }
2439  else if (degBuf0 >= j + 1)
2440  {
2441  if (degBuf1 > 0)
2442  tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
2443  else
2444  tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1], MOD);
2445  }
2446  else if (degBuf1 >= j + 1)
2447  {
2448  if (degBuf0 > 0)
2449  tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
2450  else
2451  tmp[0] += mulMod (bufFactors [0], bufFactors [1] [j + 1], MOD);
2452  }
2453  Pi [0] += tmp[0]*xToJ*F.mvar();
2454 
2455  // update Pi [l]
2456  int degPi, degBuf;
2457  for (int l= 1; l < factors.length() - 1; l++)
2458  {
2459  degPi= degree (Pi [l - 1], x);
2460  degBuf= degree (bufFactors[l + 1], x);
2461  if (degPi > 0 && degBuf > 0)
2462  {
2463  M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD);
2464  if (j + 2 <= M.rows())
2465  M (j + 2, l + 1)= mulMod (Pi [l - 1] [j + 1], bufFactors[l + 1] [j + 1],
2466  MOD);
2467  }
2468  else
2469  M (j + 1, l + 1)= 0;
2470 
2471  if (degPi > 0 && degBuf > 0)
2472  uIZeroJ= mulMod (Pi[l - 1] [0], buf[l + 1], MOD) +
2473  mulMod (uIZeroJ, bufFactors[l + 1] [0], MOD);
2474  else if (degPi > 0)
2475  uIZeroJ= mulMod (uIZeroJ, bufFactors[l + 1], MOD) +
2476  mulMod (Pi[l - 1][0], buf[l + 1], MOD);
2477  else if (degBuf > 0)
2478  uIZeroJ= mulMod (Pi[l - 1], buf[l + 1], MOD) +
2479  mulMod (uIZeroJ, bufFactors[l + 1][0], MOD);
2480  else
2481  uIZeroJ= mulMod (Pi[l - 1], buf[l + 1], MOD) +
2482  mulMod (uIZeroJ, bufFactors[l + 1], MOD);
2483 
2484  Pi [l] += xToJ*uIZeroJ;
2485 
2486  one= bufFactors [l + 1];
2487  two= Pi [l - 1];
2488  if (degBuf > 0 && degPi > 0)
2489  {
2490  while (one.hasTerms() && one.exp() > j) one++;
2491  while (two.hasTerms() && two.exp() > j) two++;
2492  for (k= 1; k <= (j+1)/2; k++)
2493  {
2494  if (k != j - k + 1)
2495  {
2496  if ((one.hasTerms() && one.exp() == j - k + 1) &&
2497  (two.hasTerms() && two.exp() == j - k + 1))
2498  {
2499  tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2500  (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
2501  M (j - k + 2, l + 1);
2502  one++;
2503  two++;
2504  }
2505  else if (one.hasTerms() && one.exp() == j - k + 1)
2506  {
2507  tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
2508  Pi[l - 1] [k], MOD) - M (k + 1, l + 1);
2509  one++;
2510  }
2511  else if (two.hasTerms() && two.exp() == j - k + 1)
2512  {
2513  tmp[l] += mulMod (bufFactors[l + 1] [k],
2514  (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
2515  two++;
2516  }
2517  }
2518  else
2519  tmp[l] += M (k + 1, l + 1);
2520  }
2521  }
2522 
2523  if (degPi >= j + 1 && degBuf >= j + 1)
2524  {
2525  if (j + 2 <= M.rows())
2526  tmp [l] += mulMod ((Pi [l - 1] [j + 1]+ Pi [l - 1] [0]),
2527  (bufFactors [l + 1] [j + 1] + bufFactors [l + 1] [0])
2528  , MOD) - M(1,l+1) - M (j + 2,l+1);
2529  }
2530  else if (degPi >= j + 1)
2531  {
2532  if (degBuf > 0)
2533  tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1] [0], MOD);
2534  else
2535  tmp[l] += mulMod (Pi [l - 1] [j+1], bufFactors [l + 1], MOD);
2536  }
2537  else if (degBuf >= j + 1)
2538  {
2539  if (degPi > 0)
2540  tmp[l] += mulMod (Pi [l - 1] [0], bufFactors [l + 1] [j + 1], MOD);
2541  else
2542  tmp[l] += mulMod (Pi [l - 1], bufFactors [l + 1] [j + 1], MOD);
2543  }
2544 
2545  Pi[l] += tmp[l]*xToJ*F.mvar();
2546  }
2547  TIMING_END_AND_PRINT (product2, "time for product in hensel step: ");
2548  return;
2549 }
2550 #endif
2551 
2552 // wrt. Variable (1)
2554 {
2555  if (degree (F, 1) <= 0)
2556  return c;
2557  else
2558  {
2559  CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
2560  result += (swapvar (c, Variable (F.level() + 1), Variable (1))
2561  - LC (result))*power (result.mvar(), degree (result));
2562  return swapvar (result, Variable (F.level() + 1), Variable (1));
2563  }
2564 }
2565 
2566 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // nonMonicHenselStep
2567 CFList
2568 nonMonicHenselLift232(const CFList& eval, const CFList& factors, int* l, CFList&
2569  diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
2570  const CFList& LCs2, bool& bad)
2571 {
2572  CFList buf= factors;
2573  int k= 0;
2574  int liftBoundBivar= l[k];
2575  CFList bufbuf= factors;
2576  Variable v= Variable (2);
2577 
2578  CFList MOD;
2579  MOD.append (power (Variable (2), liftBoundBivar));
2580  CFArray bufFactors= CFArray (factors.length());
2581  k= 0;
2583  j++;
2584  CFListIterator iter1= LCs1;
2585  CFListIterator iter2= LCs2;
2586  iter1++;
2587  iter2++;
2588  bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
2589  bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
2590 
2591  CFListIterator i= buf;
2592  i++;
2593  Variable y= j.getItem().mvar();
2594  if (y.level() != 3)
2595  y= Variable (3);
2596 
2597  Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
2598  M (1, 1)= Pi[0];
2599  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2600  Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2601  mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2602  else if (degree (bufFactors[0], y) > 0)
2603  Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2604  else if (degree (bufFactors[1], y) > 0)
2605  Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2606 
2607  CFList products;
2608  for (int i= 0; i < bufFactors.size(); i++)
2609  {
2610  if (degree (bufFactors[i], y) > 0)
2611  products.append (eval.getFirst()/bufFactors[i] [0]);
2612  else
2613  products.append (eval.getFirst()/bufFactors[i]);
2614  }
2615 
2616  for (int d= 1; d < l[1]; d++)
2617  {
2618  nonMonicHenselStep (j.getItem(), buf, bufFactors, diophant, M, Pi, products,
2619  d, MOD, bad);
2620  if (bad)
2621  return CFList();
2622  }
2623  CFList result;
2624  for (k= 0; k < factors.length(); k++)
2625  result.append (bufFactors[k]);
2626  return result;
2627 }
2628 #endif
2629 
2630 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // nonMonicHenselStep
2631 CFList
2632 nonMonicHenselLift2 (const CFList& F, const CFList& factors, const CFList& MOD,
2633  CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2634  int& lNew, const CFList& LCs1, const CFList& LCs2, bool& bad
2635  )
2636 {
2637  int k= 0;
2638  CFArray bufFactors= CFArray (factors.length());
2639  bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
2640  bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
2641  CFList buf= factors;
2642  Variable y= F.getLast().mvar();
2643  Variable x= F.getFirst().mvar();
2644  CanonicalForm xToLOld= power (x, lOld);
2645  Pi [0]= mod (Pi[0], xToLOld);
2646  M (1, 1)= Pi [0];
2647 
2648  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2649  Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2650  mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2651  else if (degree (bufFactors[0], y) > 0)
2652  Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2653  else if (degree (bufFactors[1], y) > 0)
2654  Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2655 
2656  CFList products;
2657  CanonicalForm quot;
2658  for (int i= 0; i < bufFactors.size(); i++)
2659  {
2660  if (degree (bufFactors[i], y) > 0)
2661  {
2662  if (!fdivides (bufFactors[i] [0], F.getFirst(), quot))
2663  {
2664  bad= true;
2665  return CFList();
2666  }
2667  products.append (quot);
2668  }
2669  else
2670  {
2671  if (!fdivides (bufFactors[i], F.getFirst(), quot))
2672  {
2673  bad= true;
2674  return CFList();
2675  }
2676  products.append (quot);
2677  }
2678  }
2679 
2680  for (int d= 1; d < lNew; d++)
2681  {
2682  nonMonicHenselStep (F.getLast(), buf, bufFactors, diophant, M, Pi, products,
2683  d, MOD, bad);
2684  if (bad)
2685  return CFList();
2686  }
2687 
2688  CFList result;
2689  for (k= 0; k < factors.length(); k++)
2690  result.append (bufFactors[k]);
2691  return result;
2692 }
2693 #endif
2694 
2695 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // nonMonicHenselStep
2696 CFList
2697 nonMonicHenselLift2 (const CFList& eval, const CFList& factors, int* l, int
2698  lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
2699  const CFArray& Pi, const CFList& diophant, bool& bad)
2700 {
2701  CFList bufDiophant= diophant;
2702  CFList buf= factors;
2703  if (sort)
2704  sortList (buf, Variable (1));
2705  CFArray bufPi= Pi;
2706  CFMatrix M= CFMatrix (l[1], factors.length());
2707  CFList result=
2708  nonMonicHenselLift232(eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2, bad);
2709  if (bad)
2710  return CFList();
2711 
2712  if (eval.length() == 2)
2713  return result;
2714  CFList MOD;
2715  for (int i= 0; i < 2; i++)
2716  MOD.append (power (Variable (i + 2), l[i]));
2718  j++;
2719  CFList bufEval;
2720  bufEval.append (j.getItem());
2721  j++;
2722  CFListIterator jj= LCs1;
2723  CFListIterator jjj= LCs2;
2724  CFList bufLCs1, bufLCs2;
2725  jj++, jjj++;
2726  bufLCs1.append (jj.getItem());
2727  bufLCs2.append (jjj.getItem());
2728  jj++, jjj++;
2729 
2730  for (int i= 2; i < lLength && j.hasItem(); i++, j++, jj++, jjj++)
2731  {
2732  bufEval.append (j.getItem());
2733  bufLCs1.append (jj.getItem());
2734  bufLCs2.append (jjj.getItem());
2735  M= CFMatrix (l[i], factors.length());
2736  result= nonMonicHenselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M,
2737  l[i - 1], l[i], bufLCs1, bufLCs2, bad);
2738  if (bad)
2739  return CFList();
2740  MOD.append (power (Variable (i + 2), l[i]));
2741  bufEval.removeFirst();
2742  bufLCs1.removeFirst();
2743  bufLCs2.removeFirst();
2744  }
2745  return result;
2746 }
2747 #endif
2748 
2749 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // diophantine
2750 CFList
2751 nonMonicHenselLift23 (const CanonicalForm& F, const CFList& factors, const
2752  CFList& LCs, CFList& diophant, CFArray& Pi, int liftBound,
2753  int bivarLiftBound, bool& bad)
2754 {
2755  CFList bufFactors2= factors;
2756 
2757  Variable y= Variable (2);
2758  for (CFListIterator i= bufFactors2; i.hasItem(); i++)
2759  i.getItem()= mod (i.getItem(), y);
2760 
2761  CanonicalForm bufF= F;
2762  bufF= mod (bufF, y);
2763  bufF= mod (bufF, Variable (3));
2764 
2765  TIMING_START (diotime);
2766  diophant= diophantine (bufF, bufFactors2);
2767  TIMING_END_AND_PRINT (diotime, "time for dio in 23: ");
2768 
2769  CFMatrix M= CFMatrix (liftBound, bufFactors2.length() - 1);
2770 
2771  Pi= CFArray (bufFactors2.length() - 1);
2772 
2773  CFArray bufFactors= CFArray (bufFactors2.length());
2774  CFListIterator j= LCs;
2775  int i= 0;
2776  for (CFListIterator k= factors; k.hasItem(); j++, k++, i++)
2777  bufFactors[i]= replaceLC (k.getItem(), j.getItem());
2778 
2779  //initialise Pi
2780  Variable v= Variable (3);
2781  CanonicalForm yToL= power (y, bivarLiftBound);
2782  if (degree (bufFactors[0], v) > 0 && degree (bufFactors [1], v) > 0)
2783  {
2784  M (1, 1)= mulMod2 (bufFactors [0] [0], bufFactors[1] [0], yToL);
2785  Pi [0]= M (1,1) + (mulMod2 (bufFactors [0] [1], bufFactors[1] [0], yToL) +
2786  mulMod2 (bufFactors [0] [0], bufFactors [1] [1], yToL))*v;
2787  }
2788  else if (degree (bufFactors[0], v) > 0)
2789  {
2790  M (1,1)= mulMod2 (bufFactors [0] [0], bufFactors [1], yToL);
2791  Pi [0]= M(1,1) + mulMod2 (bufFactors [0] [1], bufFactors[1], yToL)*v;
2792  }
2793  else if (degree (bufFactors[1], v) > 0)
2794  {
2795  M (1,1)= mulMod2 (bufFactors [0], bufFactors [1] [0], yToL);
2796  Pi [0]= M (1,1) + mulMod2 (bufFactors [0], bufFactors[1] [1], yToL)*v;
2797  }
2798  else
2799  {
2800  M (1,1)= mulMod2 (bufFactors [0], bufFactors [1], yToL);
2801  Pi [0]= M (1,1);
2802  }
2803 
2804  for (i= 1; i < Pi.size(); i++)
2805  {
2806  if (degree (Pi[i-1], v) > 0 && degree (bufFactors [i+1], v) > 0)
2807  {
2808  M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors[i+1] [0], yToL);
2809  Pi [i]= M (1,i+1) + (mulMod2 (Pi[i-1] [1], bufFactors[i+1] [0], yToL) +
2810  mulMod2 (Pi[i-1] [0], bufFactors [i+1] [1], yToL))*v;
2811  }
2812  else if (degree (Pi[i-1], v) > 0)
2813  {
2814  M (1,i+1)= mulMod2 (Pi[i-1] [0], bufFactors [i+1], yToL);
2815  Pi [i]= M(1,i+1) + mulMod2 (Pi[i-1] [1], bufFactors[i+1], yToL)*v;
2816  }
2817  else if (degree (bufFactors[i+1], v) > 0)
2818  {
2819  M (1,i+1)= mulMod2 (Pi[i-1], bufFactors [i+1] [0], yToL);
2820  Pi [i]= M (1,i+1) + mulMod2 (Pi[i-1], bufFactors[i+1] [1], yToL)*v;
2821  }
2822  else
2823  {
2824  M (1,i+1)= mulMod2 (Pi [i-1], bufFactors [i+1], yToL);
2825  Pi [i]= M (1,i+1);
2826  }
2827  }
2828 
2829  CFList products;
2830  bufF= mod (F, Variable (3));
2831  TIMING_START (product1);
2832  for (CFListIterator k= factors; k.hasItem(); k++)
2833  products.append (bufF/k.getItem());
2834  TIMING_END_AND_PRINT (product1, "time for product1 in 23: ");
2835 
2836  CFList MOD= CFList (power (v, liftBound));
2837  MOD.insert (yToL);
2838  for (int d= 1; d < liftBound; d++)
2839  {
2840  nonMonicHenselStep (F, factors, bufFactors, diophant, M, Pi, products, d,
2841  MOD, bad);
2842  if (bad)
2843  return CFList();
2844  }
2845 
2846  CFList result;
2847  for (i= 0; i < factors.length(); i++)
2848  result.append (bufFactors[i]);
2849  return result;
2850 }
2851 #endif
2852 
2853 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // nonMonicHenselStep
2854 CFList
2855 nonMonicHenselLift (const CFList& F, const CFList& factors, const CFList& LCs,
2856  CFList& diophant, CFArray& Pi, CFMatrix& M, int lOld,
2857  int& lNew, const CFList& MOD, bool& noOneToOne
2858  )
2859 {
2860 
2861  int k= 0;
2862  CFArray bufFactors= CFArray (factors.length());
2863  CFListIterator j= LCs;
2864  for (CFListIterator i= factors; i.hasItem(); i++, j++, k++)
2865  bufFactors [k]= replaceLC (i.getItem(), j.getItem());
2866 
2867  Variable y= F.getLast().mvar();
2868  Variable x= F.getFirst().mvar();
2869  CanonicalForm xToLOld= power (x, lOld);
2870 
2871  Pi [0]= mod (Pi[0], xToLOld);
2872  M (1, 1)= Pi [0];
2873 
2874  if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
2875  Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
2876  mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
2877  else if (degree (bufFactors[0], y) > 0)
2878  Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
2879  else if (degree (bufFactors[1], y) > 0)
2880  Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
2881 
2882  for (int i= 1; i < Pi.size(); i++)
2883  {
2884  Pi [i]= mod (Pi [i], xToLOld);
2885  M (1, i + 1)= Pi [i];
2886 
2887  if (degree (Pi[i-1], y) > 0 && degree (bufFactors [i+1], y) > 0)
2888  Pi [i] += (mulMod (Pi[i-1] [1], bufFactors[i+1] [0], MOD) +
2889  mulMod (Pi[i-1] [0], bufFactors [i+1] [1], MOD))*y;
2890  else if (degree (Pi[i-1], y) > 0)
2891  Pi [i] += mulMod (Pi[i-1] [1], bufFactors[i+1], MOD)*y;
2892  else if (degree (bufFactors[i+1], y) > 0)
2893  Pi [i] += mulMod (Pi[i-1], bufFactors[i+1] [1], MOD)*y;
2894  }
2895 
2896  CFList products;
2897  CanonicalForm quot, bufF= F.getFirst();
2898 
2899  TIMING_START (product1);
2900  for (int i= 0; i < bufFactors.size(); i++)
2901  {
2902  if (degree (bufFactors[i], y) > 0)
2903  {
2904  if (!fdivides (bufFactors[i] [0], bufF, quot))
2905  {
2906  noOneToOne= true;
2907  return factors;
2908  }
2909  products.append (quot);
2910  }
2911  else
2912  {
2913  if (!fdivides (bufFactors[i], bufF, quot))
2914  {
2915  noOneToOne= true;
2916  return factors;
2917  }
2918  products.append (quot);
2919  }
2920  }
2921  TIMING_END_AND_PRINT (product1, "time for product1 in further hensel:" );
2922 
2923  for (int d= 1; d < lNew; d++)
2924  {
2925  nonMonicHenselStep (F.getLast(), factors, bufFactors, diophant, M, Pi,
2926  products, d, MOD, noOneToOne);
2927  if (noOneToOne)
2928  return CFList();
2929  }
2930 
2931  CFList result;
2932  for (k= 0; k < factors.length(); k++)
2933  result.append (bufFactors[k]);
2934  return result;
2935 }
2936 #endif
2937 
2938 #if defined(HAVE_NTL) || defined(HAVE_FLINT) // nonMonicHenselLift23
2939 CFList
2940 nonMonicHenselLift (const CFList& eval, const CFList& factors,
2941  CFList* const& LCs, CFList& diophant, CFArray& Pi,
2942  int* liftBound, int length, bool& noOneToOne
2943  )
2944 {
2945  CFList bufDiophant= diophant;
2946  CFList buf= factors;
2947  CFArray bufPi= Pi;
2948  CFMatrix M= CFMatrix (liftBound[1], factors.length() - 1);
2949  int k= 0;
2950 
2951  TIMING_START (hensel23);
2952  CFList result=
2953  nonMonicHenselLift23 (eval.getFirst(), factors, LCs [0], diophant, bufPi,
2954  liftBound[1], liftBound[0], noOneToOne);
2955  TIMING_END_AND_PRINT (hensel23, "time for 23: ");
2956 
2957  if (noOneToOne)
2958  return CFList();
2959 
2960  if (eval.length() == 1)
2961  return result;
2962 
2963  k++;
2964  CFList MOD;
2965  for (int i= 0; i < 2; i++)
2966  MOD.append (power (Variable (i + 2), liftBound[i]));
2967 
2969  CFList bufEval;
2970  bufEval.append (j.getItem());
2971  j++;
2972 
2973  for (int i= 2; i <= length && j.hasItem(); i++, j++, k++)
2974  {
2975  bufEval.append (j.getItem());
2976  M= CFMatrix (liftBound[i], factors.length() - 1);
2977  TIMING_START (hensel);
2978  result= nonMonicHenselLift (bufEval, result, LCs [i-1], diophant, bufPi, M,
2979  liftBound[i-1], liftBound[i], MOD, noOneToOne);
2980  TIMING_END_AND_PRINT (hensel, "time for further hensel: ");
2981  if (noOneToOne)
2982  return result;
2983  MOD.append (power (Variable (i + 2), liftBound[i]));
2984  bufEval.removeFirst();
2985  }
2986 
2987  return result;
2988 }
2989 #endif
2990 #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 FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
CanonicalForm den(const CanonicalForm &f)
CanonicalForm Lc(const CanonicalForm &f)
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:4078
g
Definition: cfModGcd.cc:4090
bool equal
Definition: cfModGcd.cc:4126
CanonicalForm test
Definition: cfModGcd.cc:4096
CanonicalForm b
Definition: cfModGcd.cc:4103
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:31
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:127
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:2855
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:1369
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:1785
CFList nonMonicHenselLift23(const CanonicalForm &F, const CFList &factors, const CFList &LCs, CFList &diophant, CFArray &Pi, int liftBound, int bivarLiftBound, bool &bad)
Definition: facHensel.cc:2751
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:1274
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:2632
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:2553
CFList diophantine(const CanonicalForm &F, const CFList &factors)
Definition: facHensel.cc:1062
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:2154
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:2313
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:1343
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:2568
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:1852
CFList multiRecDiophantine(const CanonicalForm &F, const CFList &factors, const CFList &recResult, const CFList &M, int d)
Definition: facHensel.cc:1470
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:1559
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:1932
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:1070
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:1827
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