[模板][数论] 扩展欧几里得算法(ExGCD)
Templates 数论 数论 ExGCD
Lastmod: 2021-09-05 周日 19:38:49

扩展欧几里得算法(ExGCD)

// ax+by=gcd(a,b)=g
// check a>0 b>0
// if a != b, |x|<b, |y|<a
template<typename T>
void exgcd(T a,T b,T &g,T &x,T &y) {
    if(!b) { g=a; x=1; y=0; }
    else { exgcd(b,a%b,g,y,x); y-=a/b*x; }
}

求逆元

template<typename T>
void exgcd(T a,T b,T &g,T &x,T &y) {
    if(!b) { g=a; x=1; y=0; }
    else { exgcd(b,a%b,g,y,x); y-=a/b*x; }
}
template<typename T>
T inv(T a,T m) {
    T g,x,y; exgcd(a,m,g,x,y);
    if(g!=1) return -1; // no inverse element
    if(x<0) x+=m;
    return x;
}

不定方程

template<typename T, T INF=T(0x3f3f3f3f)>
struct Diophantine
{
    typedef pair<T,int> PTI;
    // floor
    T fl(T a,T b) { return (a^b)>0 ? a/b : a/b - (a%b!=0); }
    // ceiling
    T cl(T a,T b) { return (a^b)>0 ? a/b + (a%b!=0) : a/b; }
    void exgcd(T a,T b,T &g,T &x,T &y) {
        if(!b) { g=a; x=1; y=0; }
        else { exgcd(b,a%b,g,y,x); y-=a/b*x; }
    }
    // solve ax+by=m 
    // |a|x0 + |b|y0 = (a,b) = g
    // |a| * m/g * x0 + |b| * m/g * y0 = m
    // |a|*(m*x0/g + i * |b|/g) + |b|*(m*y0/g - i * |a|/g) = m
    // if the solution x,y exist return true
    // general solution is xi=x0 + i*dx_sgn*b, yi=y0-i*dy_sgn*a
    T a,b,g,x,y;
    int dx_sgn, dy_sgn;
    bool solve(T a_,T b_,T m) {
        a=a_; b=b_;
        if(a==0&&b==0) { x=y=0; return m==0; }
        if(a==0) { x=0; y=m/b; return m%b==0; }
        if(b==0) { y=0; x=m/a; return m%a==0; }
        exgcd(abs(a),abs(b),g,x,y);
        // print();
        if(m%g) return false;
        x*=m/g; y*=m/g;
        b/=g; a/=g;
        if(a<0) { a=-a; x=-x; dx_sgn=-1; } else { dx_sgn=1; }
        if(b<0) { b=-b; y=-y; dy_sgn=-1; } else { dy_sgn=1; }
        return true;
    }
    // a!=0 b!=0 
    // op = -1 : xlim <= x + i*sgn*b
    // op =  1 : xlim >= x + i*sgn*b
    // return {ilim, op} 
    // op = -1, i <= ilim
    // op =  1, i >= ilim 
    PTI x_to_i(T xlim,int op) {
        if(dx_sgn==1) {
            if(op==-1) return {cl(xlim-x,b),1}; // i>=(xlim-x)/b
            else return {fl(xlim-x,b),-1};      // i<=(xlim-x)/b
        } else { //sgn==-1
            if(op==-1) return {fl(x-xlim,b),-1}; // i<=(x-xlim)/b
            else return {cl(x-xlim,b),1};        // i>=(x-xlim)/b 
        }
    }
    // a!=0 b!=0 
    // op = -1 : ylim <= y - i*sgn*a
    // op =  1 : ylim >= y - i*sgn*a
    // return {ilim, op} 
    // op = -1, i <= ilim
    // op =  1, i >= ilim 
    PTI y_to_i(T ylim,int op) {
        if(dy_sgn==1) {
            if(op==-1) return {fl(y-ylim,a),-1}; // i<=(y-ylim)/a
            else return {cl(y-ylim,a),1};        // i>=(y-ylim)/a 
        } else {//sgn==-1
            if(op==-1) return {cl(ylim-y,a),1}; // i>=(ylim-y)/a
            else return {fl(ylim-y,a),-1};      // i<=(ylim-y)/a
        }
    }
    void print() {
        cout<<"a="<<a<<", b="<<b<<", g="<<g<<", x="<<x<<", y="<<y<<endl; }
};
Prev: [算法][数论] 扩展欧几里得算法
Next: [算法][图论] 约翰逊 Johnson 算法 全源最短路