拓展欧几里得算法求不定方程

对于
　　　　　　　　ax+by=gcd(a,b) 

这样的方程，可以用扩展欧几里得算法exgcd求出一组通解。

根据欧几里得求gcd：
　　　　　　　　gcd(a,b)=gcd(b,a%b) 

可得
　　　　　　　　bx+(a%b)y=gcd(b,a%b) 

根据
　　　　　　a%b=a−(a/b)∗b 

可得
　　　　　　　　bx+ay−(a/b)b∗y=gcd(b,a%b)  

化简得
　　　　　　　　ay+b(x−(a/b)y)=gcd(b,a%b)  

　　　　　　　x′=y,y′=(x−(a/b)y)

　　　　ax′+by′=gcd(b,a%b)<=>ax+by=gcd(a,b)

根据

　　　　　　gcd(a,0)=a

一直递归直到b为0时可得
　　　　　　ax+by=a 

可以得出一组平凡解
　　　　　　　　　　x=1,y=0 

所以一直递归下去可以得出一组平凡解，然后再往回带得出ax+by=gcd(a,b)的一组解（x′=y,y′=(x−(a/b)y)） 

求解同余方程可以用费马小定理来求也可以用拓展欧几里得来求
ax≡b mod n <==>ax+ny=b

就转变成了上述形式

费马小定理：a是上能被质数p整除的正整数，则有a^(p-1) ≡ 1(mod p)

推导：a^(p-1) = 1(mod)p = a*a^(p-2)≡1 (mod p) 因此a的逆元为 a^(p-2); 所以对于满足费马小定理的可以直接用快速幂来求

例题 杭电  B - A/B

要求(A/B)%9973，但由于A很大，我们只给出n(n=A%9973)(我们给定的A必能被B整除，且gcd(B,9973) = 1)。

Input数据的第一行是一个T，表示有T组数据。 
每组数据有两个数n(0 <= n < 9973)和B(1 <= B <= 10^9)。Output对应每组数据输出(A/B)%9973。Sample Input

2
1000 53
87 123456789

Sample Output

7922
6060

除法逆元，，b*c=1(mod p) 则c就是b的逆元 那么a/b mod p==a/b*b*c(modp)==a*c(mod p)因此直接求出b的逆元c然后再乘以a取模就可以了
两种做法 ： 
1 费马小定理

#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const int MOD=9973;
ll ksm(ll a,ll b){
    ll res=1;
    while(b){
        if(b&1) res=res*a%MOD;
        a=a*a%MOD;
        b>>=1;
    }
    return res;
}
int main(){
    ll t;
    cin>>t;
    while(t--){
        ll a,b;
        cin>>a>>b;
        cout<<(a%MOD*ksm(b,MOD-2)%MOD)%MOD<<endl;
        
    }
    
    return 0;
}

2 拓展欧几里得

//a*x≡1(mod m) 等价于a*x+m*y=1(解释：对该公式俩边同时对m取摸，m*y对m取摸一定为0),因此求得x就是a的逆元 
// 可以用扩展欧几里得求
//得一组解，（x+m)mod m即为a的逆元
//该深度的x等于其下一个深度的Y 该深度的y等于其下一个深度的x-(a/b)*y
#include<iostream>
using namespace std;
typedef long long ll;
ll x,y;
const int MOD=9973;
int exgcd(ll a,ll b){
//--------该阶段可以就理解为求a与b的最大公约数,当b等于0时说明最大公约数找到了，就是此刻的a,然后在开始回溯更新x和y的值 
    if(b==0) {
        x=1;
        y=0;
        return a; 
    }
    int r=exgcd(b,a%b);
//---------------------     
    ll t=y;
    y=x-(a/b)*y;
    x=t;
    return r;
}

int main(){
    ll k;
    cin>>k;
    while(k--){
        ll a,b;
        cin>>a>>b;
        exgcd(b,MOD);
        x=(x+MOD)%MOD;
        cout<<(a%MOD*x%MOD)%MOD<<endl;
    }    
    return 0;
}