FORM 4.3
polyfact.cc
Go to the documentation of this file.
1
5/* #[ License : */
6/*
7 * Copyright (C) 1984-2022 J.A.M. Vermaseren
8 * When using this file you are requested to refer to the publication
9 * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
10 * This is considered a matter of courtesy as the development was paid
11 * for by FOM the Dutch physics granting agency and we would like to
12 * be able to track its scientific use to convince FOM of its value
13 * for the community.
14 *
15 * This file is part of FORM.
16 *
17 * FORM is free software: you can redistribute it and/or modify it under the
18 * terms of the GNU General Public License as published by the Free Software
19 * Foundation, either version 3 of the License, or (at your option) any later
20 * version.
21 *
22 * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
23 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
24 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
25 * details.
26 *
27 * You should have received a copy of the GNU General Public License along
28 * with FORM. If not, see <http://www.gnu.org/licenses/>.
29 */
30/* #] License : */
31/*
32 #[ include :
33*/
34
35#include "poly.h"
36#include "polygcd.h"
37#include "polyfact.h"
38
39#include <cmath>
40#include <vector>
41#include <iostream>
42#include <algorithm>
43#include <climits>
44
45//#define DEBUG
46
47#ifdef DEBUG
48#include "mytime.h"
49#endif
50
51using namespace std;
52
53/*
54 #] include :
55 #[ tostring :
56*/
57
58// Turns a factorized_poly into a readable string
59const string factorized_poly::tostring () const {
60
61 // empty
62 if (factor.size()==0)
63 return "no_factors";
64
65 string res;
66
67 // polynomial
68 for (int i=0; i<(int)factor.size(); i++) {
69 if (i>0) res += "*";
70 res += "(";
71 res += poly(factor[i],0,1).to_string();
72 res += ")";
73 if (power[i]>1) {
74 res += "^";
75 char tmp[100];
76 snprintf (tmp,100,"%i",power[i]);
77 res += tmp;
78 }
79 }
80
81 // modulo p^n
82 if (factor[0].modp>0) {
83 res += " (mod ";
84 char tmp[12];
85 snprintf (tmp,12,"%i",factor[0].modp);
86 res += tmp;
87 if (factor[0].modn > 1) {
88 snprintf (tmp,12,"%i",factor[0].modn);
89 res += "^";
90 res += tmp;
91 }
92 res += ")";
93 }
94
95 return res;
96}
97
98/*
99 #] tostring :
100 #[ ostream operator :
101*/
102
103// ostream operator for outputting a factorized_poly
104ostream& operator<< (ostream &out, const factorized_poly &a) {
105 return out << a.tostring();
106}
107
108// ostream operator for outputting a vector<T>
109template<class T> ostream& operator<< (ostream &out, const vector<T> &v) {
110 out<<"{";
111 for (int i=0; i<(int)v.size(); i++) {
112 if (i>0) out<<",";
113 out<<v[i];
114 }
115 out<<"}";
116 return out;
117}
118
119/*
120 #] ostream operator :
121 #[ add_factor :
122*/
123
124// adds a factor f^p to a factorization
125void factorized_poly::add_factor(const poly &f, int p) {
126 factor.push_back(f);
127 power.push_back(p);
128}
129
130/*
131 #] add_factor :
132 #[ extended_gcd_Euclidean_lifted :
133*/
134
152const vector<poly> polyfact::extended_gcd_Euclidean_lifted (const poly &a, const poly &b) {
153
154#ifdef DEBUGALL
155 cout << "*** [" << thetime() << "] CALL: extended_Euclidean_lifted("<<a<<","<<b<<")"<<endl;
156#endif
157
158 POLY_GETIDENTITY(a);
159
160 // Calculate s,t,gcd (mod p) with the Extended Euclidean Algorithm.
161 poly s(a,a.modp,1);
162 poly t(b,b.modp,1);
163 poly sa(BHEAD 1,a.modp,1);
164 poly sb(BHEAD 0,b.modp,1);
165 poly ta(BHEAD 0,a.modp,1);
166 poly tb(BHEAD 1,b.modp,1);
167
168 while (!t.is_zero()) {
169 poly x(s/t);
170 swap(s -=x*t , t);
171 swap(sa-=x*ta, ta);
172 swap(sb-=x*tb, tb);
173 }
174
175 // Normalize the result.
176 sa /= s.integer_lcoeff();
177 sb /= s.integer_lcoeff();
178
179 // Lift the result to modolu p^n with p-adic Newton's iteration.
180 poly samodp(sa);
181 poly sbmodp(sb);
182 poly term(BHEAD 1);
183
184 sa.setmod(0,1);
185 sb.setmod(0,1);
186
187 poly amodp(a,a.modp,1);
188 poly bmodp(b,a.modp,1);
189 poly error(poly(BHEAD 1) - sa*a - sb*b);
190 poly p(BHEAD a.modp);
191
192 for (int n=1; n<a.modn && !error.is_zero(); n++) {
193 error /= p;
194 term *= p;
195 poly errormodp(error,a.modp,1);
196 poly dsa((samodp * errormodp) % bmodp);
197 // equivalent, but faster than the symmetric
198 // poly dsb((sbmodp * errormodp) % amodp);
199 poly dsb((errormodp - dsa*amodp) / bmodp);
200 sa += term * dsa;
201 sb += term * dsb;
202 error -= a*dsa + b*dsb;
203 }
204
205 sa.setmod(a.modp,a.modn);
206 sb.setmod(a.modp,a.modn);
207
208 // Output the result
209 vector<poly> res;
210 res.push_back(sa);
211 res.push_back(sb);
212
213#ifdef DEBUGALL
214 cout << "*** [" << thetime() << "] RES : extended_Euclidean_lifted("<<a<<","<<b<<") = "<<res<<endl;
215#endif
216
217 return res;
218}
219
220/*
221 #] extended_gcd_Euclidean_lifted :
222 #[ solve_Diophantine_univariate :
223*/
224
256const vector<poly> polyfact::solve_Diophantine_univariate (const vector<poly> &a, const poly &b) {
257
258#ifdef DEBUGALL
259 cout << "*** [" << thetime() << "] CALL: solve_Diophantine_univariate(" <<a<<","<<b<<")"<<endl;
260#endif
261
262 POLY_GETIDENTITY(b);
263
264 vector<poly> s(1,b);
265
266 for (int i=0; i+1<(int)a.size(); i++) {
267 poly A(BHEAD 1,b.modp,b.modn);
268 for (int j=i+1; j<(int)a.size(); j++) A *= a[j];
269
270 vector<poly> t(extended_gcd_Euclidean_lifted(a[i],A));
271 poly prev(s.back());
272 s.back() = t[1] * prev % a[i];
273 s.push_back(t[0] * prev % A);
274 }
275
276#ifdef DEBUGALL
277 cout << "*** [" << thetime() << "] RES : solve_Diophantine_univariate(" <<a<<","<<b<<") = "<<s<<endl;
278#endif
279
280 return s;
281}
282
283/*
284 #] solve_Diophantine_univariate :
285 #[ solve_Diophantine_multivariate :
286*/
287
319const vector<poly> polyfact::solve_Diophantine_multivariate (const vector<poly> &a, const poly &b, const vector<int> &x, const vector<int> &c, int d) {
320
321#ifdef DEBUGALL
322 cout << "*** [" << thetime() << "] CALL: solve_Diophantine_multivariate(" <<a<<","<<b<<","<<x<<","<<c<<","<<d<<")"<<endl;
323#endif
324
325 POLY_GETIDENTITY(b);
326
327 if (b.is_zero()) return vector<poly>(a.size(),poly(BHEAD 0));
328
329 if (x.size() == 1) return solve_Diophantine_univariate(a,b);
330
331 // Reduce the polynomial with the homomorphism <xm-c{m-1}>
332 poly simple(poly::simple_poly(BHEAD x.back(),c.back()));
333
334 vector<poly> ared (a);
335 for (int i=0; i<(int)ared.size(); i++)
336 ared[i] %= simple;
337
338 poly bred(b % simple);
339 vector<int> xred(x.begin(),x.end()-1);
340 vector<int> cred(c.begin(),c.end()-1);
341
342 // Solve the equation in one less variable
343 vector<poly> s(solve_Diophantine_multivariate(ared,bred,xred,cred,d));
344 if (s == vector<poly>()) return vector<poly>();
345
346 // Cache the Ai = product(aj | j!=i).
347 vector<poly> A(a.size(), poly(BHEAD 1,b.modp,b.modn));
348 for (int i=0; i<(int)a.size(); i++)
349 for (int j=0; j<(int)a.size(); j++)
350 if (i!=j) A[i] *= a[j];
351
352 // Add the powers (xm-c{m-1})^k with ideal-adic Newton iteration.
353 poly term(BHEAD 1,b.modp,b.modn);
354
355 poly error(b);
356 for (int i=0; i<(int)A.size(); i++)
357 error -= s[i] * A[i];
358
359 for (int deg=1; deg<=d; deg++) {
360
361 if (error.is_zero()) break;
362
363 error /= simple;
364 term *= simple;
365
366 vector<poly> ds(solve_Diophantine_multivariate(ared, error%simple, xred, cred, d));
367 if (ds == vector<poly>()) return vector<poly>();
368
369 for (int i=0; i<(int)s.size(); i++) {
370 s[i] += ds[i] * term;
371 error -= ds[i] * A[i];
372 }
373 }
374
375 if (!error.is_zero()) return vector<poly>();
376
377#ifdef DEBUGALL
378 cout << "*** [" << thetime() << "] RES : solve_Diophantine_multivariate(" <<a<<","<<b<<","<<x<<","<<c<<","<<d<<") = "<<s<<endl;
379#endif
380
381 return s;
382}
383
384/*
385 #] solve_Diophantine_multivariate :
386 #[ lift_coefficients :
387*/
388
421const vector<poly> polyfact::lift_coefficients (const poly &_A, const vector<poly> &_a) {
422
423#ifdef DEBUG
424 cout << "*** [" << thetime() << "] CALL: lift_coefficients("<<_A<<","<<_a<<")"<<endl;
425#endif
426
427 POLY_GETIDENTITY(_A);
428
429 poly A(_A);
430 vector<poly> a(_a);
431 poly term(BHEAD 1);
432
433 int x = A.first_variable();
434
435 // Replace the leading term of all factors with lterm(A) mod p
436 poly lead(A.integer_lcoeff());
437 for (int i=0; i<(int)a.size(); i++) {
438 a[i] *= lead / a[i].integer_lcoeff();
439 if (i>0) A*=lead;
440 }
441
442 // Solve Diophantine equation
443 vector<poly> s(solve_Diophantine_univariate(a,poly(BHEAD 1,A.modp,1)));
444
445 // Replace the leading term of all factors with lterm(A) mod p^n
446 for (int i=0; i<(int)a.size(); i++) {
447 a[i].setmod(A.modp,A.modn);
448 a[i] += (lead - a[i].integer_lcoeff()) * poly::simple_poly(BHEAD x,0,a[i].degree(x));
449 }
450
451 // Calculate the error, express it in terms of ai and add corrections.
452 for (int k=2; k<=A.modn; k++) {
453 term *= poly(BHEAD A.modp);
454
455 poly error(BHEAD -1);
456 for (int i=0; i<(int)a.size(); i++) error *= a[i];
457 error += A;
458 if (error.is_zero()) break;
459
460 error /= term;
461 error.setmod(A.modp,1);
462
463 for (int i=0; i<(int)a.size(); i++)
464 a[i] += term * (error * s[i] % a[i]);
465 }
466
467 // Fix leading coefficients by dividing out integer contents.
468 for (int i=0; i<(int)a.size(); i++)
469 a[i] /= polygcd::integer_content(poly(a[i],0,1));
470
471#ifdef DEBUG
472 cout << "*** [" << thetime() << "] RES : lift_coefficients("<<_A<<","<<_a<<") = "<<a<<endl;
473#endif
474
475 return a;
476}
477
478/*
479 #] lift_coefficients :
480 #[ predetermine :
481*/
482
497void polyfact::predetermine (int dep, const vector<vector<int> > &state, vector<vector<vector<int> > > &terms, vector<int> &term, int sumdeg) {
498 // store the term
499 if (dep == (int)state.size()) {
500 terms[sumdeg].push_back(term);
501 return;
502 }
503
504 // recursively create new terms
505 term.push_back(0);
506
507 for (int deg=0; sumdeg+deg<(int)state[dep].size(); deg++)
508 if (state[dep][deg] > 0) {
509 term.back() = deg;
510 predetermine(dep+1, state, terms, term, sumdeg+deg);
511 }
512
513 term.pop_back();
514}
515
516/*
517 #] predetermine :
518 #[ lift_variables :
519*/
520
558const vector<poly> polyfact::lift_variables (const poly &A, const vector<poly> &_a, const vector<int> &x, const vector<int> &c, const vector<poly> &lc) {
559
560#ifdef DEBUG
561 cout << "*** [" << thetime() << "] CALL: lift_variables("<<A<<","<<_a<<","<<x<<","<<c<<","<<lc<<")\n";
562#endif
563
564 // If univariate, don't lift
565 if (x.size()<=1) return _a;
566
567 POLY_GETIDENTITY(A);
568
569 vector<poly> a(_a);
570
571 // First method: predetermine coefficients
572
573 // check feasibility, otherwise it tries too many possibilities
574 int cnt = POLYFACT_MAX_PREDETERMINATION;
575 for (int i=0; i<(int)a.size(); i++) {
576 if (a[i].number_of_terms() == 0) return vector<poly>();
577 cnt /= a[i].number_of_terms();
578 }
579
580 if (cnt>0) {
581 // state[n][d]: coefficient of x^d in a[n] is
582 // 0: non-existent, 1: undetermined, 2: determined
583 int D = A.degree(x[0]);
584 vector<vector<int> > state(a.size(), vector<int>(D+1, 0));
585
586 for (int i=0; i<(int)a.size(); i++)
587 for (int j=1; j<a[i][0]; j+=a[i][j])
588 state[i][a[i][j+1+x[0]]] = j==1 ? 2 : 1;
589
590 // collect all products of terms
591 vector<vector<vector<int> > > terms(D+1);
592 vector<int> term;
593 predetermine(0,state,terms,term);
594
595 // count the number of undetermined coefficients
596 vector<int> num_undet(terms.size(),0);
597 for (int i=0; i<(int)terms.size(); i++)
598 for (int j=0; j<(int)terms[i].size(); j++)
599 for (int k=0; k<(int)terms[i][j].size(); k++)
600 if (state[k][terms[i][j][k]] == 1) num_undet[i]++;
601
602 // replace the current leading coefficients by the correct ones
603 for (int i=0; i<(int)a.size(); i++)
604 a[i] += (lc[i] - a[i].lcoeff_univar(x[0])) * poly::simple_poly(BHEAD x[0],0,a[i].degree(x[0]));
605
606 bool changed;
607 do {
608 changed = false;
609
610 for (int i=0; i<(int)terms.size(); i++) {
611 // is only one coefficient in a equation is undetermined, solve
612 // the equation to determine this coefficient
613 if (num_undet[i] == 1) {
614 // generate equation
615 poly lhs(BHEAD 0), rhs(A.coefficient(x[0],i), A.modp, A.modn);
616 int which_idx=-1, which_deg=-1;
617 for (int j=0; j<(int)terms[i].size(); j++) {
618 poly coeff(BHEAD 1, A.modp, A.modn);
619 bool undet=false;
620 for (int k=0; k<(int)terms[i][j].size(); k++) {
621 if (state[k][terms[i][j][k]] == 1) {
622 undet = true;
623 which_idx=k;
624 which_deg=terms[i][j][k];
625 }
626 else
627 coeff *= a[k].coefficient(x[0], terms[i][j][k]);
628 }
629 if (undet)
630 lhs = coeff;
631 else
632 rhs -= coeff;
633 }
634 // solve equation
635 if (A.modn > 1) rhs.setmod(0,1);
636 if (lhs.is_zero() || !(rhs%lhs).is_zero()) return vector<poly>();
637 a[which_idx] += (rhs / lhs - a[which_idx].coefficient(x[0],which_deg)) * poly::simple_poly(BHEAD x[0],0,which_deg);
638 state[which_idx][which_deg] = 2;
639
640 // update number of undetermined coefficients
641 for (int j=0; j<(int)terms.size(); j++)
642 for (int k=0; k<(int)terms[j].size(); k++)
643 if (terms[j][k][which_idx] == which_deg)
644 num_undet[j]--;
645
646 changed = true;
647 }
648 }
649 }
650 while (changed);
651
652 // if this is the complete result, skip lifting
653 poly check(BHEAD 1, A.modn>1?0:A.modp, 1);
654 for (int i=0; i<(int)a.size(); i++)
655 check *= a[i];
656
657 if (check == A) return a;
658 }
659
660 // Second method: Hensel lifting
661
662 // Calculate A and lc's modulo Ii = <xi-c{i-1],...,xm-c{m-1}> (for i=2,...,m)
663 vector<poly> simple(x.size(), poly(BHEAD 0));
664 for (int i=(int)x.size()-2; i>=0; i--)
665 simple[i] = poly::simple_poly(BHEAD x[i+1],c[i],1);
666
667 // Calculate the maximum degree of A in x2,...,xm
668 int maxdegA=0;
669 for (int i=1; i<(int)x.size(); i++)
670 maxdegA = MaX(maxdegA, A.degree(x[i]));
671
672 // Iteratively add the variables x2,...,xm
673 for (int xi=1; xi<(int)x.size(); xi++) {
674 // replace the current leading coefficients by the correct ones
675 for (int i=0; i<(int)a.size(); i++)
676 a[i] += (lc[i] - a[i].lcoeff_univar(x[0])) * poly::simple_poly(BHEAD x[0],0,a[i].degree(x[0]));
677
678 vector<poly> anew(a);
679 for (int i=0; i<(int)anew.size(); i++)
680 for (int j=xi-1; j<(int)c.size(); j++)
681 anew[i] %= simple[j];
682
683 vector<int> xnew(x.begin(), x.begin()+xi);
684 vector<int> cnew(c.begin(), c.begin()+xi-1);
685 poly term(BHEAD 1,A.modp,A.modn);
686
687 // Iteratively add the powers xi^k
688 for (int deg=1, maxdeg=A.degree(x[xi]); deg<=maxdeg; deg++) {
689
690 term *= simple[xi-1];
691
692 // Calculate the error, express it in terms of ai and add corrections.
693 poly error(BHEAD -1,A.modp,A.modn);
694 for (int i=0; i<(int)a.size(); i++) error *= a[i];
695 error += A;
696 for (int i=xi; i<(int)c.size(); i++) error %= simple[i];
697
698 if (error.is_zero()) break;
699
700 error /= term;
701 error %= simple[xi-1];
702
703 vector<poly> s(solve_Diophantine_multivariate(anew,error,xnew,cnew,maxdegA));
704 if (s == vector<poly>()) return vector<poly>();
705
706 for (int i=0; i<(int)a.size(); i++)
707 a[i] += s[i] * term;
708 }
709
710 // check whether PRODUCT(a[i]) = A mod <xi-c{i-1},...,xm-c{m-1]> over the integers or ZZ/p
711 poly check(BHEAD -1, A.modn>1?0:A.modp, 1);
712 for (int i=0; i<(int)a.size(); i++) check *= a[i];
713 check += A;
714 for (int i=xi; i<(int)c.size(); i++) check %= simple[i];
715 if (!check.is_zero()) return vector<poly>();
716 }
717
718#ifdef DEBUG
719 cout << "*** [" << thetime() << "] RES : lift_variables("<<A<<","<<_a<<","<<x<<","<<c<<","<<lc<<") = " << a << endl;
720#endif
721
722 return a;
723}
724
725/*
726 #] lift_variables :
727 #[ choose_prime :
728*/
729
741WORD polyfact::choose_prime (const poly &a, const vector<int> &x, WORD p) {
742
743#ifdef DEBUG
744 cout << "*** [" << thetime() << "] CALL: choose_prime("<<a<<","<<x<<","<<p<<")"<<endl;
745#endif
746
747 POLY_GETIDENTITY(a);
748
749 poly icont_lcoeff(polygcd::integer_content(a.lcoeff_univar(x[0])));
750
751 if (p==0) p = POLYFACT_FIRST_PRIME;
752
753 poly icont_lcoeff_modp(BHEAD 0);
754
755 do {
756
757 bool is_prime;
758
759 do {
760 p += 2;
761 is_prime = true;
762 for (int d=2; d*d<=p; d++)
763 if (p%d==0) { is_prime=false; break; }
764 }
765 while (!is_prime);
766
767 icont_lcoeff_modp = icont_lcoeff;
768 icont_lcoeff_modp.setmod(p,1);
769 }
770 while (icont_lcoeff_modp.is_zero());
771
772#ifdef DEBUG
773 cout << "*** [" << thetime() << "] RES : choose_prime("<<a<<","<<x<<",?) = "<<p<<endl;
774#endif
775
776 return p;
777}
778
779/*
780 #] choose_prime :
781 #[ choose_prime_power :
782*/
783
798WORD polyfact::choose_prime_power (const poly &a, WORD p) {
799
800#ifdef DEBUG
801 cout << "*** [" << thetime() << "] CALL: choose_prime_power("<<a<<","<<p<<")"<<endl;
802#endif
803
804 POLY_GETIDENTITY(a);
805
806 // analyse the polynomial for calculating the bound
807 double maxdegree=0, maxlogcoeff=0, numterms=0;
808
809 for (int i=1; i<a[0]; i+=a[i]) {
810 for (int j=0; j<AN.poly_num_vars; j++)
811 maxdegree = MaX(maxdegree, a[i+1+j]);
812
813 maxlogcoeff = MaX(maxlogcoeff,
814 log(1.0+(UWORD)a[i+a[i]-2]) + // most significant digit + 1
815 BITSINWORD*log(2.0)*(ABS(a[i+a[i]-1])-1)); // number of digits
816 numterms++;
817 }
818
819 WORD res = (WORD)ceil((log((sqrt(5.0)+1)/2)*maxdegree + maxlogcoeff + 0.5*log(numterms)) / log((double)p));
820
821#ifdef DEBUG
822 cout << "*** [" << thetime() << "] CALL: choose_prime_power("<<a<<","<<p<<") = "<<res<<endl;
823#endif
824
825 return res;
826}
827
828/*
829 #] choose_prime_power :
830 #[ choose_ideal :
831*/
832
860const vector<int> polyfact::choose_ideal (const poly &a, int p, const factorized_poly &lc, const vector<int> &x) {
861
862#ifdef DEBUG
863 cout << "*** [" << thetime() << "] CALL: polyfact::choose_ideal("
864 <<a<<","<<p<<","<<lc<<","<<x<<")"<<endl;
865#endif
866
867 if (x.size()==1) return vector<int>();
868
869 POLY_GETIDENTITY(a);
870
871 vector<int> c(x.size()-1);
872
873 int dega = a.degree(x[0]);
874 poly amodI(a);
875
876 // choose random c
877 for (int i=0; i<(int)c.size(); i++) {
878 c[i] = 1 + wranf(BHEAD0) % ((p-1) / POLYFACT_IDEAL_FRACTION);
879 amodI %= poly::simple_poly(BHEAD x[i+1],c[i],1);
880 }
881
882 poly amodIp(amodI);
883 amodIp.setmod(p,1);
884
885 // check if leading coefficient is non-zero [equivalent to degree=old_degree]
886 if (amodIp.degree(x[0]) != dega)
887 return c = vector<int>();
888
889 // check if leading coefficient is squarefree [equivalent to gcd(a,a')==const]
890 if (!polygcd::gcd_Euclidean(amodIp, amodIp.derivative(x[0])).is_integer())
891 return c = vector<int>();
892
893 if (a.modp>0 && a.modn==1) return c;
894
895 // check for unique prime factors in each factor lc[i] of the leading coefficient
896 vector<poly> d(1, polygcd::integer_content(amodI));
897
898 for (int i=0; i<(int)lc.factor.size(); i++) {
899 // constant factor
900 if (i==0 && lc.factor[i].is_integer()) {
901 d[0] *= lc.factor[i];
902 continue;
903 }
904
905 // factor modulo I
906 poly q(lc.factor[i]);
907 for (int j=0; j<(int)c.size(); j++)
908 q %= poly::simple_poly(BHEAD x[j+1],c[j]);
909 if (q.sign() == -1) q *= poly(BHEAD -1);
910
911 // divide out common factors
912 for (int j=(int)d.size()-1; j>=0; j--) {
913 poly r(d[j]);
914 while (!r.is_one()) {
915 r = polygcd::integer_gcd(r,q);
916 q /= r;
917 }
918 }
919
920 // check whether there is some factor left
921 if (q.is_one()) return vector<int>();
922 d.push_back(q);
923 }
924
925#ifdef DEBUG
926 cout << "*** [" << thetime() << "] RES : polyfact::choose_ideal("
927 <<a<<","<<p<<","<<lc<<","<<x<<") = "<<c<<endl;
928#endif
929
930 return c;
931}
932
933/*
934 #] choose_ideal :
935 #[ squarefree_factors_Yun :
936*/
937
944const factorized_poly polyfact::squarefree_factors_Yun (const poly &_a) {
945
946 factorized_poly res;
947 poly a(_a);
948
949 int pow = 1;
950 int x = a.first_variable();
951
952 poly b(a.derivative(x));
953 poly c(polygcd::gcd(a,b));
954
955 while (true) {
956 a /= c;
957 b /= c;
958 b -= a.derivative(x);
959 if (b.is_zero()) break;
960 c = polygcd::gcd(a,b);
961 if (!c.is_one()) res.add_factor(c,pow);
962 pow++;
963 }
964
965 if (!a.is_one()) res.add_factor(a,pow);
966 return res;
967}
968
969/*
970 #] squarefree_factors_Yun :
971 #[ squarefree_factors_modp :
972*/
973
980const factorized_poly polyfact::squarefree_factors_modp (const poly &_a) {
981
982 factorized_poly res;
983 poly a(_a);
984
985 int pow = 1;
986 int x = a.first_variable();
987 poly b(a.derivative(x));
988
989 // poly contains terms of the form c(x)^n (n!=c*p)
990 if (!b.is_zero()) {
991 poly c(polygcd::gcd(a,b));
992 a /= c;
993
994 while (!a.is_one()) {
995 b = polygcd::gcd(a,c);
996 a /= b;
997 if (!a.is_one()) res.add_factor(a,pow);
998 pow++;
999 a = b;
1000 c /= a;
1001 }
1002
1003 a = c;
1004 }
1005
1006 // polynomial contains terms of the form c(x)^p
1007 if (!a.is_one()) {
1008 for (int i=1; i<a[1]; i+=a[i])
1009 a[i+1+x] /= a.modp;
1010 factorized_poly res2(squarefree_factors(a));
1011 for (int i=0; i<(int)res2.factor.size(); i++) {
1012 res.factor.push_back(res2.factor[i]);
1013 res.power.push_back(a.modp*res2.power[i]);
1014 }
1015 }
1016
1017 return res;
1018}
1019
1020/*
1021 #] squarefree_factors_modp :
1022 #[ squarefree_factors :
1023*/
1024
1045const factorized_poly polyfact::squarefree_factors (const poly &a) {
1046
1047#ifdef DEBUG
1048 cout << "*** [" << thetime() << "] CALL: squarefree_factors("<<a<<")\n";
1049#endif
1050
1051 if (a.is_one()) return factorized_poly();
1052
1053 factorized_poly res;
1054
1055 if (a.modp==0)
1056 res = squarefree_factors_Yun(a);
1057 else
1058 res = squarefree_factors_modp(a);
1059
1060#ifdef DEBUG
1061 cout << "*** [" << thetime() << "] RES : squarefree_factors("<<a<<") = "<<res<<"\n";
1062#endif
1063
1064 return res;
1065}
1066
1067/*
1068 #] squarefree_factors :
1069 #[ Berlekamp_Qmatrix :
1070*/
1071
1094const vector<vector<WORD> > polyfact::Berlekamp_Qmatrix (const poly &_a) {
1095
1096#ifdef DEBUG
1097 cout << "*** [" << thetime() << "] CALL: Berlekamp_Qmatrix("<<_a<<")\n";
1098#endif
1099
1100 if (_a.all_variables() == vector<int>())
1101 return vector<vector<WORD> >(0);
1102
1103 POLY_GETIDENTITY(_a);
1104
1105 poly a(_a);
1106 int x = a.first_variable();
1107 int n = a.degree(x);
1108 int p = a.modp;
1109
1110 poly lc(a.integer_lcoeff());
1111 a /= lc;
1112
1113 vector<vector<WORD> > Q(n, vector<WORD>(n));
1114
1115 // c is the vector of coefficients of the polynomial a
1116 vector<WORD> c(n+1,0);
1117 for (int j=1; j<a[0]; j+=a[j])
1118 c[a[j+1+x]] = a[j+1+AN.poly_num_vars] * a[j+2+AN.poly_num_vars];
1119
1120 // d is the vector of coefficients of x^i mod a, starting with i=0
1121 vector<WORD> d(n,0);
1122 d[0]=1;
1123
1124 for (int i=0; i<=(n-1)*p; i++) {
1125 // store the coefficients of x^(i*p) mod a
1126 if (i%p==0) Q[i/p] = d;
1127
1128 // transform d=x^i mod a into d=x^(i+1) mod a
1129 vector<WORD> e(n);
1130 for (int j=0; j<n; j++) {
1131 e[j] = (-(LONG)d[n-1]*c[j] + (j>0?d[j-1]:0)) % p;
1132 if (e[j]<0) e[j]+=p;
1133 }
1134
1135 d=e;
1136 }
1137
1138 // Q = Q - I
1139 for (int i=0; i<n; i++)
1140 Q[i][i] = (Q[i][i] - 1 + p) % p;
1141
1142 // Gaussian elimination
1143 for (int i=0; i<n; i++) {
1144 // Find pivot
1145 int ii=i; while (ii<n && Q[i][ii]==0) ii++;
1146 if (ii==n) continue;
1147
1148 for (int k=0; k<n; k++)
1149 swap(Q[k][ii],Q[k][i]);
1150
1151 // normalize row i, which becomes the pivot
1152 WORD mul;
1153 GetModInverses (Q[i][i], p, &mul, NULL);
1154 vector<int> idx;
1155
1156 for (int k=0; k<n; k++) if (Q[k][i] != 0) {
1157 // store indices of non-zero elements for sparse matrices
1158 idx.push_back(k);
1159 Q[k][i] = ((LONG)Q[k][i] * mul) % p;
1160 }
1161
1162 // reduce
1163 for (int j=0; j<n; j++)
1164 if (j!=i && Q[i][j]!=0) {
1165 LONG mul = Q[i][j];
1166 for (int k=0; k<(int)idx.size(); k++) {
1167 Q[idx[k]][j] = (Q[idx[k]][j] - mul*Q[idx[k]][i]) % p;
1168 if (Q[idx[k]][j] < 0) Q[idx[k]][j]+=p;
1169 }
1170 }
1171 }
1172
1173 for (int i=0; i<n; i++) {
1174 // Q = Q - I
1175 Q[i][i] = Q[i][i]-1;
1176
1177 // reduce all coefficients in the range 0,1,...,p-1
1178 for (int j=0; j<n; j++) {
1179 Q[i][j] = -Q[i][j]%p;
1180 if (Q[i][j]<0) Q[i][j]+=p;
1181 }
1182 }
1183
1184#ifdef DEBUG
1185 cout << "*** [" << thetime() << "] RES : Berlekamp_Qmatrix("<<_a<<") = "<<Q<<"\n";
1186#endif
1187
1188 return Q;
1189}
1190
1191/*
1192 #] Berlekamp_Qmatrix :
1193 #[ Berlekamp_find_factors :
1194*/
1195
1219const vector<poly> polyfact::Berlekamp_find_factors (const poly &_a, const vector<vector<WORD> > &_q) {
1220
1221#ifdef DEBUG
1222 cout << "*** [" << thetime() << "] CALL: Berlekamp_find_factors("<<_a<<","<<_q<<")\n";
1223#endif
1224
1225 if (_a.all_variables() == vector<int>()) return vector<poly>(1,_a);
1226
1227 POLY_GETIDENTITY(_a);
1228
1229 vector<vector<WORD> > q=_q;
1230
1231 int rank=0;
1232 for (int i=0; i<(int)q.size(); i++)
1233 if (q[i]!=vector<WORD>(q[i].size(),0)) rank++;
1234
1235 poly a(_a);
1236 int x = a.first_variable();
1237 int n = a.degree(x);
1238 int p = a.modp;
1239
1240 a /= a.integer_lcoeff();
1241
1242 // Vector of factors, represented as dense polynomials mod p
1243 vector<vector<WORD> > fac(1, vector<WORD>(n+1,0));
1244
1245 fac[0] = poly::to_coefficient_list(a);
1246 bool finished=false;
1247
1248 // Loop over the columns of q + constant, i.e., an exhaustive list of possible factors
1249 for (int i=1; i<n && !finished; i++) {
1250 if (q[i] == vector<WORD>(n,0)) continue;
1251
1252 for (int s=0; s<p && !finished; s++) {
1253 for (int j=0; j<(int)fac.size() && !finished; j++) {
1254 vector<WORD> c = polygcd::coefficient_list_gcd(fac[j], q[i], p);
1255
1256 // If a non-trivial factor is found, add it to the list
1257 if (c.size()!=1 && c.size()!=fac[j].size()) {
1258 fac.push_back(c);
1259 fac[j] = poly::coefficient_list_divmod(fac[j], c, p, 0);
1260 if ((int)fac.size() == rank) finished=true;
1261 }
1262 }
1263
1264 // Increase the constant term by one
1265 q[i][0] = (q[i][0]+1) % p;
1266 }
1267 }
1268
1269 // Convert the densely represented polynomials to sparse ones
1270 vector<poly> res(fac.size(),poly(BHEAD 0, p));
1271 for (int i=0; i<(int)fac.size(); i++)
1272 res[i] = poly::from_coefficient_list(BHEAD fac[i],x,p);
1273
1274#ifdef DEBUG
1275 cout << "*** [" << thetime() << "] RES : Berlekamp_find_factors("<<_a<<","<<_q<<") = "<<res<<"\n";
1276#endif
1277
1278 return res;
1279}
1280
1281/*
1282 #] Berlekamp_find_factors :
1283 #[ combine_factors :
1284*/
1285
1307const vector<poly> polyfact::combine_factors (const poly &a, const vector<poly> &f) {
1308
1309#ifdef DEBUG
1310 cout << "*** [" << thetime() << "] CALL: combine_factors("<<a<<","<<f<<")\n";
1311#endif
1312
1313 POLY_GETIDENTITY(a);
1314
1315 poly a0(a,0,1);
1316 vector<poly> res;
1317
1318 int num_used = 0;
1319 vector<bool> used(f.size(), false);
1320
1321 // Loop over all bitmasks with num=1,2,...,size(factors)/2 bits
1322 // set, that contain only unused factors
1323 for (int num=1; num<=(int)(f.size() - num_used)/2; num++) {
1324 vector<int> next(f.size() - num_used, 0);
1325 for (int i=0; i<num; i++) next[next.size()-1-i] = 1;
1326
1327 do {
1328 poly fac(BHEAD 1,a.modp,a.modn);
1329 for (int i=0, j=0; i<(int)f.size(); i++)
1330 if (!used[i] && next[j++]) fac *= f[i];
1331 fac /= fac.integer_lcoeff();
1332 fac *= a.integer_lcoeff();
1333 fac /= polygcd::integer_content(poly(fac,0,1));
1334
1335 if ((a0 % fac).is_zero()) {
1336 res.push_back(fac);
1337 for (int i=0, j=0; i<(int)f.size(); i++)
1338 if (!used[i]) used[i] = next[j++];
1339 num_used += num;
1340 num--;
1341 break;
1342 }
1343 }
1344 while (next_permutation(next.begin(), next.end()));
1345 }
1346
1347 // All unused factors together form one more factor
1348 if (num_used != (int)f.size()) {
1349 poly fac(BHEAD 1,a.modp,a.modn);
1350 for (int i=0; i<(int)f.size(); i++)
1351 if (!used[i]) fac *= f[i];
1352 fac /= fac.integer_lcoeff();
1353 fac *= a.integer_lcoeff();
1354 fac /= polygcd::integer_content(poly(fac,0,1));
1355 res.push_back(fac);
1356 }
1357
1358#ifdef DEBUG
1359 cout << "*** [" << thetime() << "] RES : combine_factors("<<a<<","<<f<<") = "<<res<<"\n";
1360#endif
1361
1362 return res;
1363}
1364
1365/*
1366 #] combine_factors :
1367 #[ factorize_squarefree :
1368*/
1369
1388const vector<poly> polyfact::factorize_squarefree (const poly &a, const vector<int> &x) {
1389
1390#ifdef DEBUG
1391 cout << "*** [" << thetime() << "] CALL: factorize_squarefree("<<a<<")\n";
1392#endif
1393
1394 POLY_GETIDENTITY(a);
1395
1396 WORD p=a.modp, n=a.modn;
1397
1398 try_again:
1399
1400 int bestp=0;
1401 int min_factors = INT_MAX;
1402 poly amodI(BHEAD 0), bestamodI(BHEAD 0);
1403 vector<int> c,d,bestc,bestd;
1404 vector<vector<WORD> > q,bestq;
1405
1406 // Factorize leading coefficient
1407 factorized_poly lc(factorize(a.lcoeff_univar(x[0])));
1408
1409 // Try a number of primes
1410 int prime_tries = 0;
1411
1412 while (prime_tries<POLYFACT_NUM_CONFIRMATIONS && min_factors>1) {
1413 if (a.modp == 0) {
1414 p = choose_prime(a,x,p);
1415 n = 0;
1416 if (a.degree(x[0]) % p == 0) continue;
1417
1418 // Univariate case: check whether the polynomial mod p is squarefree
1419 // Multivariate case: this check is done after choosing I (for efficiency)
1420 if (x.size()==1) {
1421 poly amodp(a,p,1);
1422 if (polygcd::gcd_Euclidean(amodp, amodp.derivative(x[0])).degree(x[0]) != 0)
1423 continue;
1424 }
1425 }
1426
1427 // Try a number of ideals
1428 if (x.size()>1)
1429 for (int ideal_tries=0; ideal_tries<POLYFACT_MAX_IDEAL_TRIES; ideal_tries++) {
1430 c = choose_ideal(a,p,lc,x);
1431 if (c.size()>0) break;
1432 }
1433
1434 if (x.size()==1 || c.size()>0) {
1435 amodI = a;
1436 for (int i=0; i<(int)c.size(); i++)
1437 amodI %= poly::simple_poly(BHEAD x[i+1],c[i]);
1438
1439 // Determine Q-matrix and its rank. Smaller rank is better.
1440 q = Berlekamp_Qmatrix(poly(amodI,p,1));
1441 int rank=0;
1442 for (int i=0; i<(int)q.size(); i++)
1443 if (q[i]!=vector<WORD>(q[i].size(),0)) rank++;
1444
1445 if (rank<min_factors) {
1446 bestp=p;
1447 bestc=c;
1448 bestq=q;
1449 bestamodI=amodI;
1450 min_factors = rank;
1451 prime_tries = 0;
1452 }
1453
1454 if (rank==min_factors)
1455 prime_tries++;
1456
1457#ifdef DEBUG
1458 cout << "*** [" << thetime() << "] (A) : factorize_squarefree("<<a<<
1459 ") try p=" << p << " #factors=" << rank << " (min="<<min_factors<<"x"<<prime_tries<<")" << endl;
1460#endif
1461 }
1462 }
1463
1464 p=bestp;
1465 c=bestc;
1466 q=bestq;
1467 amodI=bestamodI;
1468
1469 // Determine to which power of p to lift
1470 if (n==0) {
1471 n = choose_prime_power(amodI,p);
1472 n = MaX(n, choose_prime_power(a,p));
1473 }
1474
1475 amodI.setmod(p,n);
1476
1477#ifdef DEBUG
1478 cout << "*** [" << thetime() << "] (B) : factorize_squarefree("<<a<<
1479 ") chosen c = " << c << ", p^n = "<<p<<"^"<<n<<endl;
1480 cout << "*** [" << thetime() << "] (C) : factorize_squarefree("<<a<<
1481 ") #factors = " << min_factors << endl;
1482#endif
1483
1484 // Find factors
1485 vector<poly> f(Berlekamp_find_factors(poly(amodI,p,1),q));
1486
1487 // Lift coefficients
1488 if (f.size() > 1) {
1489 f = lift_coefficients(amodI,f);
1490 if (f==vector<poly>()) {
1491#ifdef DEBUG
1492 cout << "factor_squarefree failed (lift_coeff step) : " << endl;
1493#endif
1494 goto try_again;
1495 }
1496
1497 // Combine factors
1498 if (a.modp==0) f = combine_factors(amodI,f);
1499 }
1500
1501 // Lift variables
1502 if (f.size() == 1)
1503 f = vector<poly>(1, a);
1504
1505 if (x.size() > 1 && f.size() > 1) {
1506
1507 // The correct leading coefficients of the factors can be
1508 // reconstructed from prime number factors of the leading
1509 // coefficients modulo I. This is possible since all factors of
1510 // the leading coefficient have unique prime factors for the ideal
1511 // I is chosen as such.
1512
1513 poly amodI(a);
1514 for (int i=0; i<(int)c.size(); i++)
1515 amodI %= poly::simple_poly(BHEAD x[i+1],c[i]);
1516 poly delta(polygcd::integer_content(amodI));
1517
1518 vector<poly> lcmodI(lc.factor.size(), poly(BHEAD 0));
1519 for (int i=0; i<(int)lc.factor.size(); i++) {
1520 lcmodI[i] = lc.factor[i];
1521 for (int j=0; j<(int)c.size(); j++)
1522 lcmodI[i] %= poly::simple_poly(BHEAD x[j+1],c[j]);
1523 }
1524
1525 vector<poly> correct_lc(f.size(), poly(BHEAD 1,p,n));
1526
1527 for (int j=0; j<(int)f.size(); j++) {
1528 poly lc_f(f[j].integer_lcoeff() * delta);
1529 WORD nlc_f = lc_f[lc_f[1]];
1530 poly quo(BHEAD 0),rem(BHEAD 0);
1531 WORD nquo,nrem;
1532
1533 for (int i=(int)lcmodI.size()-1; i>=0; i--) {
1534
1535 if (i==0 && lc.factor[i].is_integer()) continue;
1536
1537 do {
1538 DivLong((UWORD *)&lc_f[2+AN.poly_num_vars], nlc_f,
1539 (UWORD *)&lcmodI[i][2+AN.poly_num_vars], lcmodI[i][lcmodI[i][1]],
1540 (UWORD *)&quo[0], &nquo,
1541 (UWORD *)&rem[0], &nrem);
1542
1543 if (nrem == 0) {
1544 correct_lc[j] *= lc.factor[i];
1545 lc_f.termscopy(&quo[0], 2+AN.poly_num_vars, ABS(nquo));
1546 nlc_f = nquo;
1547 }
1548 }
1549 while (nrem == 0);
1550 }
1551 }
1552
1553 for (int i=0; i<(int)correct_lc.size(); i++) {
1554 poly correct_modI(correct_lc[i]);
1555 for (int j=0; j<(int)c.size(); j++)
1556 correct_modI %= poly::simple_poly(BHEAD x[j+1],c[j]);
1557
1558 poly d(polygcd::integer_gcd(correct_modI, f[i].integer_lcoeff()));
1559 correct_lc[i] *= f[i].integer_lcoeff() / d;
1560 delta /= correct_modI / d;
1561 f[i] *= correct_modI / d;
1562 }
1563
1564 // increase n, because of multiplying with delta
1565 if (!delta.is_one()) {
1566 poly deltapow(BHEAD 1);
1567 for (int i=1; i<(int)correct_lc.size(); i++)
1568 deltapow *= delta;
1569 while (!deltapow.is_zero()) {
1570 deltapow /= poly(BHEAD p);
1571 n++;
1572 }
1573
1574 for (int i=0; i<(int)f.size(); i++) {
1575 f[i].modn = n;
1576 correct_lc[i].modn = n;
1577 }
1578 }
1579
1580 poly aa(a,p,n);
1581
1582 for (int i=0; i<(int)correct_lc.size(); i++) {
1583 correct_lc[i] *= delta;
1584 f[i] *= delta;
1585 if (i>0) aa *= delta;
1586 }
1587
1588 f = lift_variables(aa,f,x,c,correct_lc);
1589
1590 for (int i=0; i<(int)f.size(); i++)
1591 if (a.modp == 0)
1592 f[i] /= polygcd::integer_content(poly(f[i],0,1));
1593 else
1594 f[i] /= polygcd::content_univar(f[i], x[0]);
1595
1596 if (f==vector<poly>()) {
1597#ifdef DEBUG
1598 cout << "factor_squarefree failed (lift_var step)" << endl;
1599#endif
1600 goto try_again;
1601 }
1602 }
1603
1604 // set modulus of the factors correctly
1605 if (a.modp==0)
1606 for (int i=0; i<(int)f.size(); i++)
1607 f[i].setmod(0,1);
1608
1609 // Final check (not sure if this is necessary, but it doesn't hurt)
1610 poly check(BHEAD 1,a.modp,a.modn);
1611 for (int i=0; i<(int)f.size(); i++)
1612 check *= f[i];
1613
1614 if (check != a) {
1615#ifdef DEBUG
1616 cout << "factor_squarefree failed (final check) : " << f << endl;
1617#endif
1618 goto try_again;
1619 }
1620
1621#ifdef DEBUG
1622 cout << "*** [" << thetime() << "] RES : factorize_squarefree("<<a<<","<<x<<","<<c<<") = "<<f<<"\n";
1623#endif
1624
1625 return f;
1626}
1627
1628/*
1629 #] factorize_squarefree :
1630 #[ factorize :
1631*/
1632
1641const factorized_poly polyfact::factorize (const poly &a) {
1642
1643#ifdef DEBUG
1644 cout << "*** [" << thetime() << "] CALL: factorize("<<a<<")\n";
1645#endif
1646 vector<int> x = a.all_variables();
1647
1648 // No variables, so just one factor
1649 if (x.size() == 0) {
1650 factorized_poly res;
1651 if (a.is_one()) return res;
1652 res.add_factor(a,1);
1653 return res;
1654 }
1655
1656 // Remove content
1657 poly conta(polygcd::content_univar(a,x[0]));
1658
1659 factorized_poly faca(factorize(conta));
1660
1661 poly ppa(a / conta);
1662
1663 // Find a squarefree factorization
1664 factorized_poly b(squarefree_factors(ppa));
1665
1666#ifdef DEBUG
1667 cout << "*** [" << thetime() << "] ... : factorize("<<a<<") : SFF = "<<b<<"\n";
1668#endif
1669
1670 // Factorize each squarefree factor and build the "factorized_poly"
1671 for (int i=0; i<(int)b.factor.size(); i++) {
1672
1673 vector<poly> c(factorize_squarefree(b.factor[i], x));
1674
1675 for (int j=0; j<(int)c.size(); j++) {
1676 faca.factor.push_back(c[j]);
1677 faca.power.push_back(b.power[i]);
1678 }
1679 }
1680
1681#ifdef DEBUG
1682 cout << "*** [" << thetime() << "] RES : factorize("<<a<<") = "<<faca<<"\n";
1683#endif
1684
1685 return faca;
1686}
1687
1688/*
1689 #] factorize :
1690*/
Definition poly.h:49
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition reken.c:1466