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
52using 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
61template<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
91const 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
143const 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
205const 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
263const 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
325const 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
374const 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
426const 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
488const 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)
571const 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
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;
602}
603
604// Sets all coefficients to the range 0..modp-1 and the powers of x2...xn to 0
605const 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
619const 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
677const 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
901const 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
1021const 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
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
1190const 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
1311const 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
1342const 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
1370const 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
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
1405const 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
1436const 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
1496const 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*/
Definition poly.h:49
int is_dense_univariate() const
Definition poly.cc:2278
WORD NextPrime(PHEAD WORD)
Definition reken.c:3654
int GetModInverses(WORD, WORD, WORD *, WORD *)
Definition reken.c:1466
bool gcd_heuristic_possible(const poly &a)
Definition polygcd.cc:1144