扩展欧几里得算法(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; }
};
Next: [算法][图论] 约翰逊 Johnson 算法 全源最短路