欧几里得算法和扩展欧几里得算法

欧几里得算法和扩展欧几里得算法

扩欧几里得算法：
对所有正整数a,b, b != 0, gcd(a,b) = gcd(b, a mod b);
简略证明：
若a < b; 则 a ( mod b) = a, gcd(a,b) = gcd(b,a) 显然成立。
若a >= b, 假设a = q * b + r, 则 r = a mod b; 对于a,b的任意公约数d, 由于
d | a, d | q * b ,所以 d | (a - q * b) ,即d | r,因此d也是b,r的公约数，反之也成立。所以a,b公约数的集合和b,r(a mod b)的公约数集合相同，那么最大公约数也相同。
code实现:

typedef long long ll;
ll gcd(ll a,ll b){
    return b == 0 ? a : gcd(b,a % b);
}
// T(n) = O ( log(a + b))

扩展欧几里得算法：
引入定理：
对于任意整数a,b，存在一对整数x,y，满足 ax + by = gcd(a,b);
简略证明:
在扩展欧几里得算法的最后一步，即b == 0 时候，显然存在一对整数x = 1, y = 0,使得a * 1 + 0 * 0 = gcd(a,0);
若b > 0, 则gcd(a,b) = gcd(b, a mod b)。假设存在一对整数x, y满足 b * x + (a mod b) * y = gcd(b, a mod b) 因为 a mod b = a - [a / b] * b，所以式子可以变为
bx + (a - [a / b] * b) y = ay + (x - [a / b] y) b 这时候令 x ’ = y, y ’ = ( x - [a / b] y)
就有ax’ + by’ = gcd(a,b).
再对欧几里得算法的递归过程运用数学归纳法，可知定理成立。
code：

typedef long long ll;
ll exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){ x = 1, y = 0; return a;}
    ll d = exgcd(b,a % b, x, y);
    ll z = x; x = y; y = z - (a / b) * y;
    return d;
}
//定义变量d,x,y,调用d = exgcd(a,b,x,y) 注意x, y是引用传入。
//上述程序求出的是方程ax + by = gcd(a,b)的一组特解x,y.并返回了a,b的最大公约数
//对于更为一般的方程 ax + by = c,显然它有解的充分必要条件是d | c，d = gcd(a,b)
//所以我们可以先求出ax + by = d 的一组特解x0,y0,再同时乘上c / d, 就得到了
//ax + by = c 的一组特解 (c/d)x0, (c/d)y0;
/*
	更进一步方程ax + by = c的通解可表示为：
	x = c / d * x0 + k * (b / d)  //1)
	y = c / d * y0 - k * (a / d)  //2)
	其中k是任意整数，d = gcd(a,b),x0,y0是 ax + by = d的一组特解
	简略证明是将1) 2)代入ax + by.发现 k * ab / d相互抵消
*/

扩展欧几里得算法的应用：
1） 解线性同余方程：
给定整数a,b,m，求一个整数x 使得 a * x ≡ b (mod m), 无解则给出无解。
a * x ≡ b (mod m) 等价于a * x - b 是 m的倍数，不妨设为 -y 倍，于是该方程可以变为
a * x + m * y = b.
由前面可知该线性同余方程有解当且仅当 gcd(a,m) | b.

code：

ll exgcd(ll a, ll b,ll &x,ll &y){
    if(!b){
        return x = 1, y = 0, a;
    }
    ll d = exgcd(b , a % b, x, y);
    ll t = x;
    x = y, y = t - (a / b) * y;
    return d;
}
//solve ax + by = c return x; // x > 0
ll solve(ll a,ll b,ll c){
    ll x,y;
    ll d = exgcd(a,b,x,y);
    if(c % d) return -1; //无解 c 不能整除gcd(a,b)
    x *= c / d ; // x = c / d * x0 + k * ( b / d)
    ll tmp = b / d;
    return (x % tmp + tmp) % tmp;  //求解正整数解x
}

2） 特殊的线性同余方程 a * x ≡ 1 (mod m) (等价于求解乘法逆元)

code：

ll exgcd(ll a, ll b,ll &x,ll &y){
    if(!b){
        return x = 1, y = 0, a;
    }
    ll d = exgcd(b , a % b, x, y);
    ll t = x;
    x = y, y = t - (a / b) * y;
    return d;
}

ll get_inv(ll a,ll mod){
    ll x,y;
    exgcd(a,mod,x,y);
    return (x % mod + mod)%mod; //求解最小正逆元
}

3 ) 扩展求解线性同余方程组
先考虑两个方程组的情况 t ≡ r1 ( mod a1),(1) t ≡ r2( mod a2) (2)
(a1,a2不一定互质) 要求t(t > 0)
则有 t = a1 * x1 + r1, t = a2*x2 + r2；
变形为：a1 * x1 + a2 * (-x2) = ( r2 - r1)；
方程组有解等价于： (r2 - r1) | gcd(a1, a2);
我们可用扩展欧几里得算法求出 特解x10, x20; 求出最小正整数解x10；
那么要求的 t0 = a1 * x10 + r1;
显然 t0 满足1) 式和2)式
由此推论：t0 + k * lcm(a1,a2) 是(1)，(2)两个方程组的通解，lcm(a1,a2)都可消去。
在此通解基础上可以得到新的线性同余方程 t ≡ t0 (mod lcm(a1,a2)) （3）

两个以上的时候就等价于在先求出两个线性同余方程的通解t，在和第三个方程联立又求解一次t’，不断合并。

来道例题加深理解吧：

POJ - 2891
题意：给定n 个线性同余方程组 t ≡ ri (mod ai) ,求出满足题意的最小的正整数t，无解给出-1.

code：

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long ll;
template <typename T>
void read(T &x){
    x = 0;
    char c = getchar();
    for(; c < '0' || c > '9'; c = getchar()) ;
    for(; c >='0' && c <= '9'; c = getchar())
        x = (x << 1) + (x << 3) + (c ^ 48);
}


ll exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){ x = 1, y = 0; return a;}
    ll d = exgcd(b,a % b, x, y);
    ll z = x; x = y; y = z - (a / b) * y;
    return d;
}

int main(){
    ll n;
    while(cin>>n){
        ll a1,a2,r1,r2;
        bool ok = true;
        read(a1),read(r1);
        for(int i = 2; i <= n; ++i){
            read(a2),read(r2);
            if(!ok) continue ;
            ll x,y,d;
            d = exgcd(a1,a2,x,y); //solve a1 x + a2 y = d; (d = gcd(a1,a2))
            if((r2-r1)%d) ok = false;
            else{
                x *= (r2 - r1) / d ;
                ll tmp = a2 / d;
                x = (x % tmp + tmp) % tmp;
                r1 = r1 + x * a1; // r1 就是 t0
                a1 = a1 / d * a2; // lcm(a1,a2);
                //新的同余方程 t  ≡ r1 (mod a1)
            }
        }
        printf("%lld\n",ok ? r1 :-1);
    }
    return 0;
}