欧几里得算法与逆元

（一）欧几里得与扩展欧几里得

（I）欧几里得辗转相除求gcd

描述

欧几里得辗转相除求最大公因数的表达式如下

gcd(a,b)={agcd(b,amodb)b=0b≠1

代码

int gcd(int a,int b){    if(b==0)return a;    else return gcd(b, a%b);}

（II）扩展欧几里得

描述

定理：对于任意整数a，b，都存在一组整数x、y使得ax+by=gcd(a,b)成立
于是设 bx’+(a%b)y’=gcd(b,a%b)
则ax+by=bx’+(a%b)y’
可以推导出x=y’,y=x’-[a/b]*y’
当b=0，a≠0时，gcd(a,b)=a，则x=1,y=0，
由此便可以向上推出ax+by=gcd(a,b)中x,y的值

代码

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

求解方程ax+by=c

定理：当且仅当当gcd(a,b)|c时，方程有整数解，设d=gcd(a,b)，令a’=a/d，b’=b/dc’=c/d，方程两边同时除以d就得到a’x+b’y=c’，我们可以先用扩展欧几里得求出a’x+b’y=1的一组解x0和y0，再把x0和y0同时乘c’就得到了a’x+b’y=c’的一组解，记为xx和yy，即原方程的一组解。方程的所有解可以表示为x=xx+kb’，y=yy-ka’，其中k∈Z
模板题 poj1061 青蛙的约会
http://poj.org/problem?id=1061
设青蛙的初始坐标分别为A、B，其他的字母与题中一样，设相遇时间为t
则有

A+mt=B+nt(modL)

整理得

(m−n)t+kL=B−A

令a=m-n,b=k,c=B-A，求出ax+by=c的最小正整数解x就行了
注意一个问题，此处求出xx后，如果用while(xx+kb’<0)k++;的话会超时，由于只需要最小正整数解，所以x=(xx%|b’|+|b’|)%|b’|
代码：

//poj1061 青蛙的约会 扩展欧几里得#include <iostream>#define abs(x) (x>0?x:-x)#define ll long longusing namespace std;ll m, n, A, B, L;ll exgcd(ll a, ll b, ll &x, ll &y){    if(!b){x=1;y=0;return a;}    ll t, d;    d=exgcd(b,a%b,x,y);    t=x;x=y;y=t-a/b*y;    return d;}int main(){    ll a, b, c, d, x, y;    cin >> A >> B >> m >> n >> L;    a=m-n;b=L;c=B-A;    d=exgcd(a,b,x,y);    if(c%d!=0){cout<<"Impossible"<<endl;return 0;}    exgcd(a/d,b/d,x,y);    x*=c/d;    x=(x%abs(b/d)+abs(b/d))%abs(b/d);    cout << x << endl;    return 0;}

（二）逆元

用处

为啥要用逆元？
在各种题目中，又是让你求答案对xxxxx取模的值，那么这肯定是因为答案太大，无法存储。在NOIP即以下的比赛中，凡是让你取模的题目，都不会涉及除法。
我们知道(a+b)%p=(a%p+b%p)%p，ab%p=(a%p*b%p)%p，但是绝对不能说(a/b)%p=[(a%p)/(b%p)]%p
随便试几组就会发现这样是错误的
这时就要用到逆元

ax≡1(modp)

中x的最小正整数解就是a关于p的逆元，或者说a在模p意义下的数论倒数，记为inv[a]。即是说，(a/b)%p=(a%p*inv[b])%p
这样就可以在计算过程中边取模边做除法了。

求法

（i）扩展欧几里得

我上面介绍扩展欧几里得就是为了这里用的
令ax+kp=1
明显当gcd(a,k)≠1时，方程无解，逆元不存在（这里我有点不懂）
否则就用扩展欧几里得求最小正整数x
代码在后面会有

（ii）费马小定理

费马小定理：当p为素数且a≠p时，有ap−1modp=1，继续推导

ap−2×a=1(modp)

ap−2=1a(modp)

于是inv[a]=ap−2modp，
注意：要求p必须是质数且a不等于p
ps：实践证明，无论何种求法，当把p本身的逆元设为0时，答案总是正确的
例题：hdu1576 A/B
http://acm.hdu.edu.cn/showproblem.php?pid=1576
就是给你A和B，其中A是某个数对9973取模的值，让你求A/B，A和B互质，求出B的逆元，然后输出答案就行了
代码（扩展欧几里得和费马小定理）：

//hdu1576 A/B 求逆元#include <cstdio>using namespace std;int A, B, p=9973,T;//法一：扩展欧几里得/*void exgcd(int a, int b, int &x, int &y){    int t;    if(!b){x=1;y=0;return;}    exgcd(b,a%b,x,y);    t=x;x=y;y=t-a/b*x;}int inv(int a){    int x, y;    exgcd(a,p,x,y);    return (x+p)%p;}*///法二：费玛小定理int pow(int a, int b){    int t=a, ans;    for(t=a,ans=1;b;b>>=1,t=t*t%p)if(b&1)ans=ans*t%p;    return ans;}int inv(int a){    return pow(a%p,p-2);}int main(){    while(scanf("%d",&T)!=EOF)    {        while(T--)            scanf("%d%d",&A,&B),            printf("%d\n",(A*inv(B)%p));    }    return 0;}

（iii）线性求逆元

线性求逆元的递推式如下：

inv[i]=inv[Mmodi]×(M−M/i)modM

其中i<M
证明
设

M=ki+b

那么

k=M/i,b=Mmodi

ki+b=0(modM)

ki=−b(modM)

1ki=−1b(modM)

1i=−kb(modM)

即

inv[i]=inv[Mmodi]∗(M−M/i)modM

得证
这种做法也有要求，M必须是质数
例题 bzoj2186 [Sdoi2008]沙拉公主的困惑
http://www.lydsy.com/JudgeOnline/problem.php?id=2186
题目详解（我写的）：
http://blog.csdn.net/fsahfgsadhsakndas/article/details/50995536
代码再粘一遍：

//bzoj2186 sdoi2008沙拉公主的困惑#include <cstdio>#define ll long long#define maxn 10000000#define maxk 1000000using namespace std;ll p[maxk], R, inv[maxn+100], N, M, tot, fac[maxn+100],    ans[maxn+100], T;bool no[maxn+100];void getinv(){    int i, m=(R<maxn?R:maxn); inv[1]=1;    for(i=2;i<m;i++)inv[i]=(ll)(R-R/i)*inv[R%i]%R;}void shai(){    int i, j;    for(i=2;i<=maxn+1;i++)    {        if(!no[i])p[++tot]=i;        for(j=1;i*p[j]<=maxn;j++)        {            no[i*p[j]]=true;            if(i%p[j]==0)break;        }    }}void getfac(){    int i;    fac[1]=1;    for(i=2;i<=maxn;i++)fac[i]=fac[i-1]*i%R;}void getans(){    int i, j;    ll t;    ans[1]=1;    for(i=1;i<=tot;i++)    {        t=(ll)ans[p[i]-1]*(p[i]-1)%R*inv[p[i]%R]%R;        for(j=p[i];j<p[i+1];j++)ans[j]=t;    }}int main(){    int i;    scanf("%d%d",&T,&R);    shai();    getinv();    getfac();    getans();    while(T--)    {        scanf("%d%d",&N,&M);        printf("%lld\n",fac[N]*ans[M]%R);    }    return 0;}

（iv）总结

扩展欧几里得、费马小定理+快速幂、线性求逆元是逆元的三种求法
扩展欧几里得：要求p和a必须互质，时间复杂度O(logpa)
费马小定理+快速幂：要求p必须是质数，时间复杂度O(logp2)
线性求逆元：要求p必须是质数，时间复杂度O(P)
以上可以看出，如果求一个数的逆元，扩展欧几里得是最佳的选择，如果需要N个数的逆元，线性递推式最好的方式，费马小定理+快速幂似乎没什么用武之地
还可以看出，后两种都是要求p是质数的，当p不是质数，但p和a互质时，可以选择扩展欧几里得
例题 CodeVS2301 沙拉公主的困惑
http://codevs.cn/problem/2301/
　　这道题。。。。乍看好眼熟，你发现最后一句话：“R是一个合数”
　　这下就不能用线性递推求逆元了，只能把所有素数的逆元用欧几里得求出。

//CodeVS2301 沙拉公主的困惑 欧几里得求逆元 #include <cstdio>#define min(a,b) (a<b?a:b)#define ll long long#define maxn 10000000#define maxk 1000000using namespace std;ll p[maxk], R, inv[maxn+100], N, M, tot, fac[maxn+100],    ans[maxn+100], T;int c;bool no[maxn+100];void exgcd(ll a, ll b, ll &x, ll &y){    int t;    if(!b){x=1;y=0;return;}    exgcd(b,a%b,x,y);    t=x;x=y;y=t-a/b*y;}void getinv(){    ll i, x, y, d;    for(i=1;i<tot;i++)        exgcd(p[i],R,x,y),inv[i]=(x+R)%R;}void shai(){    int i, j;    for(i=2;i<=maxn+1;i++)    {        if(!no[i])p[++tot]=i;        for(j=1;i*p[j]<=maxn&&j<=tot;j++)        {            no[i*p[j]]=true;            if(i%p[j]==0)break;        }    }}void getfac(){    int i;    fac[1]=1;    for(i=2;i<=maxn;i++)fac[i]=fac[i-1]*i%R;}void getans(){    int i, j;    ll t;    ans[1]=1;    for(i=1;i<=tot;i++)    {        t=(ll)ans[p[i]-1]*(p[i]-1)%R*inv[i]%R;        for(j=p[i];j<p[i+1];j++)ans[j]=t;    }}int main(){    int i;    scanf("%d%d",&T,&R);    shai();    getinv();    getfac();    getans();    while(T--)    {        scanf("%d%d",&N,&M);        printf("%lld\n",fac[N]*ans[M]%R);    }    return 0;}