数论 —— 逆元（费马小定理/扩展欧几里得）

拓展欧几里得：

(证明过程参考自百度百科)

原式: ax+by=gcd(a,b)(假设a≥b)

当b=0时有gcd(a,b)==a,此时x=1,y=0

当b不为0时,根据欧几里得定理gcd(a,b)=gcd(b,amodb)可得ax+by=gcd(a,b)=gcd(b,amodb)=bx′+(amodb)y′,即 

ax+by=bx′+(amodb)y′=bx′+(a−b∗⌊a/b⌋)y′

移项得 
ax+by=bx′+(amodb)y′=ay′+b(x′−⌊a/b⌋y′)

根据恒等定理,有 
{x=y′y=x′−⌊a/b⌋y′

这有什么用呢?x′和y′还是不知道呀. 
重新来看看我们得到的两个等式.x和y是gcd(a,b)=ax+by的解,而x’和y’是在对gcd(a,b)按欧几里德算法进行一步后的结果对应的贝祖等式gcd(b,amodb)=bx′+(amodb)y′的解.也就是说,gcd(a,b)对应的贝祖等式的解x,y可以由gcd(b,amodb)对应等式的解x’,y’计算得出 
由于欧几里德算法最后一步为gcd(d,0)=d,此时对应的等式的解为x=1,y=0,因此只要如上述代码,从gcd(d,0)往前处理,在进行欧几里德算法的递归的时候根据相邻两次调用间x,y和x’,y’的关系计算即可求出ax+by=gcd(a,b)的解.
更进一步,对于任意不定式ax′+by′=c,只需要在等式ax+by=gcd(a,b)=d两边乘上c/d即可得到解为x′=x∗c/d,y′=y∗c/d

如何得到所有解? 
实际上在之前的计算和证明中我们得到的只是不定方程的一组解,那么怎样得到所有解呢?对于一般形式ax+by=c有通解x=p+kb,y=q−ka(k为任意整数).(证明略,只要代入一下就知道为什么通解是这个了)

int extend_GCD(int a,int b,int &x,int &y){
int tmp,d;
if(b==0){
x=1;
y=0;
return a;
}
int res=extend_GCD(b,a%b,x,y);
tmp = x;
x = y;
y = tmp - (a/b)*y;
return res;
}

void Extend_Gcd(int a,int b,int &n,int &x,int &y){
    if(!b){
        n=a,x=1,y=0;
        printf("1:%d %d\n",x,y);
    }else{
        Extend_Gcd(b,a%b,n,y,x);
        y-=x*(a/b);
        printf("---%d %d\n",x,y);
    }
}

费马小定理：

x与p互质:

若gcd(x, p) == 1

x^(p-1) == 1(mod p)

快速幂求解即可

逆元：
a*x  == 1（mod m）
x是相对于m的a的逆元
a是相对于m的x的逆元

1.费马小定理：若gcd(x,p) == 1
x^(p-1) == 1 (mod p)
x*x^(p-2) == 1(mod p)
所以x是x^(p-2)关于p的逆元LL quick(LL x,LL n){
LL res=1;
while(n){
if(n%2)res=x*x%Mod;
x*=x%Mod;
n>>=1;
}
return res;
}
LL inv(LL x){
return quick(x,Mod-2);
}
2.扩展欧几里得求逆元：
    a*x == 1(mod p)

    a*x + kp == 1
扩展欧几里得解得extend_gcd(a,p)
求解得 a*x0+p*y0 = 1
x0是a的一个特解 
通解 x == x0+pk 

3.
当p是个质数的时候有

inv(a) = (p - p / a) * inv(p % a) % p
证明：
设x = p % a,y = p / a
于是有 x + y * a = p
(x + y * a) % p = 0
移项得 x % p = (-y) * a % p
x * inv(a) % p = (-y) % p
inv(a) = (p - y) * inv(x) % p
于是 inv(a) = (p - p / a) * inv(p % a) % p

然后一直递归到1为止，因为1的逆元就是1 

LL inv(LL t, LL p) {//求t关于p的逆元，注意:t要小于p，最好传参前先把t%p一下
return t == 1 ? 1 : (p - p / t) * inv(p % t, p) % p;
}

可以预处理出前N项逆元：

#include<cstdio>
const int N = 200000 + 5;
const int MOD = (int)1e9 + 7;
int inv
;
int init(){
inv[1] = 1;
for(int i = 2; i < N; i ++){
inv[i] = (MOD - MOD / i) * 1ll * inv[MOD % i] % MOD;
}
}
int main(){
init();
}