欧几里得算法与逆元

本文转载自 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[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 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[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 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[N]*ans[M]%R);
    }
    return 0;
}