59 const string factorized_poly::tostring ()
const {
68 for (
int i=0; i<(int)factor.size(); i++) {
71 res +=
poly(factor[i],0,1).to_string();
76 snprintf (tmp,100,
"%i",power[i]);
82 if (factor[0].modp>0) {
85 snprintf (tmp,12,
"%i",factor[0].modp);
87 if (factor[0].modn > 1) {
88 snprintf (tmp,12,
"%i",factor[0].modn);
105 return out << a.tostring();
109 template<
class T> ostream& operator<< (ostream &out, const vector<T> &v) {
111 for (
int i=0; i<(int)v.size(); i++) {
125 void factorized_poly::add_factor(
const poly &f,
int p) {
152 const vector<poly> polyfact::extended_gcd_Euclidean_lifted (
const poly &a,
const poly &b) {
155 cout <<
"*** [" << thetime() <<
"] CALL: extended_Euclidean_lifted("<<a<<
","<<b<<
")"<<endl;
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);
168 while (!t.is_zero()) {
176 sa /= s.integer_lcoeff();
177 sb /= s.integer_lcoeff();
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);
192 for (
int n=1; n<a.modn && !error.is_zero(); n++) {
195 poly errormodp(error,a.modp,1);
196 poly dsa((samodp * errormodp) % bmodp);
199 poly dsb((errormodp - dsa*amodp) / bmodp);
202 error -= a*dsa + b*dsb;
205 sa.setmod(a.modp,a.modn);
206 sb.setmod(a.modp,a.modn);
214 cout <<
"*** [" << thetime() <<
"] RES : extended_Euclidean_lifted("<<a<<
","<<b<<
") = "<<res<<endl;
256 const vector<poly> polyfact::solve_Diophantine_univariate (
const vector<poly> &a,
const poly &b) {
259 cout <<
"*** [" << thetime() <<
"] CALL: solve_Diophantine_univariate(" <<a<<
","<<b<<
")"<<endl;
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];
270 vector<poly> t(extended_gcd_Euclidean_lifted(a[i],A));
272 s.back() = t[1] * prev % a[i];
273 s.push_back(t[0] * prev % A);
277 cout <<
"*** [" << thetime() <<
"] RES : solve_Diophantine_univariate(" <<a<<
","<<b<<
") = "<<s<<endl;
319 const vector<poly> polyfact::solve_Diophantine_multivariate (
const vector<poly> &a,
const poly &b,
const vector<int> &x,
const vector<int> &c,
int d) {
322 cout <<
"*** [" << thetime() <<
"] CALL: solve_Diophantine_multivariate(" <<a<<
","<<b<<
","<<x<<
","<<c<<
","<<d<<
")"<<endl;
327 if (b.is_zero())
return vector<poly>(a.size(),
poly(BHEAD 0));
329 if (x.size() == 1)
return solve_Diophantine_univariate(a,b);
332 poly simple(poly::simple_poly(BHEAD x.back(),c.back()));
334 vector<poly> ared (a);
335 for (
int i=0; i<(int)ared.size(); i++)
338 poly bred(b % simple);
339 vector<int> xred(x.begin(),x.end()-1);
340 vector<int> cred(c.begin(),c.end()-1);
343 vector<poly> s(solve_Diophantine_multivariate(ared,bred,xred,cred,d));
344 if (s == vector<poly>())
return vector<poly>();
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];
353 poly term(BHEAD 1,b.modp,b.modn);
356 for (
int i=0; i<(int)A.size(); i++)
357 error -= s[i] * A[i];
359 for (
int deg=1; deg<=d; deg++) {
361 if (error.is_zero())
break;
366 vector<poly> ds(solve_Diophantine_multivariate(ared, error%simple, xred, cred, d));
367 if (ds == vector<poly>())
return vector<poly>();
369 for (
int i=0; i<(int)s.size(); i++) {
370 s[i] += ds[i] * term;
371 error -= ds[i] * A[i];
375 if (!error.is_zero())
return vector<poly>();
378 cout <<
"*** [" << thetime() <<
"] RES : solve_Diophantine_multivariate(" <<a<<
","<<b<<
","<<x<<
","<<c<<
","<<d<<
") = "<<s<<endl;
421 const vector<poly> polyfact::lift_coefficients (
const poly &_A,
const vector<poly> &_a) {
424 cout <<
"*** [" << thetime() <<
"] CALL: lift_coefficients("<<_A<<
","<<_a<<
")"<<endl;
427 POLY_GETIDENTITY(_A);
433 int x = A.first_variable();
436 poly lead(A.integer_lcoeff());
437 for (
int i=0; i<(int)a.size(); i++) {
438 a[i] *= lead / a[i].integer_lcoeff();
443 vector<poly> s(solve_Diophantine_univariate(a,
poly(BHEAD 1,A.modp,1)));
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));
452 for (
int k=2; k<=A.modn; k++) {
453 term *=
poly(BHEAD A.modp);
455 poly error(BHEAD -1);
456 for (
int i=0; i<(int)a.size(); i++) error *= a[i];
458 if (error.is_zero())
break;
461 error.setmod(A.modp,1);
463 for (
int i=0; i<(int)a.size(); i++)
464 a[i] += term * (error * s[i] % a[i]);
468 for (
int i=0; i<(int)a.size(); i++)
469 a[i] /= polygcd::integer_content(
poly(a[i],0,1));
472 cout <<
"*** [" << thetime() <<
"] RES : lift_coefficients("<<_A<<
","<<_a<<
") = "<<a<<endl;
497 void polyfact::predetermine (
int dep,
const vector<vector<int> > &state, vector<vector<vector<int> > > &terms, vector<int> &term,
int sumdeg) {
499 if (dep == (
int)state.size()) {
500 terms[sumdeg].push_back(term);
507 for (
int deg=0; sumdeg+deg<(int)state[dep].size(); deg++)
508 if (state[dep][deg] > 0) {
510 predetermine(dep+1, state, terms, term, sumdeg+deg);
558 const vector<poly> polyfact::lift_variables (
const poly &A,
const vector<poly> &_a,
const vector<int> &x,
const vector<int> &c,
const vector<poly> &lc) {
561 cout <<
"*** [" << thetime() <<
"] CALL: lift_variables("<<A<<
","<<_a<<
","<<x<<
","<<c<<
","<<lc<<
")\n";
565 if (x.size()<=1)
return _a;
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();
583 int D = A.degree(x[0]);
584 vector<vector<int> > state(a.size(), vector<int>(D+1, 0));
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;
591 vector<vector<vector<int> > > terms(D+1);
593 predetermine(0,state,terms,term);
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]++;
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]));
610 for (
int i=0; i<(int)terms.size(); i++) {
613 if (num_undet[i] == 1) {
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);
620 for (
int k=0; k<(int)terms[i][j].size(); k++) {
621 if (state[k][terms[i][j][k]] == 1) {
624 which_deg=terms[i][j][k];
627 coeff *= a[k].coefficient(x[0], terms[i][j][k]);
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;
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)
653 poly check(BHEAD 1, A.modn>1?0:A.modp, 1);
654 for (
int i=0; i<(int)a.size(); i++)
657 if (check == A)
return a;
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);
669 for (
int i=1; i<(int)x.size(); i++)
670 maxdegA = MaX(maxdegA, A.degree(x[i]));
673 for (
int xi=1; xi<(int)x.size(); xi++) {
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]));
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];
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);
688 for (
int deg=1, maxdeg=A.degree(x[xi]); deg<=maxdeg; deg++) {
690 term *= simple[xi-1];
693 poly error(BHEAD -1,A.modp,A.modn);
694 for (
int i=0; i<(int)a.size(); i++) error *= a[i];
696 for (
int i=xi; i<(int)c.size(); i++) error %= simple[i];
698 if (error.is_zero())
break;
701 error %= simple[xi-1];
703 vector<poly> s(solve_Diophantine_multivariate(anew,error,xnew,cnew,maxdegA));
704 if (s == vector<poly>())
return vector<poly>();
706 for (
int i=0; i<(int)a.size(); i++)
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];
714 for (
int i=xi; i<(int)c.size(); i++) check %= simple[i];
715 if (!check.is_zero())
return vector<poly>();
719 cout <<
"*** [" << thetime() <<
"] RES : lift_variables("<<A<<
","<<_a<<
","<<x<<
","<<c<<
","<<lc<<
") = " << a << endl;
741 WORD polyfact::choose_prime (
const poly &a,
const vector<int> &x, WORD p) {
744 cout <<
"*** [" << thetime() <<
"] CALL: choose_prime("<<a<<
","<<x<<
","<<p<<
")"<<endl;
749 poly icont_lcoeff(polygcd::integer_content(a.lcoeff_univar(x[0])));
751 if (p==0) p = POLYFACT_FIRST_PRIME;
753 poly icont_lcoeff_modp(BHEAD 0);
762 for (
int d=2; d*d<=p; d++)
763 if (p%d==0) { is_prime=
false;
break; }
767 icont_lcoeff_modp = icont_lcoeff;
768 icont_lcoeff_modp.setmod(p,1);
770 while (icont_lcoeff_modp.is_zero());
773 cout <<
"*** [" << thetime() <<
"] RES : choose_prime("<<a<<
","<<x<<
",?) = "<<p<<endl;
798 WORD polyfact::choose_prime_power (
const poly &a, WORD p) {
801 cout <<
"*** [" << thetime() <<
"] CALL: choose_prime_power("<<a<<
","<<p<<
")"<<endl;
807 double maxdegree=0, maxlogcoeff=0, numterms=0;
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]);
813 maxlogcoeff = MaX(maxlogcoeff,
814 log(1.0+(UWORD)a[i+a[i]-2]) +
815 BITSINWORD*log(2.0)*(ABS(a[i+a[i]-1])-1));
819 WORD res = (WORD)ceil((log((sqrt(5.0)+1)/2)*maxdegree + maxlogcoeff + 0.5*log(numterms)) / log((
double)p));
822 cout <<
"*** [" << thetime() <<
"] CALL: choose_prime_power("<<a<<
","<<p<<
") = "<<res<<endl;
860 const vector<int> polyfact::choose_ideal (
const poly &a,
int p,
const factorized_poly &lc,
const vector<int> &x) {
863 cout <<
"*** [" << thetime() <<
"] CALL: polyfact::choose_ideal(" 864 <<a<<
","<<p<<
","<<lc<<
","<<x<<
")"<<endl;
867 if (x.size()==1)
return vector<int>();
871 vector<int> c(x.size()-1);
873 int dega = a.degree(x[0]);
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);
886 if (amodIp.degree(x[0]) != dega)
887 return c = vector<int>();
890 if (!polygcd::gcd_Euclidean(amodIp, amodIp.derivative(x[0])).is_integer())
891 return c = vector<int>();
893 if (a.modp>0 && a.modn==1)
return c;
896 vector<poly> d(1, polygcd::integer_content(amodI));
898 for (
int i=0; i<(int)lc.factor.size(); i++) {
900 if (i==0 && lc.factor[i].is_integer()) {
901 d[0] *= lc.factor[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);
912 for (
int j=(
int)d.size()-1; j>=0; j--) {
914 while (!r.is_one()) {
915 r = polygcd::integer_gcd(r,q);
921 if (q.is_one())
return vector<int>();
926 cout <<
"*** [" << thetime() <<
"] RES : polyfact::choose_ideal(" 927 <<a<<
","<<p<<
","<<lc<<
","<<x<<
") = "<<c<<endl;
950 int x = a.first_variable();
952 poly b(a.derivative(x));
953 poly c(polygcd::gcd(a,b));
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);
965 if (!a.is_one()) res.add_factor(a,pow);
986 int x = a.first_variable();
987 poly b(a.derivative(x));
991 poly c(polygcd::gcd(a,b));
994 while (!a.is_one()) {
995 b = polygcd::gcd(a,c);
997 if (!a.is_one()) res.add_factor(a,pow);
1008 for (
int i=1; i<a[1]; i+=a[i])
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]);
1048 cout <<
"*** [" << thetime() <<
"] CALL: squarefree_factors("<<a<<
")\n";
1056 res = squarefree_factors_Yun(a);
1058 res = squarefree_factors_modp(a);
1061 cout <<
"*** [" << thetime() <<
"] RES : squarefree_factors("<<a<<
") = "<<res<<
"\n";
1094 const vector<vector<WORD> > polyfact::Berlekamp_Qmatrix (
const poly &_a) {
1097 cout <<
"*** [" << thetime() <<
"] CALL: Berlekamp_Qmatrix("<<_a<<
")\n";
1100 if (_a.all_variables() == vector<int>())
1101 return vector<vector<WORD> >(0);
1103 POLY_GETIDENTITY(_a);
1106 int x = a.first_variable();
1107 int n = a.degree(x);
1110 poly lc(a.integer_lcoeff());
1113 vector<vector<WORD> > Q(n, vector<WORD>(n));
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];
1121 vector<WORD> d(n,0);
1124 for (
int i=0; i<=(n-1)*p; i++) {
1126 if (i%p==0) Q[i/p] = d;
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;
1139 for (
int i=0; i<n; i++)
1140 Q[i][i] = (Q[i][i] - 1 + p) % p;
1143 for (
int i=0; i<n; i++) {
1145 int ii=i;
while (ii<n && Q[i][ii]==0) ii++;
1146 if (ii==n)
continue;
1148 for (
int k=0; k<n; k++)
1149 swap(Q[k][ii],Q[k][i]);
1156 for (
int k=0; k<n; k++)
if (Q[k][i] != 0) {
1159 Q[k][i] = ((LONG)Q[k][i] * mul) % p;
1163 for (
int j=0; j<n; j++)
1164 if (j!=i && Q[i][j]!=0) {
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;
1173 for (
int i=0; i<n; i++) {
1175 Q[i][i] = Q[i][i]-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;
1185 cout <<
"*** [" << thetime() <<
"] RES : Berlekamp_Qmatrix("<<_a<<
") = "<<Q<<
"\n";
1219 const vector<poly> polyfact::Berlekamp_find_factors (
const poly &_a,
const vector<vector<WORD> > &_q) {
1222 cout <<
"*** [" << thetime() <<
"] CALL: Berlekamp_find_factors("<<_a<<
","<<_q<<
")\n";
1225 if (_a.all_variables() == vector<int>())
return vector<poly>(1,_a);
1227 POLY_GETIDENTITY(_a);
1229 vector<vector<WORD> > q=_q;
1232 for (
int i=0; i<(int)q.size(); i++)
1233 if (q[i]!=vector<WORD>(q[i].size(),0)) rank++;
1236 int x = a.first_variable();
1237 int n = a.degree(x);
1240 a /= a.integer_lcoeff();
1243 vector<vector<WORD> > fac(1, vector<WORD>(n+1,0));
1245 fac[0] = poly::to_coefficient_list(a);
1246 bool finished=
false;
1249 for (
int i=1; i<n && !finished; i++) {
1250 if (q[i] == vector<WORD>(n,0))
continue;
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);
1257 if (c.size()!=1 && c.size()!=fac[j].size()) {
1259 fac[j] = poly::coefficient_list_divmod(fac[j], c, p, 0);
1260 if ((
int)fac.size() == rank) finished=
true;
1265 q[i][0] = (q[i][0]+1) % p;
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);
1275 cout <<
"*** [" << thetime() <<
"] RES : Berlekamp_find_factors("<<_a<<
","<<_q<<
") = "<<res<<
"\n";
1307 const vector<poly> polyfact::combine_factors (
const poly &a,
const vector<poly> &f) {
1310 cout <<
"*** [" << thetime() <<
"] CALL: combine_factors("<<a<<
","<<f<<
")\n";
1313 POLY_GETIDENTITY(a);
1319 vector<bool> used(f.size(),
false);
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;
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));
1335 if ((a0 % fac).is_zero()) {
1337 for (
int i=0, j=0; i<(int)f.size(); i++)
1338 if (!used[i]) used[i] = next[j++];
1344 while (next_permutation(next.begin(), next.end()));
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));
1359 cout <<
"*** [" << thetime() <<
"] RES : combine_factors("<<a<<
","<<f<<
") = "<<res<<
"\n";
1388 const vector<poly> polyfact::factorize_squarefree (
const poly &a,
const vector<int> &x) {
1391 cout <<
"*** [" << thetime() <<
"] CALL: factorize_squarefree("<<a<<
")\n";
1394 POLY_GETIDENTITY(a);
1396 WORD p=a.modp, n=a.modn;
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;
1410 int prime_tries = 0;
1412 while (prime_tries<POLYFACT_NUM_CONFIRMATIONS && min_factors>1) {
1414 p = choose_prime(a,x,p);
1416 if (a.degree(x[0]) % p == 0)
continue;
1422 if (polygcd::gcd_Euclidean(amodp, amodp.derivative(x[0])).degree(x[0]) != 0)
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;
1434 if (x.size()==1 || c.size()>0) {
1436 for (
int i=0; i<(int)c.size(); i++)
1437 amodI %= poly::simple_poly(BHEAD x[i+1],c[i]);
1440 q = Berlekamp_Qmatrix(
poly(amodI,p,1));
1442 for (
int i=0; i<(int)q.size(); i++)
1443 if (q[i]!=vector<WORD>(q[i].size(),0)) rank++;
1445 if (rank<min_factors) {
1454 if (rank==min_factors)
1458 cout <<
"*** [" << thetime() <<
"] (A) : factorize_squarefree("<<a<<
1459 ") try p=" << p <<
" #factors=" << rank <<
" (min="<<min_factors<<
"x"<<prime_tries<<
")" << endl;
1471 n = choose_prime_power(amodI,p);
1472 n = MaX(n, choose_prime_power(a,p));
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;
1485 vector<poly> f(Berlekamp_find_factors(
poly(amodI,p,1),q));
1489 f = lift_coefficients(amodI,f);
1490 if (f==vector<poly>()) {
1492 cout <<
"factor_squarefree failed (lift_coeff step) : " << endl;
1498 if (a.modp==0) f = combine_factors(amodI,f);
1503 f = vector<poly>(1, a);
1505 if (x.size() > 1 && f.size() > 1) {
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));
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]);
1525 vector<poly> correct_lc(f.size(),
poly(BHEAD 1,p,n));
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);
1533 for (
int i=(
int)lcmodI.size()-1; i>=0; i--) {
1535 if (i==0 && lc.factor[i].is_integer())
continue;
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);
1544 correct_lc[j] *= lc.factor[i];
1545 lc_f.termscopy(&quo[0], 2+AN.poly_num_vars, ABS(nquo));
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]);
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;
1565 if (!delta.is_one()) {
1566 poly deltapow(BHEAD 1);
1567 for (
int i=1; i<(int)correct_lc.size(); i++)
1569 while (!deltapow.is_zero()) {
1570 deltapow /=
poly(BHEAD p);
1574 for (
int i=0; i<(int)f.size(); i++) {
1576 correct_lc[i].modn = n;
1582 for (
int i=0; i<(int)correct_lc.size(); i++) {
1583 correct_lc[i] *= delta;
1585 if (i>0) aa *= delta;
1588 f = lift_variables(aa,f,x,c,correct_lc);
1590 for (
int i=0; i<(int)f.size(); i++)
1592 f[i] /= polygcd::integer_content(
poly(f[i],0,1));
1594 f[i] /= polygcd::content_univar(f[i], x[0]);
1596 if (f==vector<poly>()) {
1598 cout <<
"factor_squarefree failed (lift_var step)" << endl;
1606 for (
int i=0; i<(int)f.size(); i++)
1610 poly check(BHEAD 1,a.modp,a.modn);
1611 for (
int i=0; i<(int)f.size(); i++)
1616 cout <<
"factor_squarefree failed (final check) : " << f << endl;
1622 cout <<
"*** [" << thetime() <<
"] RES : factorize_squarefree("<<a<<
","<<x<<
","<<c<<
") = "<<f<<
"\n";
1644 cout <<
"*** [" << thetime() <<
"] CALL: factorize("<<a<<
")\n";
1646 vector<int> x = a.all_variables();
1649 if (x.size() == 0) {
1651 if (a.is_one())
return res;
1652 res.add_factor(a,1);
1657 poly conta(polygcd::content_univar(a,x[0]));
1661 poly ppa(a / conta);
1667 cout <<
"*** [" << thetime() <<
"] ... : factorize("<<a<<
") : SFF = "<<b<<
"\n";
1671 for (
int i=0; i<(int)b.factor.size(); i++) {
1673 vector<poly> c(factorize_squarefree(b.factor[i], x));
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]);
1682 cout <<
"*** [" << thetime() <<
"] RES : factorize("<<a<<
") = "<<faca<<
"\n";
int GetModInverses(WORD, WORD, WORD *, WORD *)