FORM  4.3
polygcd.cc
Go to the documentation of this file.
1 
6 /* #[ License : */
7 /*
8  * Copyright (C) 1984-2022 J.A.M. Vermaseren
9  * When using this file you are requested to refer to the publication
10  * J.A.M.Vermaseren "New features of FORM" math-ph/0010025
11  * This is considered a matter of courtesy as the development was paid
12  * for by FOM the Dutch physics granting agency and we would like to
13  * be able to track its scientific use to convince FOM of its value
14  * for the community.
15  *
16  * This file is part of FORM.
17  *
18  * FORM is free software: you can redistribute it and/or modify it under the
19  * terms of the GNU General Public License as published by the Free Software
20  * Foundation, either version 3 of the License, or (at your option) any later
21  * version.
22  *
23  * FORM is distributed in the hope that it will be useful, but WITHOUT ANY
24  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
25  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
26  * details.
27  *
28  * You should have received a copy of the GNU General Public License along
29  * with FORM. If not, see <http://www.gnu.org/licenses/>.
30  */
31 /* #] License : */
32 /*
33  #[ include :
34 */
35 
36 #include "poly.h"
37 #include "polygcd.h"
38 
39 #include <iostream>
40 #include <vector>
41 #include <cmath>
42 #include <map>
43 #include <algorithm>
44 
45 //#define DEBUG
46 //#define DEBUGALL
47 
48 #ifdef DEBUG
49 #include "mytime.h"
50 #endif
51 
52 using namespace std;
53 
54 /*
55  #] include :
56  #[ ostream operator :
57 */
58 
59 #ifdef DEBUG
60 // ostream operator for outputting vector<T>s for debugging purposes
61 template<class T> ostream& operator<< (ostream &out, const vector<T> &x) {
62  out<<"{";
63  for (int i=0; i<(int)x.size(); i++) {
64  if (i>0) out<<",";
65  out<<x[i];
66  }
67  out<<"}";
68  return out;
69 }
70 #endif
71 
72 /*
73  #] ostream operator :
74  #[ integer_gcd :
75 */
76 
91 const poly polygcd::integer_gcd (const poly &a, const poly &b) {
92 
93 #ifdef DEBUGALL
94  cout << "*** [" << thetime() << "] CALL: integer_gcd(" << a << "," << b << ")" << endl;
95 #endif
96 
97  POLY_GETIDENTITY(a);
98 
99  if (a.is_zero()) return b;
100  if (b.is_zero()) return a;
101 
102  poly c(BHEAD 0, a.modp, a.modn);
103  WORD nc;
104 
105  GcdLong(BHEAD
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);
109 
110  WORD x = 2 + AN.poly_num_vars + ABS(nc);
111  c[1] = x; // term length
112  c[0] = x+1; // total length
113  c[x] = nc; // coefficient length
114 
115  for (int i=0; i<AN.poly_num_vars; i++)
116  c[2+i] = 0; // powers
117 
118 #ifdef DEBUGALL
119  cout << "*** [" << thetime() << "] RES : integer_gcd(" << a << "," << b << ") = " << c << endl;
120 #endif
121 
122  return c;
123 }
124 
125 /*
126  #] integer_gcd :
127  #[ integer_content :
128 */
129 
143 const poly polygcd::integer_content (const poly &a) {
144 
145 #ifdef DEBUGALL
146  cout << "*** [" << thetime() << "] CALL: integer_content(" << a << ")" << endl;
147 #endif
148 
149  POLY_GETIDENTITY(a);
150 
151  if (a.modp>0) return a.integer_lcoeff();
152 
153  poly c(BHEAD 0, 0, 1);
154  WORD *d = (WORD *)NumberMalloc("polygcd::integer_content");
155  WORD nc=0;
156 
157  for (int i=0; i<AN.poly_num_vars; i++)
158  c[2+i] = 0;
159 
160  for (int i=1; i<a[0]; i+=a[i]) {
161 
162  WCOPY(d,&c[2+AN.poly_num_vars],nc);
163 
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);
167 
168  WORD x = 2 + AN.poly_num_vars + ABS(nc);
169  c[1] = x; // term length
170  c[0] = x+1; // total length
171  c[x] = nc; // coefficient length
172  }
173 
174  if (a.sign() != c.sign()) c *= poly(BHEAD -1);
175 
176  NumberFree(d,"polygcd::integer_content");
177 
178 #ifdef DEBUGALL
179  cout << "*** [" << thetime() << "] RES : integer_content(" << a << ") = " << c << endl;
180 #endif
181 
182  return c;
183 }
184 
185 /*
186  #] integer_content :
187  #[ content_univar :
188 */
189 
205 const poly polygcd::content_univar (const poly &a, int x) {
206 
207 #ifdef DEBUG
208  cout << "*** [" << thetime() << "] CALL: content_univar(" << a << "," << string(1,'a'+x) << ")" << endl;
209 #endif
210 
211  POLY_GETIDENTITY(a);
212 
213  poly res(BHEAD 0, a.modp, a.modn);
214 
215  for (int i=1; i<a[0];) {
216  poly b(BHEAD 0, a.modp, a.modn);
217  int deg = a[i+1+x];
218 
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]);
222  b[b[0]+1+x] = 0;
223  b[0] += a[i];
224  }
225 
226  res = gcd(res, b);
227 
228  if (res.is_integer()) {
229  res = integer_content(a);
230  break;
231  }
232  }
233 
234  if (a.sign() != res.sign()) res *= poly(BHEAD -1);
235 
236 #ifdef DEBUG
237  cout << "*** [" << thetime() << "] RES : content_univar(" << a << "," << string(1,'a'+x) << ") = " << res << endl;
238 #endif
239 
240  return res;
241 }
242 
243 /*
244  #] content_univar :
245  #[ content_multivar :
246 */
247 
263 const poly polygcd::content_multivar (const poly &a, int x) {
264 
265 #ifdef DEBUGALL
266  cout << "*** [" << thetime() << "] CALL: content_multivar(" << a << "," << string(1,'a'+x) << ")" << endl;
267 #endif
268 
269  POLY_GETIDENTITY(a);
270 
271  poly res(BHEAD 0, a.modp, a.modn);
272 
273  for (int i=1,j; i<a[0]; i=j) {
274  poly b(BHEAD 0, a.modp, a.modn);
275 
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]) {
280  same_powers = false;
281  break;
282  }
283  if (!same_powers) break;
284 
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;
289 
290  b[0] += a[j];
291  }
292 
293  res = gcd_Euclidean(res, b);
294 
295  if (res.is_integer()) {
296  res = poly(BHEAD a.sign(),a.modp,a.modn);
297  break;
298  }
299  }
300 
301 #ifdef DEBUGALL
302  cout << "*** [" << thetime() << "] RES : content_multivar(" << a << "," << string(1,'a'+x) << ") = " << res << endl;
303 #endif
304 
305  return res;
306 }
307 
308 /*
309  #] content_multivar :
310  #[ coefficient_list_gcd :
311 */
312 
325 const vector<WORD> polygcd::coefficient_list_gcd (const vector<WORD> &_a, const vector<WORD> &_b, WORD p) {
326 
327 #ifdef DEBUGALL
328  cout << "*** [" << thetime() << "] CALL: coefficient_list_gcd("<<_a<<","<<_b<<","<<p<<")"<<endl;
329 #endif
330 
331  vector<WORD> a(_a), b(_b);
332 
333  while (b.size() != 0) {
334  a = poly::coefficient_list_divmod(a,b,p,1);
335  swap(a,b);
336  }
337 
338  while (a.back()==0) a.pop_back();
339 
340  WORD inv;
341  GetModInverses(a.back() + (a.back()<0?p:0), p, &inv, NULL);
342 
343  for (int i=0; i<(int)a.size(); i++)
344  a[i] = (LONG)inv*a[i] % p;
345 
346 #ifdef DEBUGALL
347  cout << "*** [" << thetime() << "] RES : coefficient_list_gcd("<<_a<<","<<_b<<","<<p<<") = "<<a<<endl;
348 #endif
349 
350  return a;
351 }
352 
353 /*
354  #] coefficient_list_gcd :
355  #[ gcd_Euclidean :
356 */
357 
374 const poly polygcd::gcd_Euclidean (const poly &a, const poly &b) {
375 
376 #ifdef DEBUG
377  cout << "*** [" << thetime() << "] CALL: gcd_Euclidean("<<a<<","<<b<<")"<<endl;
378 #endif
379 
380  POLY_GETIDENTITY(a);
381 
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);
385 
386  poly res(BHEAD 0);
387 
388  if (a.is_dense_univariate()>=-1 && b.is_dense_univariate()>=-1) {
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);
392  }
393  else {
394  res = a;
395  poly rem(b);
396  while (!rem.is_zero())
397  swap(res%=rem, rem);
398  res /= res.integer_lcoeff();
399  }
400 
401 #ifdef DEBUG
402  cout << "*** [" << thetime() << "] RES : gcd_Euclidean("<<a<<","<<b<<") = "<<res<<endl;
403 #endif
404 
405  return res;
406 }
407 
408 /*
409  #] gcd_Euclidean :
410  #[ chinese_remainder :
411 */
412 
426 const poly polygcd::chinese_remainder (const poly &a1, const poly &m1, const poly &a2, const poly &m2) {
427 
428 #ifdef DEBUGALL
429  cout << "*** [" << thetime() << "] CALL: chinese_remainder(" << a1 << "," << m1 << "," << a2 << "," << m2 << ")" << endl;
430 #endif
431 
432  POLY_GETIDENTITY(a1);
433 
434  WORD nx,ny,nz;
435  UWORD *x = (UWORD *)NumberMalloc("polygcd::chinese_remainder");
436  UWORD *y = (UWORD *)NumberMalloc("polygcd::chinese_remainder");
437  UWORD *z = (UWORD *)NumberMalloc("polygcd::chinese_remainder");
438 
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);
442 
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]],
445  y, &ny);
446 
447  MulLong (x,nx,y,ny,z,&nz);
448  MulLong (z,nz,(UWORD *)&m1[2+AN.poly_num_vars],m1[m1[1]],x,&nx);
449 
450  AddLong (x,nx,(UWORD *)&a1[2+AN.poly_num_vars], a1.is_zero() ? 0 : a1[a1[1]],y,&ny);
451 
452  MulLong ((UWORD *)&m1[2+AN.poly_num_vars], m1[m1[1]],
453  (UWORD *)&m2[2+AN.poly_num_vars], m2[m2[1]],
454  (UWORD *)z,&nz);
455 
456  TakeNormalModulus (y,&ny,(UWORD *)z,nz,NOUNPACK);
457 
458  poly res(BHEAD y,ny);
459 
460  NumberFree(x,"polygcd::chinese_remainder");
461  NumberFree(y,"polygcd::chinese_remainder");
462  NumberFree(z,"polygcd::chinese_remainder");
463 
464 #ifdef DEBUGALL
465  cout << "*** [" << thetime() << "] RES : chinese_remainder(" << a1 << "," << m1 << "," << a2 << "," << m2 << ") = " << res << endl;
466 #endif
467 
468  return res;
469 }
470 
471 /*
472  #] chinese_remainder :
473  #[ substitute :
474 */
475 
488 const poly polygcd::substitute(const poly &a, int x, int c) {
489 
490  POLY_GETIDENTITY(a);
491 
492  poly b(BHEAD 0);
493 
494  if (a.is_zero()) {
495  return b;
496  }
497 
498  bool zero=true;
499  int bi=1;
500 
501  // cache size is bounded by the degree in x, twice the number of terms of
502  // the polynomial and a constant
503  vector<WORD> cache(min(a.degree(x)+1,min(2*a.number_of_terms(),
504  POLYGCD_RAISPOWMOD_CACHE_SIZE)), 0);
505 
506  for (int ai=1; ai<=a[0]; ai+=a[ai]) {
507  // last term or different power, then add term to b iff non-zero
508  if (!zero) {
509  bool add=false;
510  if (ai==a[0])
511  add=true;
512  else {
513  for (int i=0; i<AN.poly_num_vars; i++)
514  if (i!=x && a[ai+1+i]!=b[bi+1+i]) {
515  zero=true;
516  add=true;
517  break;
518  }
519  }
520 
521  if (add) {
522  if (b[bi+AN.poly_num_vars+1] < 0)
523  b[bi+AN.poly_num_vars+1] += a.modp;
524  bi+=b[bi];
525  }
526 
527  if (ai==a[0]) break;
528  }
529 
530  b.check_memory(bi);
531 
532  // create new term in b
533  if (zero) {
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];
537  b[bi+1+x] = 0;
538  b[bi+AN.poly_num_vars+1] = 0;
539  b[bi+AN.poly_num_vars+2] = 1;
540  }
541 
542  // add term of a to the current term in b
543  LONG coeff = a[ai+1+AN.poly_num_vars] * a[ai+2+AN.poly_num_vars];
544  int pow = a[ai+1+x];
545 
546  if (pow<(int)cache.size()) {
547  if (cache[pow]==0)
548  cache[pow] = RaisPowMod(c, pow, a.modp);
549  coeff = (coeff * cache[pow]) % a.modp;
550  }
551  else {
552  coeff = (coeff * RaisPowMod(c, pow, a.modp)) % a.modp;
553  }
554 
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;
557  }
558 
559  b[0]=bi;
560  b.setmod(a.modp);
561 
562  return b;
563 }
564 
565 /*
566  #] substitute :
567  #[ sparse_interpolation helper functions :
568 */
569 
570 // Returns a list of size #terms(a) with entries PROD(ci^powi, i=2..n)
571 const vector<int> polygcd::sparse_interpolation_get_mul_list (const poly &a, const vector<int> &x, const vector<int> &c) {
572  // cache size for variable x is bounded by the degree in x, twice
573  // the number of terms of the polynomial and a constant
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);
578 
579  vector<int> res;
580  for (int i=1; i<a[0]; i+=a[i]) {
581  LONG coeff=1;
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;
588  }
589  else {
590  coeff = (coeff * RaisPowMod(c[j], pow, a.modp)) % a.modp;
591  }
592  }
593  res.push_back(coeff);
594  }
595  return res;
596 }
597 
598 // Multiplies the coefficients of a with the entries of mul
599 void 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;
602 }
603 
604 // Sets all coefficients to the range 0..modp-1 and the powers of x2...xn to 0
605 const poly polygcd::sparse_interpolation_reduce_poly (const poly &a, const vector<int> &x) {
606  poly res(a);
607  for (int i=1; i<a[0]; i+=a[i]) {
608  for (int j=1; j<(int)x.size(); j++)
609  res[i+1+x[j]]=0;
610  if (res[i+a[i]-1]==-1) {
611  res[i+a[i]-1]=1;
612  res[i+a[i]-2]=a.modp-res[i+a[i]-2];
613  }
614  }
615  return res;
616 }
617 
618 // Collects entries with equal powers, so that the result is a proper polynomial
619 const poly polygcd::sparse_interpolation_fix_poly (const poly &a, int x) {
620 
621  POLY_GETIDENTITY(a);
622  poly res(BHEAD 0,a.modp,1);
623 
624  int j=1;
625  bool newterm=true;
626 
627  for (int i=1; i<a[0]; i+=a[i]) {
628  if (newterm)
629  res.termscopy(&a[i], j, a[i]);
630  else
631  res[j+res[j]-2] = ((LONG)res[j+res[j]-2] + a[i+a[i]-2]) % a.modp;
632 
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];
635  }
636 
637  res[0]=j;
638  return res;
639 }
640 
641 /*
642  #] sparse_interpolation helper functions :
643  #[ gcd_modular_sparse_interpolation :
644 */
645 
677 const poly polygcd::gcd_modular_sparse_interpolation (const poly &origa, const poly &origb, const vector<int> &x, const poly &s) {
678 
679 #ifdef DEBUG
680  cout << "*** [" << thetime() << "] CALL: gcd_modular_sparse_interpolation("
681  << origa << "," << origb << "," << x << "," << "," << s <<")" << endl;
682 #endif
683 
684  POLY_GETIDENTITY(origa);
685 
686  // strip multivariate content
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;
692 
693  // for non-monic cases, we need to normalize with the gcd of the lcoeffs of a poly in x[0]
694  // or else the shape fitting does not work.
695  // FIXME: the current implementation still rejects some valid shapes.
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));
699  }
700 
701  // reduce polynomials
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));
706 
707  // set all coefficients to 1
708  int N=0;
709  for (int i=1; i<sred[0]; i+=sred[i]) {
710  sred[i+sred[i]-2] = sred[i+sred[i]-1] = 1;
711  N++;
712  }
713 
714  // generate random numbers and check there this set doesn't result
715  // in a singular matrix
716  vector<int> c(x.size()-1);
717  vector<int> smul;
718 
719  bool duplicates;
720  do {
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);
724 
725  duplicates = false;
726 
727  int fr=0,to=0;
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])
734  duplicates = true;
735  fr=to;
736  }
737  }
738  while (duplicates);
739 
740  // get the lists to multiply the polynomials with every iteration
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));
744 
745  vector<vector<vector<LONG> > > M;
746  vector<vector<LONG> > V;
747 
748  int maxMsize=0;
749 
750  // create (empty) matrices
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>());
755  }
756  M.back().push_back(vector<LONG>());
757  V.back().push_back(0);
758  maxMsize = max(maxMsize, (int)M.back().size());
759  }
760 
761  // generate linear equations
762  for (int numg=0; numg<maxMsize; numg++) {
763 
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]));
767 
768  // A fix for non-monic gcds. This could be slow if lmodI has many terms,
769  // since it overfits the gcd now. Another gcd has to be run to remove the
770  // extra terms.
771  poly gcd(lmodI * gcd_Euclidean(amodI,bmodI));
772 
773  // if correct gcd
774  if (!gcd.is_zero() && gcd[2+x[0]]==sred[2+x[0]]) {
775 
776  // for each power in the gcd, generate an equation if needed
777  int gi=1, midx=0;
778 
779  for (int si=1; si<s[0];) {
780  // if the term exists, set Vi=coeff, otherwise Vi remains 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];
784  gi += gcd[gi];
785  }
786 
787  // add the coefficients of s to the matrix M
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]);
791  si += s[si];
792  }
793 
794  midx++;
795  }
796  }
797  else {
798  // incorrect gcd
799  if (!gcd.is_zero() && gcd[2+x[0]]<sred[2+x[0]])
800  return poly(BHEAD 0);
801  numg--;
802  }
803 
804  // multiply polynomials by the lists to obtain new ones
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);
809  }
810 
811  // solve the linear equations
812  for (int i=0; i<(int)M.size(); i++) {
813  int n = M[i].size();
814 
815  // Gaussian elimination
816  for (int j=0; j<n; j++) {
817  for (int k=0; k<j; k++) {
818  LONG x = M[i][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;
822  }
823 
824  // normalize row
825  WORD x = M[i][j][j]; // WORD for GetModInverses
826  GetModInverses(x + (x<0?a.modp:0), a.modp, &x, NULL);
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;
830  }
831 
832  // solve
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;
836  }
837 
838  // create coefficient list
839  vector<LONG> coeff;
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]);
843 
844  // create resulting polynomial
845  poly res(BHEAD 0);
846  int ri=1, i=0;
847  for (int si=1; si<s[0]; si+=s[si]) {
848  res.check_memory(ri);
849  res[ri] = 3 + AN.poly_num_vars; // term length
850  for (int j=0; j<AN.poly_num_vars; j++)
851  res[ri+1+j] = s[si+1+j]; // powers
852  res[ri+1+AN.poly_num_vars] = ABS(coeff[i]); // coefficient
853  res[ri+2+AN.poly_num_vars] = SGN(coeff[i]); // coefficient length
854  i++;
855  ri += res[ri];
856  }
857  res[0]=ri; // total length
858  res.setmod(a.modp,1);
859 
860  if (!poly::divides(res, lcgcd * a) || !poly::divides(res, lcgcd * b)) {
861  return poly(BHEAD 0); // bad shape
862  } else {
863  // refine gcd
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));
868  }
869 
870 #ifdef DEBUG
871  cout << "*** [" << thetime() << "] RES : gcd_modular_sparse_interpolation("
872  << a << "," << b << "," << x << "," << "," << s <<") = " << res << endl;
873 #endif
874 
875  return gcdconts * res;
876 }
877 
878 /*
879  #] gcd_modular_sparse_interpolation :
880  #[ gcd_modular_dense_interpolation :
881 */
882 
901 const poly polygcd::gcd_modular_dense_interpolation (const poly &a, const poly &b, const vector<int> &x, const poly &s) {
902 
903 #ifdef DEBUG
904  cout << "*** [" << thetime() << "] CALL: gcd_modular_dense_interpolation(" << a << "," << b << "," << x << "," << s <<")" << endl;
905 #endif
906 
907  POLY_GETIDENTITY(a);
908 
909  // if univariate, then use Euclidean algorithm
910  if (x.size() == 1) {
911  return gcd_Euclidean(a,b);
912  }
913 
914  // if shape is known, use sparse interpolation
915  if (!s.is_zero()) {
916  poly res = gcd_modular_sparse_interpolation (a,b,x,s);
917  if (!res.is_zero()) return res;
918  // apparently the shape was not correct. continue.
919  }
920 
921  // divide out multivariate content in last variable
922  int X = x.back();
923 
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);
929 
930  // gcd of leading coefficients
931  poly lcoeffa(ppa.lcoeff_multivar(X));
932  poly lcoeffb(ppb.lcoeff_multivar(X));
933  poly gcdlcoeffs(gcd_Euclidean(lcoeffa,lcoeffb));
934 
935  // calculate the degree bound for each variable
936  int m = MiN(ppa.degree(x[x.size() - 2]),ppb.degree(x[x.size() - 2]));
937 
938  poly res(BHEAD 0);
939  poly oldres(BHEAD 0);
940  poly newshape(BHEAD 0);
941  poly modpoly(BHEAD 1,a.modp);
942 
943  while (true) {
944  // generate random constants and substitute it
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;
948 
949  poly amodc(substitute(ppa,X,c));
950  poly bmodc(substitute(ppb,X,c));
951 
952  // calculate gcd recursively
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]);
955 
956  // normalize
957  gcdmodc = (gcdmodc * substitute(gcdlcoeffs,X,c)) / gcdmodc.integer_lcoeff();
958  poly simple(poly::simple_poly(BHEAD X,c,1,a.modp)); // (X-c) mod p
959 
960  // if power is smaller, the old one was wrong
961  if ((res.is_zero() && n == m) || n < m) {
962  m = n;
963  res = gcdmodc;
964  newshape = gcdmodc; // set a new shape (interpolation does not change it)
965  modpoly = simple;
966  }
967  else if (n == m) {
968  oldres = res;
969  // equal powers, so interpolate results
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;
973 
974  GetModInverses(coeff_word, a.modp, &coeff_word, NULL);
975 
976  res.setmod(a.modp); // make sure the mod is set before substituting
977  res += poly(BHEAD coeff_word, a.modp, 1) * modpoly * (gcdmodc - substitute(res,X,c));
978  modpoly *= simple;
979  }
980 
981  // check whether this is the complete gcd
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)) {
985 #ifdef DEBUG
986  cout << "*** [" << thetime() << "] RES : gcd_modular_dense_interpolation(" << a << "," << b << ","
987  << x << "," << "," << s <<") = " << gcdconts * nres << endl;
988 #endif
989  return gcdconts * nres;
990  }
991 
992  // At this point, the gcd may be too large due to bad luck
993  // TODO: create an efficient fail state that tries to find a smaller
994  // polynomial without interpolating bad ones?
995  newshape = poly(BHEAD 0); // reset the shape, important!
996  }
997  }
998 }
999 
1000 /*
1001  #] gcd_modular_dense_interpolation : :
1002  #[ gcd_modular :
1003 */
1004 
1021 const poly polygcd::gcd_modular (const poly &origa, const poly &origb, const vector<int> &x) {
1022 
1023 #ifdef DEBUG
1024  cout << "*** [" << thetime() << "] CALL: gcd_modular(" << origa << "," << origb << "," << x << ")" << endl;
1025 #endif
1026 
1027  POLY_GETIDENTITY(origa);
1028 
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();
1032 
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());
1039 
1040  int pnum=0;
1041 
1042  poly d(BHEAD 0);
1043  poly m1(BHEAD 1);
1044  int mindeg=MAXPOSITIVE;
1045 
1046  while (true) {
1047  // choose a prime and solve modulo the prime
1048  WORD p = NextPrime(BHEAD pnum++);
1049  if (poly(a.integer_lcoeff(),p).is_zero()) continue;
1050  if (poly(b.integer_lcoeff(),p).is_zero()) continue;
1051 
1052  poly c(gcd_modular_dense_interpolation(poly(a,p),poly(b,p),x,poly(d,p)));
1053  c = (c * poly(g,p)) / c.integer_lcoeff(); // normalize so that lcoeff(c) = g mod p
1054 
1055  if (c.is_zero()) {
1056  // unlucky choices somewhere, so start all over again
1057  d = poly(BHEAD 0);
1058  m1 = poly(BHEAD 1);
1059  mindeg = MAXPOSITIVE;
1060  continue;
1061  }
1062 
1063  if (!(poly(a,p)%c).is_zero()) continue;
1064  if (!(poly(b,p)%c).is_zero()) continue;
1065 
1066  int deg = c.degree(x[0]);
1067 
1068  if (deg < mindeg) {
1069  // small degree, so the old one is wrong
1070  d=c;
1071  d.modp=a.modp;
1072  d.modn=a.modn;
1073  m1 = poly(BHEAD p);
1074  mindeg=deg;
1075  }
1076  else if (deg == mindeg) {
1077  // same degree, so use Chinese Remainder Algorithm
1078  poly newd(BHEAD 0);
1079 
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);
1083 
1084  newd.check_memory(newd[0]);
1085 
1086  if (comp <= 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]);
1089  di+=d[di];
1090  }
1091  if (comp >= 0) {
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]);
1094  ci+=c[ci];
1095  }
1096 
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]];
1101  }
1102 
1103  m1 *= poly(BHEAD p);
1104  d=newd;
1105  }
1106 
1107  // divide out spurious integer content
1108  poly ppd(d / integer_content(d));
1109 
1110  // check whether this is the complete gcd
1111  if (poly::divides(ppd,a) && poly::divides(ppd,b)) {
1112 #ifdef DEBUG
1113  cout << "*** [" << thetime() << "] RES : gcd_modular(" << origa << "," << origb << "," << x << ") = "
1114  << ic * ppd << endl;
1115 #endif
1116  return ic * ppd;
1117  }
1118 #ifdef DEBUG
1119  cout << "*** [" << thetime() << "] Retrying modular_gcd with new prime" << endl;
1120 #endif
1121  }
1122 }
1123 
1124 /*
1125  #] gcd_modular :
1126  #[ gcd_heuristic_possible :
1127 */
1128 
1144 bool gcd_heuristic_possible (const poly &a) {
1145 
1146  POLY_GETIDENTITY(a);
1147 
1148  double prod_deg = 1;
1149  for (int j=0; j<AN.poly_num_vars; j++)
1150  prod_deg *= a[2+j]+1;
1151 
1152  double digits = ABS(a[1+a[1]-1]);
1153  double lead = a[1+1+AN.poly_num_vars];
1154 
1155  return prod_deg*(digits-1+log(2*ABS(lead))/log(2.0)/(BITSINWORD/2)) < POLYGCD_HEURISTIC_MAX_DIGITS;
1156 }
1157 
1158 /*
1159  #] gcd_heuristic_possible :
1160  #[ gcd_heuristic :
1161 */
1162 
1163 
1190 const poly polygcd::gcd_heuristic (const poly &a, const poly &b, const vector<int> &x, int max_tries) {
1191 
1192 #ifdef DEBUG
1193  cout << "*** [" << thetime() << "] CALL: gcd_heuristic("<<a<<","<<b<<","<<x<<")\n";
1194 #endif
1195 
1196  if (a.is_integer()) return integer_gcd(a,integer_content(b));
1197  if (b.is_integer()) return integer_gcd(integer_content(a),b);
1198 
1199  POLY_GETIDENTITY(a);
1200 
1201  // Calculate xi = 2*min(max(coefficients a),max(coefficients b))+2
1202 
1203  UWORD *pxi=NULL;
1204  WORD nxi=0;
1205 
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];
1210  nxi = na;
1211  }
1212  }
1213 
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];
1218  nxi = nb;
1219  }
1220  }
1221 
1222  poly xi(BHEAD pxi,nxi);
1223 
1224  // Addition of another random factor gives better performance
1225  xi = xi*poly(BHEAD 2) + poly(BHEAD 2 + wranf(BHEAD0)%POLYGCD_HEURISTIC_MAX_ADD_RANDOM);
1226 
1227  // If degree*digits(xi) is too large, throw exception
1228  if (max(a.degree(x[0]),b.degree(x[0])) * xi[xi[1]] >= MiN(AM.MaxTal, POLYGCD_HEURISTIC_MAX_DIGITS)) {
1229 #ifdef DEBUG
1230  cout << "*** [" << thetime() << "] RES : gcd_heuristic("<<a<<","<<b<<","<<x<<") = overflow\n";
1231 #endif
1232  throw(gcd_heuristic_failed());
1233  }
1234 
1235  for (int times=0; times<max_tries; times++) {
1236 
1237  // Recursively calculate the gcd
1238 
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));
1242 
1243  // If a gcd is found, reconstruct the powers of x
1244  if (!gamma.is_zero()) {
1245  // res is construct is reverse order. idx/len are for reversing
1246  // it in the correct order
1247  poly res(BHEAD 0), c(BHEAD 0);
1248  vector<int> idx, len;
1249 
1250  for (int power=0; !gamma.is_zero(); power++) {
1251 
1252  // calculate c = gamma % xi (c and gamma are polynomials, xi is integer)
1253  c = gamma;
1254  c.coefficients_modulo((UWORD *)&xi[2+AN.poly_num_vars], xi[xi[0]-1], false);
1255 
1256  // Add the terms c * x^power to res
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;
1261 
1262  // Store idx/len for reversing
1263  if (!c.is_zero()) {
1264  idx.push_back(res[0]);
1265  len.push_back(c[0]-1);
1266  }
1267 
1268  res[0] += c[0]-1;
1269 
1270  // Divide gamma by xi
1271  gamma = (gamma - c) / xi;
1272  }
1273 
1274  // Reverse the resulting polynomial
1275  poly rev(BHEAD 0);
1276  rev.check_memory(res[0]);
1277 
1278  rev[0] = 1;
1279  for (int i=idx.size()-1; i>=0; i--) {
1280  rev.termscopy(&res[idx[i]], rev[0], len[i]);
1281  rev[0] += len[i];
1282  }
1283  res = rev;
1284 
1285  poly ppres = res / integer_content(res);
1286 
1287  if ((a%ppres).is_zero() && (b%ppres).is_zero()) {
1288 #ifdef DEBUG
1289  cout << "*** [" << thetime() << "] RES : gcd_heuristic("<<a<<","<<b<<","<<x<<") = "<<res<<"\n";
1290 #endif
1291  return res;
1292  }
1293  }
1294 
1295  // Next xi by multiplying with the golden ratio to avoid correlated errors
1296  xi = xi * poly(BHEAD 28657) / poly(BHEAD 17711) + poly(BHEAD wranf(BHEAD0) % POLYGCD_HEURISTIC_MAX_ADD_RANDOM);
1297  }
1298 
1299 #ifdef DEBUG
1300  cout << "*** [" << thetime() << "] RES : gcd_heuristic("<<a<<","<<b<<","<<x<<") = failed\n";
1301 #endif
1302 
1303  return poly(BHEAD 0);
1304 }
1305 
1306 /*
1307  #] gcd_heuristic :
1308  #[ bracket :
1309 */
1310 
1311 const map<vector<int>,poly> polygcd::full_bracket(const poly &a, const vector<int>& filter) {
1312  POLY_GETIDENTITY(a);
1313 
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];
1320 
1321  // create monomial
1322  poly mon(BHEAD 1);
1323  mon.setmod(a.modp);
1324  mon[0] = a[ai] + 1;
1325  for (int i=0; i<a[ai]; i++)
1326  if (i > 0 && i <= AN.poly_num_vars && varpattern[i - 1])
1327  mon[1+i] = 0;
1328  else
1329  mon[1+i] = a[ai+i];
1330 
1331  map<vector<int>,poly>::iterator i = bracket.find(varpattern);
1332  if (i == bracket.end()) {
1333  bracket.insert(std::make_pair(varpattern, mon));
1334  } else {
1335  i->second += mon;
1336  }
1337  }
1338 
1339  return bracket;
1340 }
1341 
1342 const poly polygcd::bracket(const poly &a, const vector<int>& pattern, const vector<int>& filter) {
1343  POLY_GETIDENTITY(a);
1344 
1345  poly bracket(BHEAD 0);
1346  for (int ai=1; ai<a[0]; ai+=a[ai]) {
1347  bool ispat = true;
1348  for (int i=0; i<AN.poly_num_vars; i++)
1349  if (filter[i] == 1 && pattern[i] != a[ai + i + 1]) {
1350  ispat = false;
1351  break;
1352  }
1353 
1354  if (ispat) {
1355  poly mon(BHEAD 1);
1356  mon.setmod(a.modp);
1357  mon[0] = a[ai] + 1;
1358  for (int i=0; i<a[ai]; i++)
1359  if (i > 0 && i <= AN.poly_num_vars && pattern[i - 1])
1360  mon[1+i] = 0;
1361  else
1362  mon[1+i] = a[ai+i];
1363  bracket += mon;
1364  }
1365  }
1366 
1367  return bracket;
1368 }
1369 
1370 const map<vector<int>,int> polygcd::bracket_count(const poly &a, const vector<int>& filter) {
1371  POLY_GETIDENTITY(a);
1372 
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];
1379 
1380  map<vector<int>,int>::iterator i = bracket.find(varpattern);
1381  if (i == bracket.end()) {
1382  bracket.insert(std::make_pair(varpattern, 0));
1383  } else {
1384  i->second++;
1385  }
1386  }
1387 
1388  return bracket;
1389 }
1390 
1391 struct BracketInfo {
1392  std::vector<int> pattern;
1393  int num_terms, dummy;
1394  const poly* p;
1395 
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; } // biggest should be first!
1398 };
1399 
1400 /*
1401  #] bracket :
1402  #[ gcd_linear:
1403 */
1404 
1405 const poly gcd_linear_helper (const poly &a, const poly &b) {
1406  POLY_GETIDENTITY(a);
1407 
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);
1411  filter[i] = 1;
1412 
1413  // bracket the linear variable
1414  map<vector<int>,poly> ba = polygcd::full_bracket(a, filter);
1415 
1416  poly subgcd(BHEAD 1);
1417  if (ba.size() == 2)
1418  subgcd = polygcd::gcd_linear(ba.begin()->second, (++ba.begin())->second);
1419  else
1420  subgcd = ba.begin()->second;
1421 
1422  poly linfac = a / subgcd;
1423  if (poly::divides(linfac,b))
1424  return linfac * polygcd::gcd_linear(subgcd, b / linfac);
1425 
1426  return polygcd::gcd_linear(subgcd, b);
1427  }
1428 
1429  return poly(BHEAD 0);
1430 }
1431 
1436 const poly polygcd::gcd_linear (const poly &a, const poly &b) {
1437 #ifdef DEBUG
1438  cout << "*** [" << thetime() << "] CALL: gcd_linear("<<a<<","<<b<<")\n";
1439 #endif
1440 
1441  POLY_GETIDENTITY(a);
1442 
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();
1446 
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);
1450  }
1451 
1452  poly h = gcd_linear_helper(a, b);
1453  if (!h.is_zero()) return h;
1454 
1455  h = gcd_linear_helper(b, a);
1456  if (!h.is_zero()) return h;
1457 
1458  vector<int> xa = a.all_variables();
1459  vector<int> xb = b.all_variables();
1460 
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]]++;
1464  vector<int> x;
1465  for (int i=0; i<AN.poly_num_vars; i++)
1466  if (used[i]) x.push_back(i);
1467 
1468  return gcd_modular(a,b,x);
1469 }
1470 
1471 /*
1472  #] gcd_linear:
1473  #[ gcd :
1474 */
1475 
1496 const poly polygcd::gcd (const poly &a, const poly &b) {
1497 
1498 #ifdef DEBUG
1499  cout << "*** [" << thetime() << "] CALL: gcd("<<a<<","<<b<<")\n";
1500 #endif
1501 
1502  POLY_GETIDENTITY(a);
1503 
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();
1507 
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);
1511  }
1512 
1513  // Generate a list of variables of a and b
1514  vector<int> xa = a.all_variables();
1515  vector<int> xb = b.all_variables();
1516 
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]]++;
1520  vector<int> x;
1521  for (int i=0; i<AN.poly_num_vars; i++)
1522  if (used[i]) x.push_back(i);
1523 
1524  if (a.is_monomial() || b.is_monomial()) {
1525 
1526  poly res(BHEAD 1,a.modp,a.modn);
1527  if (a.modp == 0) res = integer_gcd(integer_content(a),integer_content(b));
1528 
1529  for (int i=0; i<(int)x.size(); i++)
1530  res[2+x[i]] = 1<<(BITSINWORD-2);
1531 
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]]);
1535 
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]]);
1539 
1540  return res;
1541  }
1542 
1543  // Calculate the contents, their gcd and the primitive parts
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);
1549 
1550  if (ppa == ppb)
1551  return ppa * gcdconts;
1552 
1553  poly res(BHEAD 0);
1554 
1555 #ifdef POLYGCD_USE_HEURISTIC
1556  // Try the heuristic gcd algorithm
1557  if (a.modp==0 && gcd_heuristic_possible(a) && gcd_heuristic_possible(b)) {
1558  try {
1559  res = gcd_heuristic(ppa,ppb,x);
1560  if (!res.is_zero()) res /= integer_content(res);
1561  }
1562  catch (gcd_heuristic_failed) {}
1563  }
1564 #endif
1565 
1566  // If gcd==0, the heuristic algorithm failed, so we do more extensive checks.
1567  // First, we filter out variables that appear in only one of the expressions.
1568  if (res.is_zero()) {
1569  bool unusedVars = false;
1570  for (unsigned int i = 0; i < used.size(); i++) {
1571  if (used[i] == 1) {
1572  unusedVars = true;
1573  break;
1574  }
1575  }
1576 
1577  // if there are no unused variables, go to the linear routine directly
1578  if (!unusedVars) {
1579  res = gcd_linear(ppa,ppb);
1580 #ifdef DEBUG
1581  cout << "New GCD attempt (unused vars): " << res << endl;
1582 #endif
1583  }
1584 
1585  // if res is not the gcd, it is 0 or larger than the gcd.
1586  // we bracket the expression in all the variables that appear only in one expr.
1587  // and we refine the gcd.
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;
1592 
1593  if (!diva) {
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));
1597  }
1598 
1599  if (!divb) {
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));
1603  }
1604 
1605  // sort so that the smallest bracket will be last
1606  sort(bracketinfo.begin(), bracketinfo.end());
1607 
1608  if (res.is_zero()) {
1609  res = bracket(*bracketinfo.back().p, bracketinfo.back().pattern, used);
1610  bracketinfo.pop_back();
1611  }
1612 
1613  while (bracketinfo.size() > 0) {
1614  poly subpoly(bracket(*bracketinfo.back().p, bracketinfo.back().pattern, used));
1615  if (!poly::divides(res,subpoly)) {
1616  // if we can filter out more variables, call gcd again
1617  if (res.all_variables() != subpoly.all_variables())
1618  res = gcd(subpoly,res);
1619  else
1620  res = gcd_linear(subpoly,res);
1621  }
1622 
1623  bracketinfo.pop_back();
1624  }
1625  }
1626 
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;
1630  Terminate(1);
1631  }
1632  }
1633 
1634  res *= gcdconts * poly(BHEAD res.sign());
1635 
1636 #ifdef DEBUG
1637  cout << "*** [" << thetime() << "] RES : gcd("<<a<<","<<b<<") = "<<res<<endl;
1638 #endif
1639 
1640  return res;
1641 }
1642 
1643 /*
1644  #] gcd :
1645 */
int is_dense_univariate() const
Definition: poly.cc:2278
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition: reken.c:1466
Definition: poly.h:49
WORD NextPrime(PHEAD WORD)
Definition: reken.c:3654
bool gcd_heuristic_possible(const poly &a)
Definition: polygcd.cc:1144