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