61template<
class T> ostream& operator<< (ostream &out,
const vector<T> &x) {
63 for (
int i=0; i<(int)x.size(); i++) {
91const poly polygcd::integer_gcd (
const poly &a,
const poly &b) {
94 cout <<
"*** [" << thetime() <<
"] CALL: integer_gcd(" << a <<
"," << b <<
")" << endl;
99 if (a.is_zero())
return b;
100 if (b.is_zero())
return a;
102 poly c(BHEAD 0, a.modp, a.modn);
106 (UWORD *)&a[AN.poly_num_vars+2],a[a[0]-1],
107 (UWORD *)&b[AN.poly_num_vars+2],b[b[0]-1],
108 (UWORD *)&c[AN.poly_num_vars+2],&nc);
110 WORD x = 2 + AN.poly_num_vars + ABS(nc);
115 for (
int i=0; i<AN.poly_num_vars; i++)
119 cout <<
"*** [" << thetime() <<
"] RES : integer_gcd(" << a <<
"," << b <<
") = " << c << endl;
143const poly polygcd::integer_content (
const poly &a) {
146 cout <<
"*** [" << thetime() <<
"] CALL: integer_content(" << a <<
")" << endl;
151 if (a.modp>0)
return a.integer_lcoeff();
153 poly c(BHEAD 0, 0, 1);
154 WORD *d = (WORD *)NumberMalloc(
"polygcd::integer_content");
157 for (
int i=0; i<AN.poly_num_vars; i++)
160 for (
int i=1; i<a[0]; i+=a[i]) {
162 WCOPY(d,&c[2+AN.poly_num_vars],nc);
164 GcdLong(BHEAD (UWORD *)d, nc,
165 (UWORD *)&a[i+1+AN.poly_num_vars], a[i+a[i]-1],
166 (UWORD *)&c[2+AN.poly_num_vars], &nc);
168 WORD x = 2 + AN.poly_num_vars + ABS(nc);
174 if (a.sign() != c.sign()) c *=
poly(BHEAD -1);
176 NumberFree(d,
"polygcd::integer_content");
179 cout <<
"*** [" << thetime() <<
"] RES : integer_content(" << a <<
") = " << c << endl;
205const poly polygcd::content_univar (
const poly &a,
int x) {
208 cout <<
"*** [" << thetime() <<
"] CALL: content_univar(" << a <<
"," << string(1,
'a'+x) <<
")" << endl;
213 poly res(BHEAD 0, a.modp, a.modn);
215 for (
int i=1; i<a[0];) {
216 poly b(BHEAD 0, a.modp, a.modn);
219 for (; i<a[0] && a[i+1+x]==deg; i+=a[i]) {
220 b.check_memory(b[0]+a[i]);
221 b.termscopy(&a[i],b[0],a[i]);
228 if (res.is_integer()) {
229 res = integer_content(a);
234 if (a.sign() != res.sign()) res *=
poly(BHEAD -1);
237 cout <<
"*** [" << thetime() <<
"] RES : content_univar(" << a <<
"," << string(1,
'a'+x) <<
") = " << res << endl;
263const poly polygcd::content_multivar (
const poly &a,
int x) {
266 cout <<
"*** [" << thetime() <<
"] CALL: content_multivar(" << a <<
"," << string(1,
'a'+x) <<
")" << endl;
271 poly res(BHEAD 0, a.modp, a.modn);
273 for (
int i=1,j; i<a[0]; i=j) {
274 poly b(BHEAD 0, a.modp, a.modn);
276 for (j=i; j<a[0]; j+=a[j]) {
277 bool same_powers =
true;
278 for (
int k=0; k<AN.poly_num_vars; k++)
279 if (k!=x && a[i+1+k]!=a[j+1+k]) {
283 if (!same_powers)
break;
285 b.check_memory(b[0]+a[j]);
286 b.termscopy(&a[j],b[0],a[j]);
287 for (
int k=0; k<AN.poly_num_vars; k++)
288 if (k!=x) b[b[0]+1+k]=0;
293 res = gcd_Euclidean(res, b);
295 if (res.is_integer()) {
296 res =
poly(BHEAD a.sign(),a.modp,a.modn);
302 cout <<
"*** [" << thetime() <<
"] RES : content_multivar(" << a <<
"," << string(1,
'a'+x) <<
") = " << res << endl;
325const vector<WORD> polygcd::coefficient_list_gcd (
const vector<WORD> &_a,
const vector<WORD> &_b, WORD p) {
328 cout <<
"*** [" << thetime() <<
"] CALL: coefficient_list_gcd("<<_a<<
","<<_b<<
","<<p<<
")"<<endl;
331 vector<WORD> a(_a), b(_b);
333 while (b.size() != 0) {
334 a = poly::coefficient_list_divmod(a,b,p,1);
338 while (a.back()==0) a.pop_back();
343 for (
int i=0; i<(int)a.size(); i++)
344 a[i] = (LONG)inv*a[i] % p;
347 cout <<
"*** [" << thetime() <<
"] RES : coefficient_list_gcd("<<_a<<
","<<_b<<
","<<p<<
") = "<<a<<endl;
374const poly polygcd::gcd_Euclidean (
const poly &a,
const poly &b) {
377 cout <<
"*** [" << thetime() <<
"] CALL: gcd_Euclidean("<<a<<
","<<b<<
")"<<endl;
382 if (a.is_zero())
return b;
383 if (b.is_zero())
return a;
384 if (a.is_integer() || b.is_integer())
return integer_gcd(a,b);
389 vector<WORD> coeff = coefficient_list_gcd(poly::to_coefficient_list(a),
390 poly::to_coefficient_list(b), a.modp);
391 res = poly::from_coefficient_list(BHEAD coeff, a.first_variable(), a.modp);
396 while (!rem.is_zero())
398 res /= res.integer_lcoeff();
402 cout <<
"*** [" << thetime() <<
"] RES : gcd_Euclidean("<<a<<
","<<b<<
") = "<<res<<endl;
426const poly polygcd::chinese_remainder (
const poly &a1,
const poly &m1,
const poly &a2,
const poly &m2) {
429 cout <<
"*** [" << thetime() <<
"] CALL: chinese_remainder(" << a1 <<
"," << m1 <<
"," << a2 <<
"," << m2 <<
")" << endl;
432 POLY_GETIDENTITY(a1);
435 UWORD *x = (UWORD *)NumberMalloc(
"polygcd::chinese_remainder");
436 UWORD *y = (UWORD *)NumberMalloc(
"polygcd::chinese_remainder");
437 UWORD *z = (UWORD *)NumberMalloc(
"polygcd::chinese_remainder");
439 GetLongModInverses(BHEAD (UWORD *)&m1[2+AN.poly_num_vars], m1[m1[1]],
440 (UWORD *)&m2[2+AN.poly_num_vars], m2[m2[1]],
441 (UWORD *)x, &nx, NULL, NULL);
443 AddLong((UWORD *)&a2[2+AN.poly_num_vars], a2.is_zero() ? 0 : a2[a2[1]],
444 (UWORD *)&a1[2+AN.poly_num_vars], a1.is_zero() ? 0 : -a1[a1[1]],
447 MulLong (x,nx,y,ny,z,&nz);
448 MulLong (z,nz,(UWORD *)&m1[2+AN.poly_num_vars],m1[m1[1]],x,&nx);
450 AddLong (x,nx,(UWORD *)&a1[2+AN.poly_num_vars], a1.is_zero() ? 0 : a1[a1[1]],y,&ny);
452 MulLong ((UWORD *)&m1[2+AN.poly_num_vars], m1[m1[1]],
453 (UWORD *)&m2[2+AN.poly_num_vars], m2[m2[1]],
456 TakeNormalModulus (y,&ny,(UWORD *)z,nz,NOUNPACK);
458 poly res(BHEAD y,ny);
460 NumberFree(x,
"polygcd::chinese_remainder");
461 NumberFree(y,
"polygcd::chinese_remainder");
462 NumberFree(z,
"polygcd::chinese_remainder");
465 cout <<
"*** [" << thetime() <<
"] RES : chinese_remainder(" << a1 <<
"," << m1 <<
"," << a2 <<
"," << m2 <<
") = " << res << endl;
488const poly polygcd::substitute(
const poly &a,
int x,
int c) {
503 vector<WORD> cache(min(a.degree(x)+1,min(2*a.number_of_terms(),
504 POLYGCD_RAISPOWMOD_CACHE_SIZE)), 0);
506 for (
int ai=1; ai<=a[0]; ai+=a[ai]) {
513 for (
int i=0; i<AN.poly_num_vars; i++)
514 if (i!=x && a[ai+1+i]!=b[bi+1+i]) {
522 if (b[bi+AN.poly_num_vars+1] < 0)
523 b[bi+AN.poly_num_vars+1] += a.modp;
534 b[bi] = 3+AN.poly_num_vars;
535 for (
int i=0; i<AN.poly_num_vars; i++)
536 b[bi+1+i] = a[ai+1+i];
538 b[bi+AN.poly_num_vars+1] = 0;
539 b[bi+AN.poly_num_vars+2] = 1;
543 LONG coeff = a[ai+1+AN.poly_num_vars] * a[ai+2+AN.poly_num_vars];
546 if (pow<(
int)cache.size()) {
548 cache[pow] = RaisPowMod(c, pow, a.modp);
549 coeff = (coeff * cache[pow]) % a.modp;
552 coeff = (coeff * RaisPowMod(c, pow, a.modp)) % a.modp;
555 b[bi+AN.poly_num_vars+1] = (coeff + b[bi+AN.poly_num_vars+1]) % a.modp;
556 if (b[bi+AN.poly_num_vars+1] != 0) zero=
false;
571const vector<int> polygcd::sparse_interpolation_get_mul_list (
const poly &a,
const vector<int> &x,
const vector<int> &c) {
574 vector<vector<WORD> > cache(c.size());
575 int max_cache_size = min(2*a.number_of_terms(),POLYGCD_RAISPOWMOD_CACHE_SIZE);
576 for (
int i=0; i<(int)c.size(); i++)
577 cache[i] = vector<WORD>(min(a.degree(x[i+1])+1,max_cache_size), 0);
580 for (
int i=1; i<a[0]; i+=a[i]) {
582 for (
int j=0; j<(int)c.size(); j++) {
583 int pow = a[i+1+x[j+1]];
584 if (pow<(
int)cache[j].size()) {
585 if (cache[j][pow]==0)
586 cache[j][pow] = RaisPowMod(c[j], pow, a.modp);
587 coeff = (coeff * cache[j][pow]) % a.modp;
590 coeff = (coeff * RaisPowMod(c[j], pow, a.modp)) % a.modp;
593 res.push_back(coeff);
599void polygcd::sparse_interpolation_mul_poly (
poly &a,
const vector<int> &mul) {
600 for (
int i=1,j=0; i<a[0]; i+=a[i],j++)
601 a[i+a[i]-2] = ((LONG)a[i+a[i]-2]*mul[j]) % a.modp;
605const poly polygcd::sparse_interpolation_reduce_poly (
const poly &a,
const vector<int> &x) {
607 for (
int i=1; i<a[0]; i+=a[i]) {
608 for (
int j=1; j<(int)x.size(); j++)
610 if (res[i+a[i]-1]==-1) {
612 res[i+a[i]-2]=a.modp-res[i+a[i]-2];
619const poly polygcd::sparse_interpolation_fix_poly (
const poly &a,
int x) {
622 poly res(BHEAD 0,a.modp,1);
627 for (
int i=1; i<a[0]; i+=a[i]) {
629 res.termscopy(&a[i], j, a[i]);
631 res[j+res[j]-2] = ((LONG)res[j+res[j]-2] + a[i+a[i]-2]) % a.modp;
633 newterm = i+a[i] == a[0] || res[j+1+x] != a[i+a[i]+1+x];
634 if (newterm && res[j+res[j]-2]!=0) j += res[j];
677const poly polygcd::gcd_modular_sparse_interpolation (
const poly &origa,
const poly &origb,
const vector<int> &x,
const poly &s) {
680 cout <<
"*** [" << thetime() <<
"] CALL: gcd_modular_sparse_interpolation("
681 << origa <<
"," << origb <<
"," << x <<
"," <<
"," << s <<
")" << endl;
684 POLY_GETIDENTITY(origa);
687 poly conta(content_multivar(origa,x.back()));
688 poly contb(content_multivar(origb,x.back()));
689 poly gcdconts(gcd_Euclidean(conta,contb));
690 const poly& a = conta.is_one() ? origa : origa/conta;
691 const poly& b = contb.is_one() ? origb : origb/contb;
696 poly lcgcd(BHEAD 1, a.modp);
697 if (!s.lcoeff_univar(x[0]).is_integer()) {
698 lcgcd = gcd_modular_dense_interpolation(a.lcoeff_univar(x[0]), b.lcoeff_univar(x[0]), x,
poly(BHEAD 0));
702 poly ared(sparse_interpolation_reduce_poly(a,x));
703 poly bred(sparse_interpolation_reduce_poly(b,x));
704 poly sred(sparse_interpolation_reduce_poly(s,x));
705 poly lred(sparse_interpolation_reduce_poly(lcgcd,x));
709 for (
int i=1; i<sred[0]; i+=sred[i]) {
710 sred[i+sred[i]-2] = sred[i+sred[i]-1] = 1;
716 vector<int> c(x.size()-1);
721 for (
int i=0; i<(int)c.size(); i++)
722 c[i] = 1 + wranf(BHEAD0) % (a.modp-1);
723 smul = sparse_interpolation_get_mul_list(s,x,c);
728 for (
int i=1; i<s[0];) {
729 int pow = s[i+1+x[0]];
730 while (i<s[0] && s[i+1+x[0]]==pow) i+=s[i], to++;
731 for (
int j=fr; j<to; j++)
732 for (
int k=fr; k<j; k++)
733 if (smul[j] == smul[k])
741 vector<int> amul(sparse_interpolation_get_mul_list(a,x,c));
742 vector<int> bmul(sparse_interpolation_get_mul_list(b,x,c));
743 vector<int> lmul(sparse_interpolation_get_mul_list(lcgcd,x,c));
745 vector<vector<vector<LONG> > > M;
746 vector<vector<LONG> > V;
751 for (
int i=1; i<s[0]; i+=s[i]) {
752 if (i==1 || s[i+1+x[0]]!=s[i+1+x[0]-s[i]]) {
753 M.push_back(vector<vector<LONG> >());
754 V.push_back(vector<LONG>());
756 M.back().push_back(vector<LONG>());
757 V.back().push_back(0);
758 maxMsize = max(maxMsize, (
int)M.back().size());
762 for (
int numg=0; numg<maxMsize; numg++) {
764 poly amodI(sparse_interpolation_fix_poly(ared,x[0]));
765 poly bmodI(sparse_interpolation_fix_poly(bred,x[0]));
766 poly lmodI(sparse_interpolation_fix_poly(lred,x[0]));
771 poly gcd(lmodI * gcd_Euclidean(amodI,bmodI));
774 if (!gcd.is_zero() && gcd[2+x[0]]==sred[2+x[0]]) {
779 for (
int si=1; si<s[0];) {
781 if (gi<gcd[0] && gcd[gi+1+x[0]]==sred[si+1+x[0]]) {
782 if (numg < (
int)V[midx].size())
783 V[midx][numg] = gcd[gi+gcd[gi]-1]*gcd[gi+gcd[gi]-2];
788 for (
int i=0; i<(int)M[midx].size(); i++) {
789 if (numg < (
int)M[midx].size())
790 M[midx][numg].push_back(sred[si+1+AN.poly_num_vars]);
799 if (!gcd.is_zero() && gcd[2+x[0]]<sred[2+x[0]])
800 return poly(BHEAD 0);
805 sparse_interpolation_mul_poly(ared,amul);
806 sparse_interpolation_mul_poly(bred,bmul);
807 sparse_interpolation_mul_poly(sred,smul);
808 sparse_interpolation_mul_poly(lred,lmul);
812 for (
int i=0; i<(int)M.size(); i++) {
816 for (
int j=0; j<n; j++) {
817 for (
int k=0; k<j; k++) {
819 for (
int l=k; l<n; l++)
820 M[i][j][l] = (M[i][j][l] - M[i][k][l]*x) % a.modp;
821 V[i][j] = (V[i][j] - V[i][k]*x) % a.modp;
827 for (
int k=0; k<n; k++)
828 M[i][j][k] = (M[i][j][k]*x) % a.modp;
829 V[i][j] = (V[i][j]*x) % a.modp;
833 for (
int j=n-1; j>=0; j--)
834 for (
int k=j+1; k<n; k++)
835 V[i][j] = (V[i][j] - M[i][j][k]*V[i][k]) % a.modp;
840 for (
int i=0; i<(int)V.size(); i++)
841 for (
int j=0; j<(int)V[i].size(); j++)
842 coeff.push_back(V[i][j]);
847 for (
int si=1; si<s[0]; si+=s[si]) {
848 res.check_memory(ri);
849 res[ri] = 3 + AN.poly_num_vars;
850 for (
int j=0; j<AN.poly_num_vars; j++)
851 res[ri+1+j] = s[si+1+j];
852 res[ri+1+AN.poly_num_vars] = ABS(coeff[i]);
853 res[ri+2+AN.poly_num_vars] = SGN(coeff[i]);
858 res.setmod(a.modp,1);
860 if (!poly::divides(res, lcgcd * a) || !poly::divides(res, lcgcd * b)) {
861 return poly(BHEAD 0);
864 if (!poly::divides(res, a))
865 res = gcd_modular_dense_interpolation(res, a, x,
poly(BHEAD 0));
866 if (!poly::divides(res, b))
867 res = gcd_modular_dense_interpolation(res, b, x,
poly(BHEAD 0));
871 cout <<
"*** [" << thetime() <<
"] RES : gcd_modular_sparse_interpolation("
872 << a <<
"," << b <<
"," << x <<
"," <<
"," << s <<
") = " << res << endl;
875 return gcdconts * res;
901const poly polygcd::gcd_modular_dense_interpolation (
const poly &a,
const poly &b,
const vector<int> &x,
const poly &s) {
904 cout <<
"*** [" << thetime() <<
"] CALL: gcd_modular_dense_interpolation(" << a <<
"," << b <<
"," << x <<
"," << s <<
")" << endl;
911 return gcd_Euclidean(a,b);
916 poly res = gcd_modular_sparse_interpolation (a,b,x,s);
917 if (!res.is_zero())
return res;
924 poly conta(content_multivar(a,X));
925 poly contb(content_multivar(b,X));
926 poly gcdconts(gcd_Euclidean(conta,contb));
927 const poly& ppa = conta.is_one() ? a :
poly(a/conta);
928 const poly& ppb = contb.is_one() ? b :
poly(b/contb);
931 poly lcoeffa(ppa.lcoeff_multivar(X));
932 poly lcoeffb(ppb.lcoeff_multivar(X));
933 poly gcdlcoeffs(gcd_Euclidean(lcoeffa,lcoeffb));
936 int m = MiN(ppa.degree(x[x.size() - 2]),ppb.degree(x[x.size() - 2]));
939 poly oldres(BHEAD 0);
940 poly newshape(BHEAD 0);
941 poly modpoly(BHEAD 1,a.modp);
945 int c = 1 + wranf(BHEAD0) % (a.modp-1);
946 if (substitute(gcdlcoeffs,X,c).is_zero())
continue;
947 if (substitute(modpoly,X,c).is_zero())
continue;
949 poly amodc(substitute(ppa,X,c));
950 poly bmodc(substitute(ppb,X,c));
953 poly gcdmodc(gcd_modular_dense_interpolation(amodc,bmodc,vector<int>(x.begin(),x.end()-1), newshape));
954 int n = gcdmodc.degree(x[x.size() - 2]);
957 gcdmodc = (gcdmodc * substitute(gcdlcoeffs,X,c)) / gcdmodc.integer_lcoeff();
958 poly simple(poly::simple_poly(BHEAD X,c,1,a.modp));
961 if ((res.is_zero() && n == m) || n < m) {
970 poly coeff_poly(substitute(modpoly,X,c));
971 WORD coeff_word = coeff_poly[2+AN.poly_num_vars] * coeff_poly[3+AN.poly_num_vars];
972 if (coeff_word < 0) coeff_word += a.modp;
977 res +=
poly(BHEAD coeff_word, a.modp, 1) * modpoly * (gcdmodc - substitute(res,X,c));
982 if (!res.is_zero() && res == oldres) {
983 poly nres = res / content_multivar(res, X);
984 if (poly::divides(nres,ppa) && poly::divides(nres,ppb)) {
986 cout <<
"*** [" << thetime() <<
"] RES : gcd_modular_dense_interpolation(" << a <<
"," << b <<
","
987 << x <<
"," <<
"," << s <<
") = " << gcdconts * nres << endl;
989 return gcdconts * nres;
995 newshape =
poly(BHEAD 0);
1021const poly polygcd::gcd_modular (
const poly &origa,
const poly &origb,
const vector<int> &x) {
1024 cout <<
"*** [" << thetime() <<
"] CALL: gcd_modular(" << origa <<
"," << origb <<
"," << x <<
")" << endl;
1027 POLY_GETIDENTITY(origa);
1029 if (origa.is_zero())
return origa.modp==0 ? origb : origb / origb.integer_lcoeff();
1030 if (origb.is_zero())
return origa.modp==0 ? origa : origa / origa.integer_lcoeff();
1031 if (origa==origb)
return origa.modp==0 ? origa : origa / origa.integer_lcoeff();
1033 poly ac = integer_content(origa);
1034 poly bc = integer_content(origb);
1035 const poly& a = ac.is_one() ? origa :
poly(origa/ac);
1036 const poly& b = bc.is_one() ? origb :
poly(origb/bc);
1037 poly ic = integer_gcd(ac, bc);
1038 poly g = integer_gcd(a.integer_lcoeff(), b.integer_lcoeff());
1044 int mindeg=MAXPOSITIVE;
1049 if (
poly(a.integer_lcoeff(),p).is_zero())
continue;
1050 if (
poly(b.integer_lcoeff(),p).is_zero())
continue;
1053 c = (c *
poly(g,p)) / c.integer_lcoeff();
1059 mindeg = MAXPOSITIVE;
1063 if (!(
poly(a,p)%c).is_zero())
continue;
1064 if (!(
poly(b,p)%c).is_zero())
continue;
1066 int deg = c.degree(x[0]);
1076 else if (deg == mindeg) {
1080 for (
int ci=1,di=1; ci<c[0]||di<d[0]; ) {
1081 int comp = ci==c[0] ? -1 : di==d[0] ? +1 : poly::monomial_compare(BHEAD &c[ci],&d[di]);
1082 poly a1(BHEAD 0),a2(BHEAD 0);
1084 newd.check_memory(newd[0]);
1087 newd.termscopy(&d[di],newd[0],1+AN.poly_num_vars);
1088 a1 =
poly(BHEAD (UWORD *)&d[di+1+AN.poly_num_vars],d[di+d[di]-1]);
1092 newd.termscopy(&c[ci],newd[0],1+AN.poly_num_vars);
1093 a2 =
poly(BHEAD (UWORD *)&c[ci+1+AN.poly_num_vars],c[ci+c[ci]-1]);
1097 poly e(chinese_remainder(a1,m1,a2,
poly(BHEAD p)));
1098 newd.termscopy(&e[2+AN.poly_num_vars], newd[0]+1+AN.poly_num_vars, ABS(e[e[1]])+1);
1099 newd[newd[0]] = 2 + AN.poly_num_vars + ABS(e[e[1]]);
1100 newd[0] += newd[newd[0]];
1103 m1 *=
poly(BHEAD p);
1108 poly ppd(d / integer_content(d));
1111 if (poly::divides(ppd,a) && poly::divides(ppd,b)) {
1113 cout <<
"*** [" << thetime() <<
"] RES : gcd_modular(" << origa <<
"," << origb <<
"," << x <<
") = "
1114 << ic * ppd << endl;
1119 cout <<
"*** [" << thetime() <<
"] Retrying modular_gcd with new prime" << endl;
1146 POLY_GETIDENTITY(a);
1148 double prod_deg = 1;
1149 for (
int j=0; j<AN.poly_num_vars; j++)
1150 prod_deg *= a[2+j]+1;
1152 double digits = ABS(a[1+a[1]-1]);
1153 double lead = a[1+1+AN.poly_num_vars];
1155 return prod_deg*(digits-1+log(2*ABS(lead))/log(2.0)/(BITSINWORD/2)) < POLYGCD_HEURISTIC_MAX_DIGITS;
1190const poly polygcd::gcd_heuristic (
const poly &a,
const poly &b,
const vector<int> &x,
int max_tries) {
1193 cout <<
"*** [" << thetime() <<
"] CALL: gcd_heuristic("<<a<<
","<<b<<
","<<x<<
")\n";
1196 if (a.is_integer())
return integer_gcd(a,integer_content(b));
1197 if (b.is_integer())
return integer_gcd(integer_content(a),b);
1199 POLY_GETIDENTITY(a);
1206 for (
int i=1; i<a[0]; i+=a[i]) {
1207 WORD na = ABS(a[i+a[i]-1]);
1208 if (BigLong((UWORD *)&a[i+a[i]-1-na], na, pxi, nxi) > 0) {
1209 pxi = (UWORD *)&a[i+a[i]-1-na];
1214 for (
int i=1; i<b[0]; i+=b[i]) {
1215 WORD nb = ABS(b[i+b[i]-1]);
1216 if (BigLong((UWORD *)&b[i+b[i]-1-nb], nb, pxi, nxi) > 0) {
1217 pxi = (UWORD *)&b[i+b[i]-1-nb];
1222 poly xi(BHEAD pxi,nxi);
1225 xi = xi*
poly(BHEAD 2) +
poly(BHEAD 2 + wranf(BHEAD0)%POLYGCD_HEURISTIC_MAX_ADD_RANDOM);
1228 if (max(a.degree(x[0]),b.degree(x[0])) * xi[xi[1]] >= MiN(AM.MaxTal, POLYGCD_HEURISTIC_MAX_DIGITS)) {
1230 cout <<
"*** [" << thetime() <<
"] RES : gcd_heuristic("<<a<<
","<<b<<
","<<x<<
") = overflow\n";
1235 for (
int times=0; times<max_tries; times++) {
1239 poly gamma(gcd_heuristic(a % poly::simple_poly(BHEAD x[0],xi),
1240 b % poly::simple_poly(BHEAD x[0],xi),
1241 vector<int>(x.begin()+1,x.end()),1));
1244 if (!gamma.is_zero()) {
1247 poly res(BHEAD 0), c(BHEAD 0);
1248 vector<int> idx, len;
1250 for (
int power=0; !gamma.is_zero(); power++) {
1254 c.coefficients_modulo((UWORD *)&xi[2+AN.poly_num_vars], xi[xi[0]-1],
false);
1257 res.check_memory(res[0]+c[0]);
1258 res.termscopy(&c[1],res[0],c[0]-1);
1259 for (
int i=1; i<c[0]; i+=c[i])
1260 res[res[0]-1+i+1+x[0]] = power;
1264 idx.push_back(res[0]);
1265 len.push_back(c[0]-1);
1271 gamma = (gamma - c) / xi;
1276 rev.check_memory(res[0]);
1279 for (
int i=idx.size()-1; i>=0; i--) {
1280 rev.termscopy(&res[idx[i]], rev[0], len[i]);
1285 poly ppres = res / integer_content(res);
1287 if ((a%ppres).is_zero() && (b%ppres).is_zero()) {
1289 cout <<
"*** [" << thetime() <<
"] RES : gcd_heuristic("<<a<<
","<<b<<
","<<x<<
") = "<<res<<
"\n";
1296 xi = xi *
poly(BHEAD 28657) /
poly(BHEAD 17711) +
poly(BHEAD wranf(BHEAD0) % POLYGCD_HEURISTIC_MAX_ADD_RANDOM);
1300 cout <<
"*** [" << thetime() <<
"] RES : gcd_heuristic("<<a<<
","<<b<<
","<<x<<
") = failed\n";
1303 return poly(BHEAD 0);
1311const map<vector<int>,
poly> polygcd::full_bracket(
const poly &a,
const vector<int>& filter) {
1312 POLY_GETIDENTITY(a);
1314 map<vector<int>,
poly> bracket;
1315 for (
int ai=1; ai<a[0]; ai+=a[ai]) {
1316 vector<int> varpattern(AN.poly_num_vars);
1317 for (
int i=0; i<AN.poly_num_vars; i++)
1318 if (filter[i] == 1 && a[ai + i + 1] > 0)
1319 varpattern[i] = a[ai + i + 1];
1325 for (
int i=0; i<a[ai]; i++)
1326 if (i > 0 && i <= AN.poly_num_vars && varpattern[i - 1])
1331 map<vector<int>,
poly>::iterator i = bracket.find(varpattern);
1332 if (i == bracket.end()) {
1333 bracket.insert(std::make_pair(varpattern, mon));
1342const poly polygcd::bracket(
const poly &a,
const vector<int>& pattern,
const vector<int>& filter) {
1343 POLY_GETIDENTITY(a);
1345 poly bracket(BHEAD 0);
1346 for (
int ai=1; ai<a[0]; ai+=a[ai]) {
1348 for (
int i=0; i<AN.poly_num_vars; i++)
1349 if (filter[i] == 1 && pattern[i] != a[ai + i + 1]) {
1358 for (
int i=0; i<a[ai]; i++)
1359 if (i > 0 && i <= AN.poly_num_vars && pattern[i - 1])
1370const map<vector<int>,
int> polygcd::bracket_count(
const poly &a,
const vector<int>& filter) {
1371 POLY_GETIDENTITY(a);
1373 map<vector<int>,
int> bracket;
1374 for (
int ai=1; ai<a[0]; ai+=a[ai]) {
1375 vector<int> varpattern(AN.poly_num_vars);
1376 for (
int i=0; i<AN.poly_num_vars; i++)
1377 if (filter[i] == 1 && a[ai + i + 1] > 0)
1378 varpattern[i] = a[ai + i + 1];
1380 map<vector<int>,
int>::iterator i = bracket.find(varpattern);
1381 if (i == bracket.end()) {
1382 bracket.insert(std::make_pair(varpattern, 0));
1392 std::vector<int> pattern;
1393 int num_terms, dummy;
1396 BracketInfo(
const std::vector<int>& pattern,
int num_terms,
const poly* p) : pattern(pattern), num_terms(num_terms), p(p) {}
1397 bool operator<(
const BracketInfo& rhs)
const {
return num_terms > rhs.num_terms; }
1405const poly gcd_linear_helper (
const poly &a,
const poly &b) {
1406 POLY_GETIDENTITY(a);
1408 for (
int i = 0; i < AN.poly_num_vars; i++)
1409 if (a.degree(i) == 1) {
1410 vector<int> filter(AN.poly_num_vars);
1414 map<vector<int>,
poly> ba = polygcd::full_bracket(a, filter);
1416 poly subgcd(BHEAD 1);
1418 subgcd = polygcd::gcd_linear(ba.begin()->second, (++ba.begin())->second);
1420 subgcd = ba.begin()->second;
1422 poly linfac = a / subgcd;
1423 if (poly::divides(linfac,b))
1424 return linfac * polygcd::gcd_linear(subgcd, b / linfac);
1426 return polygcd::gcd_linear(subgcd, b);
1429 return poly(BHEAD 0);
1436const poly polygcd::gcd_linear (
const poly &a,
const poly &b) {
1438 cout <<
"*** [" << thetime() <<
"] CALL: gcd_linear("<<a<<
","<<b<<
")\n";
1441 POLY_GETIDENTITY(a);
1443 if (a.is_zero())
return a.modp==0 ? b : b / b.integer_lcoeff();
1444 if (b.is_zero())
return a.modp==0 ? a : a / a.integer_lcoeff();
1445 if (a==b)
return a.modp==0 ? a : a / a.integer_lcoeff();
1447 if (a.is_integer() || b.is_integer()) {
1448 if (a.modp > 0)
return poly(BHEAD 1,a.modp,a.modn);
1449 return poly(integer_gcd(integer_content(a),integer_content(b)),0,1);
1452 poly h = gcd_linear_helper(a, b);
1453 if (!h.is_zero())
return h;
1455 h = gcd_linear_helper(b, a);
1456 if (!h.is_zero())
return h;
1458 vector<int> xa = a.all_variables();
1459 vector<int> xb = b.all_variables();
1461 vector<int> used(AN.poly_num_vars,0);
1462 for (
int i=0; i<(int)xa.size(); i++) used[xa[i]]++;
1463 for (
int i=0; i<(int)xb.size(); i++) used[xb[i]]++;
1465 for (
int i=0; i<AN.poly_num_vars; i++)
1466 if (used[i]) x.push_back(i);
1468 return gcd_modular(a,b,x);
1496const poly polygcd::gcd (
const poly &a,
const poly &b) {
1499 cout <<
"*** [" << thetime() <<
"] CALL: gcd("<<a<<
","<<b<<
")\n";
1502 POLY_GETIDENTITY(a);
1504 if (a.is_zero())
return a.modp==0 ? b : b / b.integer_lcoeff();
1505 if (b.is_zero())
return a.modp==0 ? a : a / a.integer_lcoeff();
1506 if (a==b)
return a.modp==0 ? a : a / a.integer_lcoeff();
1508 if (a.is_integer() || b.is_integer()) {
1509 if (a.modp > 0)
return poly(BHEAD 1,a.modp,a.modn);
1510 return poly(integer_gcd(integer_content(a),integer_content(b)),0,1);
1514 vector<int> xa = a.all_variables();
1515 vector<int> xb = b.all_variables();
1517 vector<int> used(AN.poly_num_vars,0);
1518 for (
int i=0; i<(int)xa.size(); i++) used[xa[i]]++;
1519 for (
int i=0; i<(int)xb.size(); i++) used[xb[i]]++;
1521 for (
int i=0; i<AN.poly_num_vars; i++)
1522 if (used[i]) x.push_back(i);
1524 if (a.is_monomial() || b.is_monomial()) {
1526 poly res(BHEAD 1,a.modp,a.modn);
1527 if (a.modp == 0) res = integer_gcd(integer_content(a),integer_content(b));
1529 for (
int i=0; i<(int)x.size(); i++)
1530 res[2+x[i]] = 1<<(BITSINWORD-2);
1532 for (
int i=1; i<a[0]; i+=a[i])
1533 for (
int j=0; j<(int)x.size(); j++)
1534 res[2+x[j]] = MiN(res[2+x[j]], a[i+1+x[j]]);
1536 for (
int i=1; i<b[0]; i+=b[i])
1537 for (
int j=0; j<(int)x.size(); j++)
1538 res[2+x[j]] = MiN(res[2+x[j]], b[i+1+x[j]]);
1544 poly conta(x.size()==1 ? integer_content(a) : content_univar(a,x[0]));
1545 poly contb(x.size()==1 ? integer_content(b) : content_univar(b,x[0]));
1546 poly gcdconts(x.size()==1 ? integer_gcd(conta,contb) : gcd(conta,contb));
1547 const poly& ppa = conta.is_one() ? a :
poly(a/conta);
1548 const poly& ppb = contb.is_one() ? b :
poly(b/contb);
1551 return ppa * gcdconts;
1555#ifdef POLYGCD_USE_HEURISTIC
1559 res = gcd_heuristic(ppa,ppb,x);
1560 if (!res.is_zero()) res /= integer_content(res);
1568 if (res.is_zero()) {
1569 bool unusedVars =
false;
1570 for (
unsigned int i = 0; i < used.size(); i++) {
1579 res = gcd_linear(ppa,ppb);
1581 cout <<
"New GCD attempt (unused vars): " << res << endl;
1588 bool diva = !res.is_zero() && poly::divides(res,ppa);
1589 bool divb = !res.is_zero() && poly::divides(res,ppb);
1590 if (!diva || !divb) {
1591 vector<BracketInfo> bracketinfo;
1594 map<vector<int>,
int> ba = bracket_count(ppa, used);
1595 for(map<vector<int>,
int>::iterator it = ba.begin(); it != ba.end(); it++)
1596 bracketinfo.push_back(
BracketInfo(it->first, it->second, &ppa));
1600 map<vector<int>,
int> bb = bracket_count(ppb, used);
1601 for(map<vector<int>,
int>::iterator it = bb.begin(); it != bb.end(); it++)
1602 bracketinfo.push_back(
BracketInfo(it->first, it->second, &ppb));
1606 sort(bracketinfo.begin(), bracketinfo.end());
1608 if (res.is_zero()) {
1609 res = bracket(*bracketinfo.back().p, bracketinfo.back().pattern, used);
1610 bracketinfo.pop_back();
1613 while (bracketinfo.size() > 0) {
1614 poly subpoly(bracket(*bracketinfo.back().p, bracketinfo.back().pattern, used));
1615 if (!poly::divides(res,subpoly)) {
1617 if (res.all_variables() != subpoly.all_variables())
1618 res = gcd(subpoly,res);
1620 res = gcd_linear(subpoly,res);
1623 bracketinfo.pop_back();
1627 if (res.is_zero() || !poly::divides(res,ppa) || !poly::divides(res,ppb)) {
1628 MesPrint(
"Bad gcd found.");
1629 std::cout <<
"Bad gcd:" << res <<
" for " << ppa <<
" " << ppb << std::endl;
1634 res *= gcdconts *
poly(BHEAD res.sign());
1637 cout <<
"*** [" << thetime() <<
"] RES : gcd("<<a<<
","<<b<<
") = "<<res<<endl;
int is_dense_univariate() const
WORD NextPrime(PHEAD WORD)
int GetModInverses(WORD, WORD, WORD *, WORD *)
bool gcd_heuristic_possible(const poly &a)