polygcd.cc

Go to the documentation of this file.

 
/* #[ License : */
/*
 *   Copyright (C) 1984-2026 J.A.M. Vermaseren
 *   When using this file you are requested to refer to the publication
 *   J.A.M.Vermaseren "New features of FORM" math-ph/0010025
 *   This is considered a matter of courtesy as the development was paid
 *   for by FOM the Dutch physics granting agency and we would like to
 *   be able to track its scientific use to convince FOM of its value
 *   for the community.
 *
 *   This file is part of FORM.
 *
 *   FORM is free software: you can redistribute it and/or modify it under the
 *   terms of the GNU General Public License as published by the Free Software
 *   Foundation, either version 3 of the License, or (at your option) any later
 *   version.
 *
 *   FORM is distributed in the hope that it will be useful, but WITHOUT ANY
 *   WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 *   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
 *   details.
 *
 *   You should have received a copy of the GNU General Public License along
 *   with FORM.  If not, see <http://www.gnu.org/licenses/>.
 */
/* #] License : */ 
/*
    #[ include :
*/
 
#include "poly.h"
#include "polygcd.h"
 
#include <iostream>
#include <vector>
#include <cmath>
#include <map>
#include <algorithm>
 
//#define DEBUG
//#define DEBUGALL
 
#ifdef DEBUG
#include "mytime.h"
#endif
 
using namespace std;
 
/*
    #] include : 
    #[ ostream operator :
*/
 
#ifdef DEBUG
// ostream operator for outputting vector<T>s   for debugging purposes
template<class T> ostream& operator<< (ostream &out, const vector<T> &x) {
    out<<"{";
    for (int i=0; i<(int)x.size(); i++) {
        if (i>0) out<<",";
        out<<x[i];
    }
    out<<"}";
    return out;
}
#endif
 
/*
    #] ostream operator : 
    #[ integer_gcd :
*/
 
const poly polygcd::integer_gcd (const poly &a, const poly &b) {
 
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  CALL: integer_gcd(" << a << "," << b << ")" << endl;
#endif
 
    POLY_GETIDENTITY(a);
    
    if (a.is_zero()) return b;
    if (b.is_zero()) return a;
 
    poly c(BHEAD 0, a.modp, a.modn);
    WORD nc;
    
    GcdLong(BHEAD
                    (UWORD *)&a[AN.poly_num_vars+2],a[a[0]-1],
                    (UWORD *)&b[AN.poly_num_vars+2],b[b[0]-1],
                    (UWORD *)&c[AN.poly_num_vars+2],&nc);
 
    WORD x = 2 + AN.poly_num_vars + ABS(nc);
    c[1] = x;     // term length
    c[0] = x+1;   // total length
    c[x] = nc;    // coefficient length
 
    for (int i=0; i<AN.poly_num_vars; i++)
        c[2+i] = 0; // powers
 
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  RES : integer_gcd(" << a << "," << b << ") = " << c << endl;
#endif
 
    return c;
}
 
/*
    #] integer_gcd : 
    #[ integer_content :
*/
 
const poly polygcd::integer_content (const poly &a) {
 
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  CALL: integer_content(" << a << ")" << endl;
#endif
 
    POLY_GETIDENTITY(a);
 
    if (a.modp>0) return a.integer_lcoeff();    
 
    poly c(BHEAD 0, 0, 1);
    WORD *d = (WORD *)NumberMalloc("polygcd::integer_content");
    WORD nc=0;
 
    for (int i=0; i<AN.poly_num_vars; i++)
        c[2+i] = 0;
 
    for (int i=1; i<a[0]; i+=a[i]) {
 
        WCOPY(d,&c[2+AN.poly_num_vars],nc);
        
        GcdLong(BHEAD (UWORD *)d, nc,
                        (UWORD *)&a[i+1+AN.poly_num_vars], a[i+a[i]-1],
                        (UWORD *)&c[2+AN.poly_num_vars], &nc);
 
        WORD x = 2 + AN.poly_num_vars + ABS(nc);
        c[1] = x;   // term length
        c[0] = x+1; // total length
        c[x] = nc;  // coefficient length
    }   
 
    if (a.sign() != c.sign()) c *= poly(BHEAD -1);
    
    NumberFree(d,"polygcd::integer_content");
    
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  RES : integer_content(" << a << ") = " << c << endl;
#endif
    
    return c;
}
 
/*
    #] integer_content : 
    #[ content_univar :
*/
 
const poly polygcd::content_univar (const poly &a, int x) {
    
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  CALL: content_univar(" << a << "," << string(1,'a'+x) << ")" << endl;
#endif
    
    POLY_GETIDENTITY(a);
    
    poly res(BHEAD 0, a.modp, a.modn);
 
    for (int i=1; i<a[0];) {
        poly b(BHEAD 0, a.modp, a.modn);
        int deg = a[i+1+x];
        
        for (; i<a[0] && a[i+1+x]==deg; i+=a[i]) {
            b.check_memory(b[0]+a[i]);
            b.termscopy(&a[i],b[0],a[i]);
            b[b[0]+1+x] = 0;
            b[0] += a[i];
        }           
        
        res = gcd(res, b);
 
        if (res.is_integer()) {
            res = integer_content(a);
            break;
        }
    }   
 
    if (a.sign() != res.sign()) res *= poly(BHEAD -1);
    
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  RES : content_univar(" << a << "," << string(1,'a'+x) << ") = " << res << endl;
#endif
 
    return res;
}
 
/*
    #] content_univar : 
        #[ content_multivar :
*/
 
const poly polygcd::content_multivar (const poly &a, int x) {
    
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  CALL: content_multivar(" << a << "," << string(1,'a'+x) << ")" << endl;
#endif
 
    POLY_GETIDENTITY(a);
 
    poly res(BHEAD 0, a.modp, a.modn);
 
    for (int i=1,j; i<a[0]; i=j) {
        poly b(BHEAD 0, a.modp, a.modn);
 
        for (j=i; j<a[0]; j+=a[j]) {
            bool same_powers = true;
            for (int k=0; k<AN.poly_num_vars; k++)
                if (k!=x && a[i+1+k]!=a[j+1+k]) {
                    same_powers = false;
                    break;
                }
            if (!same_powers) break;
            
            b.check_memory(b[0]+a[j]);
            b.termscopy(&a[j],b[0],a[j]);
            for (int k=0; k<AN.poly_num_vars; k++)
                if (k!=x) b[b[0]+1+k]=0;
            
            b[0] += a[j];
        }           
 
        res = gcd_Euclidean(res, b);
        
        if (res.is_integer()) {
            res = poly(BHEAD a.sign(),a.modp,a.modn); 
            break;
        }
    }   
    
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  RES : content_multivar(" << a << "," << string(1,'a'+x) << ") = " << res << endl;
#endif
 
    return res;
}
 
/*
    #] content_multivar : 
        #[ coefficient_list_gcd :
*/
 
const vector<WORD> polygcd::coefficient_list_gcd (const vector<WORD> &_a, const vector<WORD> &_b, WORD p) {
 
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  CALL: coefficient_list_gcd("<<_a<<","<<_b<<","<<p<<")"<<endl;
#endif
    
    vector<WORD> a(_a), b(_b);
    
    while (b.size() != 0) {
      a = poly::coefficient_list_divmod(a,b,p,1);
        swap(a,b);
    }
 
    while (a.back()==0) a.pop_back();
    
    WORD inv;
    GetModInverses(a.back() + (a.back()<0?p:0), p, &inv, NULL);
    
    for (int i=0; i<(int)a.size(); i++)
        a[i] = (LONG)inv*a[i] % p;
 
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  RES : coefficient_list_gcd("<<_a<<","<<_b<<","<<p<<") = "<<a<<endl;
#endif
    
    return a;
}
 
/*      
        #] coefficient_list_gcd : 
    #[ gcd_Euclidean :
*/
 
const poly polygcd::gcd_Euclidean (const poly &a, const poly &b) {
 
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  CALL: gcd_Euclidean("<<a<<","<<b<<")"<<endl;
#endif
 
    POLY_GETIDENTITY(a);
    
    if (a.is_zero()) return b;
    if (b.is_zero()) return a;
    if (a.is_integer() || b.is_integer())   return integer_gcd(a,b);
 
    poly res(BHEAD 0);
    
    if (a.is_dense_univariate()>=-1 && b.is_dense_univariate()>=-1) {
        vector<WORD> coeff = coefficient_list_gcd(poly::to_coefficient_list(a),
                                                                                            poly::to_coefficient_list(b), a.modp);
        res = poly::from_coefficient_list(BHEAD coeff, a.first_variable(), a.modp);
    }
    else {
        res = a;
        poly rem(b);
        while (!rem.is_zero()) 
            swap(res%=rem, rem);
        res /= res.integer_lcoeff();
    }
 
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  RES : gcd_Euclidean("<<a<<","<<b<<") = "<<res<<endl;
#endif
 
    return res;
}
 
/*
    #] gcd_Euclidean : 
        #[ chinese_remainder :
*/
 
const poly polygcd::chinese_remainder (const poly &a1, const poly &m1, const poly &a2, const poly &m2) {
 
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  CALL: chinese_remainder(" << a1 << "," << m1 << "," << a2 << "," << m2 << ")" << endl;
#endif
 
    POLY_GETIDENTITY(a1);
    
    WORD nx,ny,nz;
    UWORD *x = (UWORD *)NumberMalloc("polygcd::chinese_remainder");
    UWORD *y = (UWORD *)NumberMalloc("polygcd::chinese_remainder");
    UWORD *z = (UWORD *)NumberMalloc("polygcd::chinese_remainder");
 
    GetLongModInverses(BHEAD (UWORD *)&m1[2+AN.poly_num_vars], m1[m1[1]],
                                         (UWORD *)&m2[2+AN.poly_num_vars], m2[m2[1]],
                                         (UWORD *)x, &nx, NULL, NULL);
    
    AddLong((UWORD *)&a2[2+AN.poly_num_vars], a2.is_zero() ? 0 :  a2[a2[1]],
                    (UWORD *)&a1[2+AN.poly_num_vars], a1.is_zero() ? 0 : -a1[a1[1]],
                    y, &ny);
 
    MulLong (x,nx,y,ny,z,&nz);
    MulLong (z,nz,(UWORD *)&m1[2+AN.poly_num_vars],m1[m1[1]],x,&nx);
 
    AddLong (x,nx,(UWORD *)&a1[2+AN.poly_num_vars], a1.is_zero() ? 0 : a1[a1[1]],y,&ny);
    
    MulLong ((UWORD *)&m1[2+AN.poly_num_vars], m1[m1[1]],
                     (UWORD *)&m2[2+AN.poly_num_vars], m2[m2[1]],
                     (UWORD *)z,&nz);
 
    TakeNormalModulus (y,&ny,(UWORD *)z,nz,NOUNPACK);
    
    poly res(BHEAD y,ny);
 
    NumberFree(x,"polygcd::chinese_remainder");
    NumberFree(y,"polygcd::chinese_remainder");
    NumberFree(z,"polygcd::chinese_remainder");
 
#ifdef DEBUGALL
    cout << "*** [" << thetime() << "]  RES : chinese_remainder(" << a1 << "," << m1 << "," << a2 << "," << m2 << ") = " << res << endl;
#endif
 
    return res;
}
 
/*
    #] chinese_remainder : 
        #[ substitute :
*/
 
const poly polygcd::substitute(const poly &a, int x, int c) {
 
    POLY_GETIDENTITY(a);
 
    poly b(BHEAD 0);
 
    if (a.is_zero()) {
        return b;
    }
 
    bool zero=true;
    int bi=1;
 
    // cache size is bounded by the degree in x, twice the number of terms of
    // the polynomial and a constant
    vector<WORD> cache(min(a.degree(x)+1,min(2*a.number_of_terms(),
                                                                                     POLYGCD_RAISPOWMOD_CACHE_SIZE)), 0);
 
    for (int ai=1; ai<=a[0]; ai+=a[ai]) {
        // last term or different power, then add term to b iff non-zero
        if (!zero) {
            bool add=false;
            if (ai==a[0])
                add=true;
            else {
                for (int i=0; i<AN.poly_num_vars; i++)
                    if (i!=x && a[ai+1+i]!=b[bi+1+i]) {
                        zero=true;
                        add=true;
                        break;
                    }
            }
 
            if (add) {
                if (b[bi+AN.poly_num_vars+1] < 0)
                    b[bi+AN.poly_num_vars+1] += a.modp;
                bi+=b[bi];
            }
 
            if (ai==a[0]) break;
        }
        
        b.check_memory(bi);
 
        // create new term in b
        if (zero) {
            b[bi] = 3+AN.poly_num_vars;
            for (int i=0; i<AN.poly_num_vars; i++)
                b[bi+1+i] = a[ai+1+i];
            b[bi+1+x] = 0;
            b[bi+AN.poly_num_vars+1] = 0;
            b[bi+AN.poly_num_vars+2] = 1;
        }
 
        // add term of a to the current term in b
        LONG coeff = a[ai+1+AN.poly_num_vars] * a[ai+2+AN.poly_num_vars];
        int pow = a[ai+1+x];
        
        if (pow<(int)cache.size()) {
            if (cache[pow]==0) 
                cache[pow] = RaisPowMod(c, pow, a.modp);
            coeff = (coeff * cache[pow]) % a.modp;
        }
        else {
            coeff = (coeff * RaisPowMod(c, pow, a.modp)) % a.modp;
        }
        
        b[bi+AN.poly_num_vars+1] = (coeff + b[bi+AN.poly_num_vars+1]) % a.modp;
        if (b[bi+AN.poly_num_vars+1] != 0) zero=false;
    }
 
    b[0]=bi;
    b.setmod(a.modp);
 
    return b;   
}
 
/*
        #] substitute : 
        #[ sparse_interpolation helper functions :
*/
 
// Returns a list of size #terms(a) with entries PROD(ci^powi, i=2..n)
const vector<int> polygcd::sparse_interpolation_get_mul_list (const poly &a, const vector<int> &x, const vector<int> &c) {
    // cache size for variable x is bounded by the degree in x, twice
    // the number of terms of the polynomial and a constant
    vector<vector<WORD> > cache(c.size());
    int max_cache_size = min(2*a.number_of_terms(),POLYGCD_RAISPOWMOD_CACHE_SIZE);
    for (int i=0; i<(int)c.size(); i++)
        cache[i] = vector<WORD>(min(a.degree(x[i+1])+1,max_cache_size), 0);
    
    vector<int> res;
    for (int i=1; i<a[0]; i+=a[i]) {
        LONG coeff=1;
        for (int j=0; j<(int)c.size(); j++) {
            int pow = a[i+1+x[j+1]];
            if (pow<(int)cache[j].size()) {
                if (cache[j][pow]==0) 
                    cache[j][pow] = RaisPowMod(c[j], pow, a.modp);
                coeff = (coeff * cache[j][pow]) % a.modp;
            }
            else {
                coeff = (coeff * RaisPowMod(c[j], pow, a.modp)) % a.modp;
            }
        }
        res.push_back(coeff);
    }
    return res;
}
 
// Multiplies the coefficients of a with the entries of mul
void polygcd::sparse_interpolation_mul_poly (poly &a, const vector<int> &mul) {
    for (int i=1,j=0; i<a[0]; i+=a[i],j++) 
        a[i+a[i]-2] = ((LONG)a[i+a[i]-2]*mul[j]) % a.modp;
}
 
// Sets all coefficients to the range 0..modp-1 and the powers of x2...xn to 0
const poly polygcd::sparse_interpolation_reduce_poly (const poly &a, const vector<int> &x) {
    poly res(a);
    for (int i=1; i<a[0]; i+=a[i]) {
        for (int j=1; j<(int)x.size(); j++)
            res[i+1+x[j]]=0;
        if (res[i+a[i]-1]==-1) {
            res[i+a[i]-1]=1;
            res[i+a[i]-2]=a.modp-res[i+a[i]-2];
        }
    }
    return res;
}
 
// Collects entries with equal powers, so that the result is a proper polynomial
const poly polygcd::sparse_interpolation_fix_poly (const poly &a, int x) {
    
    POLY_GETIDENTITY(a);
    poly res(BHEAD 0,a.modp,1);
 
    int j=1;
    bool newterm=true;
        
    for (int i=1; i<a[0]; i+=a[i]) {
        if (newterm)
            res.termscopy(&a[i], j, a[i]);
        else 
            res[j+res[j]-2] = ((LONG)res[j+res[j]-2] + a[i+a[i]-2]) % a.modp;
        
        newterm = i+a[i] == a[0] || res[j+1+x] != a[i+a[i]+1+x];
        if (newterm && res[j+res[j]-2]!=0) j += res[j];
    }
 
    res[0]=j;
    return res;
}
 
/*
        #] sparse_interpolation helper functions : 
        #[ gcd_modular_sparse_interpolation :
*/
 
const poly polygcd::gcd_modular_sparse_interpolation (const poly &origa, const poly &origb, const vector<int> &x, const poly &s) {
 
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  CALL: gcd_modular_sparse_interpolation("
             << origa << "," << origb << "," << x << "," << "," << s <<")" << endl;
#endif
 
    POLY_GETIDENTITY(origa);
 
    // strip multivariate content
    poly conta(content_multivar(origa,x.back()));
    poly contb(content_multivar(origb,x.back()));
    poly gcdconts(gcd_Euclidean(conta,contb));
    const poly& a = conta.is_one() ? origa : origa/conta;
    const poly& b = contb.is_one() ? origb : origb/contb;
 
    // for non-monic cases, we need to normalize with the gcd of the lcoeffs of a poly in x[0]
    // or else the shape fitting does not work.
    // FIXME: the current implementation still rejects some valid shapes.
    poly lcgcd(BHEAD 1, a.modp);
    if (!s.lcoeff_univar(x[0]).is_integer()) {
        lcgcd = gcd_modular_dense_interpolation(a.lcoeff_univar(x[0]), b.lcoeff_univar(x[0]), x, poly(BHEAD 0));
    }
 
    // reduce polynomials
    poly ared(sparse_interpolation_reduce_poly(a,x));
    poly bred(sparse_interpolation_reduce_poly(b,x));
    poly sred(sparse_interpolation_reduce_poly(s,x));
    poly lred(sparse_interpolation_reduce_poly(lcgcd,x));
 
    // set all coefficients to 1
    for (int i=1; i<sred[0]; i+=sred[i]) {
        sred[i+sred[i]-2] = sred[i+sred[i]-1] = 1;
    }
    
    // generate random numbers and check there this set doesn't result
    // in a singular matrix
    vector<int> c(x.size()-1);
    vector<int> smul;
 
    bool duplicates;
    do {
        for (int i=0; i<(int)c.size(); i++)
            c[i] = 1 + wranf(BHEAD0) % (a.modp-1);
        smul = sparse_interpolation_get_mul_list(s,x,c);
 
        duplicates = false;
 
        int fr=0,to=0;
        for (int i=1; i<s[0];) {
            int pow = s[i+1+x[0]];
            while (i<s[0] && s[i+1+x[0]]==pow) i+=s[i], to++;
            for (int j=fr; j<to; j++)
                for (int k=fr; k<j; k++)
                    if (smul[j] == smul[k]) 
                        duplicates = true;
            fr=to;
        }       
    }
    while (duplicates);
 
    // get the lists to multiply the polynomials with every iteration
    vector<int> amul(sparse_interpolation_get_mul_list(a,x,c));
    vector<int> bmul(sparse_interpolation_get_mul_list(b,x,c));
    vector<int> lmul(sparse_interpolation_get_mul_list(lcgcd,x,c));
 
    vector<vector<vector<LONG> > > M;
    vector<vector<LONG> > V;
 
    int maxMsize=0;
    
    // create (empty) matrices
    for (int i=1; i<s[0]; i+=s[i]) {
        if (i==1 || s[i+1+x[0]]!=s[i+1+x[0]-s[i]]) {
            M.push_back(vector<vector<LONG> >());
            V.push_back(vector<LONG>());
        }
        M.back().push_back(vector<LONG>());
        V.back().push_back(0);
        maxMsize = max(maxMsize, (int)M.back().size());
    }
 
    // generate linear equations
    for (int numg=0; numg<maxMsize; numg++) {
 
        poly amodI(sparse_interpolation_fix_poly(ared,x[0]));
        poly bmodI(sparse_interpolation_fix_poly(bred,x[0]));
        poly lmodI(sparse_interpolation_fix_poly(lred,x[0]));
 
        // A fix for non-monic gcds. This could be slow if lmodI has many terms,
        // since it overfits the gcd now. Another gcd has to be run to remove the
        // extra terms.
        poly gcd(lmodI * gcd_Euclidean(amodI,bmodI));
 
        // if correct gcd
        if (!gcd.is_zero() && gcd[2+x[0]]==sred[2+x[0]]) {
 
            // for each power in the gcd, generate an equation if needed
            int gi=1, midx=0;
            
            for (int si=1; si<s[0];) {
                // if the term exists, set Vi=coeff, otherwise Vi remains 0
                if (gi<gcd[0] && gcd[gi+1+x[0]]==sred[si+1+x[0]]) {
                    if (numg < (int)V[midx].size()) 
                        V[midx][numg] = gcd[gi+gcd[gi]-1]*gcd[gi+gcd[gi]-2];
                    gi += gcd[gi];
                }
 
                // add the coefficients of s to the matrix M
                for (int i=0; i<(int)M[midx].size(); i++) {
                    if (numg < (int)M[midx].size())
                        M[midx][numg].push_back(sred[si+1+AN.poly_num_vars]);
                    si += s[si];
                }
                
                midx++;
            }
        }
        else {
            // incorrect gcd
            if (!gcd.is_zero() && gcd[2+x[0]]<sred[2+x[0]])
                return poly(BHEAD 0);
            numg--;
        }
        
        // multiply polynomials by the lists to obtain new ones
        sparse_interpolation_mul_poly(ared,amul);
        sparse_interpolation_mul_poly(bred,bmul);
        sparse_interpolation_mul_poly(sred,smul);
        sparse_interpolation_mul_poly(lred,lmul);
    }
 
    // solve the linear equations
    for (int i=0; i<(int)M.size(); i++) {
        int n = M[i].size();
 
        // Gaussian elimination
        for (int j=0; j<n; j++) {
            for (int k=0; k<j; k++) {
                LONG x = M[i][j][k];
                for (int l=k; l<n; l++) 
                    M[i][j][l] = (M[i][j][l] - M[i][k][l]*x) % a.modp;
                V[i][j] = (V[i][j] - V[i][k]*x) % a.modp;
            }
            
            // normalize row
            WORD x = M[i][j][j]; // WORD for GetModInverses
            GetModInverses(x + (x<0?a.modp:0), a.modp, &x, NULL);
            for (int k=0; k<n; k++) 
                M[i][j][k] = (M[i][j][k]*x) % a.modp;
            V[i][j] = (V[i][j]*x) % a.modp;
        }
 
        // solve
        for (int j=n-1; j>=0; j--)
            for (int k=j+1; k<n; k++) 
                V[i][j] = (V[i][j] - M[i][j][k]*V[i][k]) % a.modp;
    }
 
    // create coefficient list
    vector<LONG> coeff;
    for (int i=0; i<(int)V.size(); i++)
        for (int j=0; j<(int)V[i].size(); j++) 
            coeff.push_back(V[i][j]);
    
    // create resulting polynomial
    poly res(BHEAD 0);
    int ri=1, i=0;
    for (int si=1; si<s[0]; si+=s[si]) {
        res.check_memory(ri);
        res[ri] = 3 + AN.poly_num_vars;             // term length
        for (int j=0; j<AN.poly_num_vars; j++) 
            res[ri+1+j] = s[si+1+j];                  // powers
        res[ri+1+AN.poly_num_vars] = ABS(coeff[i]); // coefficient
        res[ri+2+AN.poly_num_vars] = SGN(coeff[i]); // coefficient length
        i++;
        ri += res[ri];
    }
    res[0]=ri;                                    // total length
    res.setmod(a.modp,1);
 
    if (!poly::divides(res, lcgcd * a) || !poly::divides(res, lcgcd * b)) {
        return poly(BHEAD 0); // bad shape
    } else {
        // refine gcd
        if (!poly::divides(res, a))
            res = gcd_modular_dense_interpolation(res, a, x, poly(BHEAD 0));
        if (!poly::divides(res, b))
            res = gcd_modular_dense_interpolation(res, b, x, poly(BHEAD 0));
    }
 
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  RES : gcd_modular_sparse_interpolation("
             << a << "," << b << "," << x << "," << "," << s <<") = " << res << endl;
#endif
    
    return gcdconts * res;
}
 
/*
    #] gcd_modular_sparse_interpolation : 
        #[ gcd_modular_dense_interpolation :
*/
 
const poly polygcd::gcd_modular_dense_interpolation (const poly &a, const poly &b, const vector<int> &x, const poly &s) {
    
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  CALL: gcd_modular_dense_interpolation(" << a << "," << b << "," << x << "," << s <<")" << endl;
#endif
 
    POLY_GETIDENTITY(a);
 
    // if univariate, then use Euclidean algorithm
    if (x.size() == 1) {
        return gcd_Euclidean(a,b);
    }
 
    // if shape is known, use sparse interpolation
    if (!s.is_zero()) {
        poly res = gcd_modular_sparse_interpolation (a,b,x,s);
        if (!res.is_zero()) return res;
        // apparently the shape was not correct. continue.
    }
 
    // divide out multivariate content in last variable
    int X = x.back();
 
    poly conta(content_multivar(a,X));
    poly contb(content_multivar(b,X));
    poly gcdconts(gcd_Euclidean(conta,contb));
    const poly& ppa = conta.is_one() ? a : poly(a/conta);
    const poly& ppb = contb.is_one() ? b : poly(b/contb);
 
    // gcd of leading coefficients
    poly lcoeffa(ppa.lcoeff_multivar(X));
    poly lcoeffb(ppb.lcoeff_multivar(X));
    poly gcdlcoeffs(gcd_Euclidean(lcoeffa,lcoeffb));
 
    // calculate the degree bound for each variable
    int m = MiN(ppa.degree(x[x.size() - 2]),ppb.degree(x[x.size() - 2]));
 
    poly res(BHEAD 0);
    poly oldres(BHEAD 0);
    poly newshape(BHEAD 0);
    poly modpoly(BHEAD 1,a.modp);
    
    while (true) {
        // generate random constants and substitute it
        int c = 1 + wranf(BHEAD0) % (a.modp-1);
        if (substitute(gcdlcoeffs,X,c).is_zero()) continue;
        if (substitute(modpoly,X,c).is_zero()) continue;
        
        poly amodc(substitute(ppa,X,c));
        poly bmodc(substitute(ppb,X,c));
 
        // calculate gcd recursively
        poly gcdmodc(gcd_modular_dense_interpolation(amodc,bmodc,vector<int>(x.begin(),x.end()-1), newshape));
        int n = gcdmodc.degree(x[x.size() - 2]);
 
        // normalize
        gcdmodc = (gcdmodc * substitute(gcdlcoeffs,X,c)) / gcdmodc.integer_lcoeff();
        poly simple(poly::simple_poly(BHEAD X,c,1,a.modp)); // (X-c) mod p
 
        // if power is smaller, the old one was wrong
        if ((res.is_zero() && n == m) || n < m) {
            m = n;
            res = gcdmodc;
            newshape = gcdmodc; // set a new shape (interpolation does not change it)
            modpoly = simple;
        }
        else if (n == m) {
            oldres = res;
            // equal powers, so interpolate results
            poly coeff_poly(substitute(modpoly,X,c));
            WORD coeff_word = coeff_poly[2+AN.poly_num_vars] * coeff_poly[3+AN.poly_num_vars];
            if (coeff_word < 0) coeff_word += a.modp;
 
            GetModInverses(coeff_word, a.modp, &coeff_word, NULL);
            
            res.setmod(a.modp); // make sure the mod is set before substituting
            res += poly(BHEAD coeff_word, a.modp, 1) * modpoly * (gcdmodc - substitute(res,X,c));
            modpoly *= simple;
        }
 
        // check whether this is the complete gcd
        if (!res.is_zero() && res == oldres) {
            poly nres = res / content_multivar(res, X);
            if (poly::divides(nres,ppa) && poly::divides(nres,ppb)) {
#ifdef DEBUG
                cout << "*** [" << thetime() << "]  RES : gcd_modular_dense_interpolation(" << a << "," << b << ","
                         << x << "," << "," << s <<") = " << gcdconts * nres << endl;
#endif
                return gcdconts * nres;
            }
 
            // At this point, the gcd may be too large due to bad luck
            // TODO: create an efficient fail state that tries to find a smaller
            // polynomial without interpolating bad ones?
            newshape = poly(BHEAD 0); // reset the shape, important!
        }
    }
}
 
/*
    #] gcd_modular_dense_interpolation : : 
        #[ gcd_modular :
*/
 
const poly polygcd::gcd_modular (const poly &origa, const poly &origb, const vector<int> &x) {
 
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  CALL: gcd_modular(" << origa << "," << origb << "," << x << ")" << endl;
#endif
 
    POLY_GETIDENTITY(origa);
 
    if (origa.is_zero()) return origa.modp==0 ? origb : origb / origb.integer_lcoeff();
    if (origb.is_zero()) return origa.modp==0 ? origa : origa / origa.integer_lcoeff();
    if (origa==origb) return origa.modp==0 ? origa : origa / origa.integer_lcoeff();
 
    poly ac = integer_content(origa);
    poly bc = integer_content(origb);
    const poly& a = ac.is_one() ? origa : poly(origa/ac);
    const poly& b = bc.is_one() ? origb : poly(origb/bc);
    poly ic = integer_gcd(ac, bc);
    poly g = integer_gcd(a.integer_lcoeff(), b.integer_lcoeff());
 
    int pnum=0;
    
    poly d(BHEAD 0);
    poly m1(BHEAD 1);
    int mindeg=MAXPOSITIVE;
 
    while (true) {
        // choose a prime and solve modulo the prime
        WORD p = NextPrime(BHEAD pnum++);
        if (poly(a.integer_lcoeff(),p).is_zero()) continue;
        if (poly(b.integer_lcoeff(),p).is_zero()) continue;
 
        poly c(gcd_modular_dense_interpolation(poly(a,p),poly(b,p),x,poly(d,p)));
        c = (c * poly(g,p)) / c.integer_lcoeff(); // normalize so that lcoeff(c) = g mod p
 
        if (c.is_zero()) {
            // unlucky choices somewhere, so start all over again
            d = poly(BHEAD 0);
            m1 = poly(BHEAD 1);
            mindeg = MAXPOSITIVE;
            continue;
        }
        
        if (!(poly(a,p)%c).is_zero()) continue;
        if (!(poly(b,p)%c).is_zero()) continue;
 
        int deg = c.degree(x[0]);
 
        if (deg < mindeg) {
            // small degree, so the old one is wrong
            d=c;
            d.modp=a.modp;
            d.modn=a.modn;
            m1 = poly(BHEAD p);
            mindeg=deg;
        }
        else if (deg == mindeg) {
            // same degree, so use Chinese Remainder Algorithm
            poly newd(BHEAD 0);
            
            for (int ci=1,di=1; ci<c[0]||di<d[0]; ) {
                int comp = ci==c[0] ? -1 : di==d[0] ? +1 : poly::monomial_compare(BHEAD &c[ci],&d[di]);
                poly a1(BHEAD 0),a2(BHEAD 0);
                
                newd.check_memory(newd[0]);
                
                if (comp <= 0) {
                    newd.termscopy(&d[di],newd[0],1+AN.poly_num_vars);
                    a1 = poly(BHEAD (UWORD *)&d[di+1+AN.poly_num_vars],d[di+d[di]-1]);
                    di+=d[di];
                }
                if (comp >= 0) {
                    newd.termscopy(&c[ci],newd[0],1+AN.poly_num_vars);
                    a2 = poly(BHEAD (UWORD *)&c[ci+1+AN.poly_num_vars],c[ci+c[ci]-1]);
                    ci+=c[ci];
                }
 
                poly e(chinese_remainder(a1,m1,a2,poly(BHEAD p)));
                newd.termscopy(&e[2+AN.poly_num_vars], newd[0]+1+AN.poly_num_vars, ABS(e[e[1]])+1);
                newd[newd[0]] = 2 + AN.poly_num_vars + ABS(e[e[1]]);
                newd[0] += newd[newd[0]];
            }
 
            m1 *= poly(BHEAD p);
            d=newd;
        }
 
        // divide out spurious integer content
        poly ppd(d / integer_content(d));
 
        // check whether this is the complete gcd
        if (poly::divides(ppd,a) && poly::divides(ppd,b)) {
#ifdef DEBUG
            cout << "*** [" << thetime() << "]  RES : gcd_modular(" << origa << "," << origb << "," << x << ") = "
                     << ic * ppd << endl;
#endif
            return ic * ppd;
        }
#ifdef DEBUG
        cout << "*** [" << thetime() << "] Retrying modular_gcd with new prime" << endl;
#endif
    }
}
 
/*
    #] gcd_modular : 
    #[ gcd_heuristic_possible :
*/
 
bool gcd_heuristic_possible (const poly &a) {
 
    POLY_GETIDENTITY(a);
    
    double prod_deg = 1;
    for (int j=0; j<AN.poly_num_vars; j++)
        prod_deg *= a[2+j]+1;
 
    double digits = ABS(a[1+a[1]-1]);
  double lead = a[1+1+AN.poly_num_vars];
        
    return prod_deg*(digits-1+log(2*ABS(lead))/log(2.0)/(BITSINWORD/2)) < POLYGCD_HEURISTIC_MAX_DIGITS;
}
 
/*
    #] gcd_heuristic_possible : 
    #[ gcd_heuristic :
*/
 
 
const poly polygcd::gcd_heuristic (const poly &a, const poly &b, const vector<int> &x, int max_tries) {
 
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  CALL: gcd_heuristic("<<a<<","<<b<<","<<x<<")\n";
#endif
 
    if (a.is_integer()) return integer_gcd(a,integer_content(b));
    if (b.is_integer()) return integer_gcd(integer_content(a),b);
 
    POLY_GETIDENTITY(a);
 
    // Calculate xi = 2*min(max(coefficients a),max(coefficients b))+2
 
    UWORD *pxi=NULL;
    WORD nxi=0;
 
    for (int i=1; i<a[0]; i+=a[i]) {
        WORD na = ABS(a[i+a[i]-1]);
        if (BigLong((UWORD *)&a[i+a[i]-1-na], na, pxi, nxi) > 0) {
            pxi = (UWORD *)&a[i+a[i]-1-na];
            nxi = na;
        }
    }
    
    for (int i=1; i<b[0]; i+=b[i]) {        
        WORD nb = ABS(b[i+b[i]-1]);
        if (BigLong((UWORD *)&b[i+b[i]-1-nb], nb, pxi, nxi) > 0) {
            pxi = (UWORD *)&b[i+b[i]-1-nb];
            nxi = nb;
        }
    }
 
    poly xi(BHEAD pxi,nxi);
    
    // Addition of another random factor gives better performance
    xi = xi*poly(BHEAD 2) + poly(BHEAD 2 + wranf(BHEAD0)%POLYGCD_HEURISTIC_MAX_ADD_RANDOM);
 
    // If degree*digits(xi) is too large, throw exception
    if (max(a.degree(x[0]),b.degree(x[0])) * xi[xi[1]] >= MiN(AM.MaxTal, POLYGCD_HEURISTIC_MAX_DIGITS)) {
#ifdef DEBUG
        cout << "*** [" << thetime() << "]  RES : gcd_heuristic("<<a<<","<<b<<","<<x<<") = overflow\n";
#endif
        throw(gcd_heuristic_failed());
    }
 
    for (int times=0; times<max_tries; times++) {
 
        // Recursively calculate the gcd
 
        poly gamma(gcd_heuristic(a % poly::simple_poly(BHEAD x[0],xi),
                                                         b % poly::simple_poly(BHEAD x[0],xi),
                                                         vector<int>(x.begin()+1,x.end()),1));
                                                                    
        // If a gcd is found, reconstruct the powers of x
        if (!gamma.is_zero()) {
            // res is construct is reverse order. idx/len are for reversing
            // it in the correct order
            poly res(BHEAD 0), c(BHEAD 0);
            vector<int> idx, len;
            
            for (int power=0; !gamma.is_zero(); power++) {
 
                // calculate c = gamma % xi (c and gamma are polynomials, xi is integer)
                c = gamma;
                c.coefficients_modulo((UWORD *)&xi[2+AN.poly_num_vars], xi[xi[0]-1], false);
                
                // Add the terms c * x^power to res
                res.check_memory(res[0]+c[0]);
                res.termscopy(&c[1],res[0],c[0]-1);
                for (int i=1; i<c[0]; i+=c[i])
                    res[res[0]-1+i+1+x[0]] = power;
 
                // Store idx/len for reversing
                if (!c.is_zero()) {
                    idx.push_back(res[0]);
                    len.push_back(c[0]-1);
                }
                
                res[0] += c[0]-1;
 
                // Divide gamma by xi
                gamma = (gamma - c) / xi;
            }
            
            // Reverse the resulting polynomial
            poly rev(BHEAD 0);
            rev.check_memory(res[0]);
            
            rev[0] = 1;
            for (int i=idx.size()-1; i>=0; i--) {
                rev.termscopy(&res[idx[i]], rev[0], len[i]);
                rev[0] += len[i];
            }
            res = rev;
 
            poly ppres = res / integer_content(res);
 
            if ((a%ppres).is_zero() && (b%ppres).is_zero()) {
#ifdef DEBUG
                cout << "*** [" << thetime() << "]  RES : gcd_heuristic("<<a<<","<<b<<","<<x<<") = "<<res<<"\n";
#endif
                return res;
            }
        }
 
        // Next xi by multiplying with the golden ratio to avoid correlated errors
        xi = xi * poly(BHEAD 28657) / poly(BHEAD 17711) + poly(BHEAD wranf(BHEAD0) % POLYGCD_HEURISTIC_MAX_ADD_RANDOM);
    }
    
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  RES : gcd_heuristic("<<a<<","<<b<<","<<x<<") = failed\n";
#endif
 
    return poly(BHEAD 0);
}
 
/*
    #] gcd_heuristic : 
    #[ bracket :
*/
 
const map<vector<int>,poly> polygcd::full_bracket(const poly &a, const vector<int>& filter) {
    POLY_GETIDENTITY(a);
 
    map<vector<int>,poly> bracket;
    for (int ai=1; ai<a[0]; ai+=a[ai]) {
        vector<int> varpattern(AN.poly_num_vars);
        for (int i=0; i<AN.poly_num_vars; i++)
            if (filter[i] == 1 && a[ai + i + 1] > 0)
                varpattern[i] = a[ai + i + 1];
 
        // create monomial
        poly mon(BHEAD 1);
        mon.setmod(a.modp);
        mon[0] = a[ai] + 1;
        for (int i=0; i<a[ai]; i++)
            if (i > 0 && i <= AN.poly_num_vars && varpattern[i - 1])
                mon[1+i] = 0;
            else
                mon[1+i] = a[ai+i];
 
        map<vector<int>,poly>::iterator i = bracket.find(varpattern);
        if (i == bracket.end()) {
            bracket.insert(std::make_pair(varpattern, mon));
        } else {
            i->second += mon;
        }
    }
 
    return bracket;
}
 
const poly polygcd::bracket(const poly &a, const vector<int>& pattern, const vector<int>& filter) {
    POLY_GETIDENTITY(a);
 
    poly bracket(BHEAD 0);
    for (int ai=1; ai<a[0]; ai+=a[ai]) {
        bool ispat = true;
        for (int i=0; i<AN.poly_num_vars; i++)
            if (filter[i] == 1 && pattern[i] != a[ai + i + 1]) {
                ispat = false;
                break;
            }
 
        if (ispat) {
            poly mon(BHEAD 1);
            mon.setmod(a.modp);
            mon[0] = a[ai] + 1;
            for (int i=0; i<a[ai]; i++)
                if (i > 0 && i <= AN.poly_num_vars && pattern[i - 1])
                    mon[1+i] = 0;
                else
                    mon[1+i] = a[ai+i];
            bracket += mon;
        }
    }
 
    return bracket;
}
 
const map<vector<int>,int> polygcd::bracket_count(const poly &a, const vector<int>& filter) {
    POLY_GETIDENTITY(a);
 
    map<vector<int>,int> bracket;
    for (int ai=1; ai<a[0]; ai+=a[ai]) {
        vector<int> varpattern(AN.poly_num_vars);
        for (int i=0; i<AN.poly_num_vars; i++)
            if (filter[i] == 1 && a[ai + i + 1] > 0)
                varpattern[i] = a[ai + i + 1];
 
        map<vector<int>,int>::iterator i = bracket.find(varpattern);
        if (i == bracket.end()) {
            bracket.insert(std::make_pair(varpattern, 0));
        } else {
            i->second++;
        }
    }
 
    return bracket;
}
 
struct BracketInfo {
    std::vector<int> pattern;
    int num_terms, dummy = 0;
    const poly* p;
 
    BracketInfo(const std::vector<int>& pattern, int num_terms, const poly* p) : pattern(pattern), num_terms(num_terms), p(p) {}
    bool operator<(const BracketInfo& rhs) const { return num_terms > rhs.num_terms; } // biggest should be first!
};
 
/*
    #] bracket : 
    #[ gcd_linear:
*/
 
const poly gcd_linear_helper (const poly &a, const poly &b) {
    POLY_GETIDENTITY(a);
 
    for (int i = 0; i < AN.poly_num_vars; i++)
        if (a.degree(i) == 1) {
            vector<int> filter(AN.poly_num_vars);
            filter[i] = 1;
 
            // bracket the linear variable
            map<vector<int>,poly> ba = polygcd::full_bracket(a, filter);
 
            poly subgcd(BHEAD 1);
            if (ba.size() == 2)
                subgcd = polygcd::gcd_linear(ba.begin()->second, (++ba.begin())->second);
            else
                subgcd = ba.begin()->second;
 
            poly linfac = a / subgcd;
            if (poly::divides(linfac,b))
                return linfac * polygcd::gcd_linear(subgcd, b / linfac);
 
            return polygcd::gcd_linear(subgcd, b);
        }
 
    return poly(BHEAD 0);
}
 
const poly polygcd::gcd_linear (const poly &a, const poly &b) {
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  CALL: gcd_linear("<<a<<","<<b<<")\n";
#endif
 
    POLY_GETIDENTITY(a);
 
    if (a.is_zero()) return a.modp==0 ? b : b / b.integer_lcoeff();
    if (b.is_zero()) return a.modp==0 ? a : a / a.integer_lcoeff();
    if (a==b) return a.modp==0 ? a : a / a.integer_lcoeff();
 
    if (a.is_integer() || b.is_integer()) {
        if (a.modp > 0) return poly(BHEAD 1,a.modp,a.modn);
        return poly(integer_gcd(integer_content(a),integer_content(b)),0,1);
    }
 
    poly h = gcd_linear_helper(a, b);
    if (!h.is_zero()) return h;
 
    h = gcd_linear_helper(b, a);
    if (!h.is_zero()) return h;
 
    vector<int> xa = a.all_variables();
    vector<int> xb = b.all_variables();
 
    vector<int> used(AN.poly_num_vars,0);
    for (int i=0; i<(int)xa.size(); i++) used[xa[i]]++;
    for (int i=0; i<(int)xb.size(); i++) used[xb[i]]++;
    vector<int> x;
    for (int i=0; i<AN.poly_num_vars; i++)
        if (used[i]) x.push_back(i);
 
    return gcd_modular(a,b,x);
}
 
/*
    #] gcd_linear: 
    #[ gcd :
*/
 
const poly polygcd::gcd (const poly &a, const poly &b) {
 
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  CALL: gcd("<<a<<","<<b<<")\n";
#endif
 
    POLY_GETIDENTITY(a);
    
    if (a.is_zero()) return a.modp==0 ? b : b / b.integer_lcoeff();
    if (b.is_zero()) return a.modp==0 ? a : a / a.integer_lcoeff();
    if (a==b) return a.modp==0 ? a : a / a.integer_lcoeff();
 
    if (a.is_integer() || b.is_integer()) {
        if (a.modp > 0) return poly(BHEAD 1,a.modp,a.modn);
        return poly(integer_gcd(integer_content(a),integer_content(b)),0,1);
    }
 
    // Generate a list of variables of a and b
    vector<int> xa = a.all_variables();
    vector<int> xb = b.all_variables();
 
    vector<int> used(AN.poly_num_vars,0);
    for (int i=0; i<(int)xa.size(); i++) used[xa[i]]++;
    for (int i=0; i<(int)xb.size(); i++) used[xb[i]]++;
    vector<int> x;
    for (int i=0; i<AN.poly_num_vars; i++)
        if (used[i]) x.push_back(i);
 
    if (a.is_monomial() || b.is_monomial()) {
 
        poly res(BHEAD 1,a.modp,a.modn);
        if (a.modp == 0) res = integer_gcd(integer_content(a),integer_content(b));
        
        for (int i=0; i<(int)x.size(); i++)
            res[2+x[i]] = 1<<(BITSINWORD-2);
 
        for (int i=1; i<a[0]; i+=a[i]) 
            for (int j=0; j<(int)x.size(); j++) 
                res[2+x[j]] = MiN(res[2+x[j]], a[i+1+x[j]]);
        
        for (int i=1; i<b[0]; i+=b[i]) 
            for (int j=0; j<(int)x.size(); j++) 
                res[2+x[j]] = MiN(res[2+x[j]], b[i+1+x[j]]);
 
        return res;
    }
 
    // Calculate the contents, their gcd and the primitive parts
    poly conta(x.size()==1 ? integer_content(a) : content_univar(a,x[0]));
    poly contb(x.size()==1 ? integer_content(b) : content_univar(b,x[0]));
    poly gcdconts(x.size()==1 ? integer_gcd(conta,contb) : gcd(conta,contb));
    const poly& ppa = conta.is_one() ? a : poly(a/conta);
    const poly& ppb = contb.is_one() ? b : poly(b/contb);
    
    if (ppa == ppb) 
        return ppa * gcdconts;
 
    poly res(BHEAD 0);
    
#ifdef POLYGCD_USE_HEURISTIC
    // Try the heuristic gcd algorithm
    if (a.modp==0 && gcd_heuristic_possible(a) && gcd_heuristic_possible(b)) {
        try {
            res = gcd_heuristic(ppa,ppb,x);
            if (!res.is_zero()) res /= integer_content(res);
        }
        catch (gcd_heuristic_failed) {}
    }
#endif
    
    // If gcd==0, the heuristic algorithm failed, so we do more extensive checks.
    // First, we filter out variables that appear in only one of the expressions.
    if (res.is_zero()) {
        bool unusedVars = false;
        for (unsigned int i = 0; i < used.size(); i++) {
            if (used[i] == 1) {
                unusedVars = true;
                break;
            }
        }
 
        // if there are no unused variables, go to the linear routine directly
        if (!unusedVars) {
            res = gcd_linear(ppa,ppb);
#ifdef DEBUG
            cout << "New GCD attempt (unused vars): " << res << endl;
#endif
        }
 
        // if res is not the gcd, it is 0 or larger than the gcd.
        // we bracket the expression in all the variables that appear only in one expr.
        // and we refine the gcd.
        bool diva = !res.is_zero() && poly::divides(res,ppa);
        bool divb = !res.is_zero() && poly::divides(res,ppb);
        if (!diva || !divb) {
            vector<BracketInfo> bracketinfo;
 
            if (!diva) {
                map<vector<int>,int> ba = bracket_count(ppa, used);
                for(map<vector<int>,int>::iterator it = ba.begin(); it != ba.end(); it++)
                    bracketinfo.push_back(BracketInfo(it->first, it->second, &ppa));
            }
 
            if (!divb) {
                map<vector<int>,int> bb = bracket_count(ppb, used);
                for(map<vector<int>,int>::iterator it = bb.begin(); it != bb.end(); it++)
                    bracketinfo.push_back(BracketInfo(it->first, it->second, &ppb));
            }
 
            // sort so that the smallest bracket will be last
            sort(bracketinfo.begin(), bracketinfo.end());
 
            if (res.is_zero()) {
                res = bracket(*bracketinfo.back().p, bracketinfo.back().pattern, used);
                bracketinfo.pop_back();
            }
 
            while (bracketinfo.size() > 0) {
                poly subpoly(bracket(*bracketinfo.back().p, bracketinfo.back().pattern, used));
                if (!poly::divides(res,subpoly)) {
                    // if we can filter out more variables, call gcd again
                    if (res.all_variables() != subpoly.all_variables())
                        res = gcd(subpoly,res);
                    else
                        res = gcd_linear(subpoly,res);
                }
 
                bracketinfo.pop_back();
            }
        }
 
        if (res.is_zero() || !poly::divides(res,ppa) || !poly::divides(res,ppb)) {
            MesPrint("Bad gcd found.");
            std::cout << "Bad gcd:" << res << " for " << ppa << " " << ppb << std::endl;
            Terminate(1);
        }
    }
 
    res *= gcdconts * poly(BHEAD res.sign());
 
#ifdef DEBUG
    cout << "*** [" << thetime() << "]  RES : gcd("<<a<<","<<b<<") = "<<res<<endl;
#endif
 
    return res;
}
 
/*
    #] gcd : 
*/