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