欧几里得算法与逆元

本文转载自http://blog.csdn.net/FSAHFGSADHSAKNDAS/article/details/51027338

欧几里得算法与逆元

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

（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 long
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;
}

求解方程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 long
using 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)%p

这样就可以在计算过程中边取模边做除法了。

求法

[b]（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 1000000
using 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
*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 1000000
using 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
*ans[M]%R);
}
return 0;
}