同余方程 ax≡1(mod b) & POJ 1061 青蛙的约会

---

题目：求ax%b=c最小正整数x解，题目中的c=1。
先感谢两位大犇ngncmh和笑巧。

---

对于一般的问题，我们通常有两种做法：

1) Baby Step Giant Step（BSGS）

定义Block为一个适中的常数，假设我们知道了[0,Block−1]内i∗a%b的值，就可以直接把区间分成：[Block,2∗Block−1]，[2∗Block,3∗Block−1]……
然后我们就可以对题目进行推导(ﾉ*･ω･)ﾉ：

p=i∗M+k(0<=k<M)

c=a∗p%b=a∗(i∗M+k)%b=1

a∗k%b=1−a∗i∗M

接下来就开始枚举i，如果出现符合题意的值，就可以得到最小整数解x了。而上述求得的[0,Block−1]内i∗a%b的值，可以用map映射。

#include<iostream>#include<algorithm>#include<cstring>#include<map>using namespace std;const long long M=100000;const long long P=2e9;map<long long,int>Map;int main(){    long long a,b,c;    cin>>a>>b;    Map.clear();    for(long long i=0;i<M;i++){        c=a*i%b;        if(c==1){            cout<<i<<endl;            return 0;        }        if(Map.find(c)==Map.end())Map[c]=i;    }/*    p=i*M+k(0<=k<M),    c=a*p%b=a*(i*M+k)%b=1;    a*k%b=1-a*i*M;*/    for(long long i=M;i<=P;i+=M){        c=(1LL-a*i)%b;        if(c<0)c+=b;        if(Map.find(c)!=Map.end()){            cout<<Map[c]+i<<endl;            return 0;        }    }    return 0;}

2)扩展欧几里得算法（extended gcd）

首先变形原式：

ax≡d%b−>ax+by=d(x,y∈Z)−>adx+bdy=1

即adx≡1(%bd)。
如果ad%bd≠1，说明不可能会有1的解（例如ad，bd均含有因子2，则取模后的结果只可能是2的倍数），所以ad%bd=1，由上述推论也可知gcd(ad,bd)∣1（整除）。
接着利用欧几里得算法可以求出一组特解x,y，所有的解属于y={x+bdk∣k∈Z}，根据求出的特解可以得到最小正整数解。

#include<iostream>#include<algorithm>using namespace std;typedef long long LL;LL a,b,t,x,y;//ax+by=tvoid gcd(LL a,LL b,LL &t,LL &x,LL &y){    if(!b){t=a;x=1;y=0;}//边界：gcd(a,0)=1*a+0*0=a    else{gcd(b,a%b,t,y,x);y-=x*(a/b);}}int main(){    cin>>a>>b;    gcd(a,b,t,x,y);    if(t!=1)cout<<"No answer"<<endl;    else cout<<(x%b+b)%b<<endl;    return 0;}

---

同余方程有三条定理：
定理一：若gcd(a,b)=d，则必能找到整数集中的k和l，使d=a∗x+b∗y；
定理二：若gcd(a,b)=1，则方程ax≡c(%b)在[0,b−1]上有唯一解；
定理三：若gcd(a,b)=d，则方程ax≡c(%b)在[0,bd−1]上有唯一解。

---

同余方程运用： 青蛙的约会

按照上述方法推出：(n−m)t%L=(x−y)，如果(x−y)%gcd(n−m,L)=0，方程有解。

其实听起来还是挺难懂。对于欧几里得函数的理解建议看这篇同类型的文章：NOI Openjudge 4975 两只鼹鼠。（事实证明理解清楚需要哪些部分后再写一遍真TM水）

#include<cstdlib>#include<iostream>#include<algorithm>using namespace std;typedef long long LL;LL gcd(LL a,LL b,LL &x,LL &y){    if(!b){x=1;y=0;return a;}//gcd(a,0)=1*a+0*0=a    else{        LL ans=gcd(b,a%b,x,y);        LL t=x;x=y;y=t-(a/b)*y;        return ans;    }}//扩展欧几里得算法 int main(){    LL x,y,m,n,L,x0,y0;    cin>>x>>y>>m>>n>>L;    //(x+tm)-(y+tn)=kL->(n-m)*t+L*k=(x-y)=>(n-m)t%L=(x-y)    //if (x-y)%gcd(n-m,L)==0,方程有解    int ans=gcd((n-m),L,x0,y0);//不能写成gcd((n-m)/(x-y),L/(x-y),x0,y0)    if((x-y)%ans!=0)cout<<"Impossible"<<endl;    else{//      cout<<x0<<' '<<y0<<endl;        x0=x0*(x-y)/ans;        LL r=L/ans;        cout<<(x0%r+r)%r<<endl;    }    return 0;}// 16/10/31 更新一版写法#include <cstdio>template <class temp>inline temp _abs(temp a){    return a<0?-a:a;}template <class temp>inline void swap(temp &a,temp &b){    temp t=a;a=b;b=t;}int exgcd(int a,int b,long long &x,long long &y){    if(!b){x=1,y=0;return a;}    int gcd=exgcd(b,a%b,y,x);    y-=(a/b)*x;    return gcd;}int query(int a,int b,int c){    long long x,y;    int gcd=exgcd(a,b,x,y);    if(c%gcd)return -1;    x*=c/gcd;    b=_abs(b/gcd);    return (int)(x%b+b)%b;}int main(){    int a,b,n,m,L;    scanf("%d %d %d %d %d",&a,&b,&n,&m,&L);    if(n<m)swap(n,m),swap(a,b);    int x=query(n-m,L,b-a);    if(!~x)puts("Impossible");    else printf("%d\n",x);}

---

总结：同余方程的重要部分是公式变形、欧几里得函数使用两处，解决的问题应能得到ax+by=c的变形。