扩展gcd
普通的gcd(a,b)算法就是求a和b的最大公约数。
int gcd(int a,int b) {return b?gcd(b,a%b):0;}
什么是扩展gcd,这是一种可以求一种特别的方程,已知整数a,b,c,求一组整数(x,y),满足 a x + b y = c ax+by=c ax+by=c
那么……我们先解决这个问题, g = g c d ( a , b ) , a > b g=gcd(a,b),a>b g=gcd(a,b),a>b,求一组整数解(x,y),满足 a x + b y = g ax+by=g ax+by=g。
首先,gcd(a,0)=0,此时必有x=1,y一般取0。
然后利用递归,假设我们已经求出 b ∗ x + ( a m o d b ) ∗ y = g c d ( a , b ) b*x+(a\ mod\ b)*y=gcd(a,b) b∗x+(a mod b)∗y=gcd(a,b)的解为(x0,y0),现在推a,b的答案。
首先,有 a ∗ x + b ∗ y = b ∗ x 0 + ( a m o d b ) ∗ y 0 a*x+b*y=b*x0+(a\ mod\ b)*y0 a∗x+b∗y=b∗x0+(a mod b)∗y0
又因为有 a m o d b = a − b ∗ ( a / b ) a\ mod\ b=a-b*(a/b) a mod b=a−b∗(a/b)(/代表整除,下同)
所以有$ax+by=bx0+ay0-b*(a/b)*y0 $
然后,就可以,移项 a ∗ ( x − y 0 ) = b ∗ ( x 0 − y − y 0 ∗ ( a / b ) ) a*(x-y0)=b*(x0-y-y0*(a/b)) a∗(x−y0)=b∗(x0−y−y0∗(a/b))
然后解决这个方程的最快方式,就是让左边=右边=0,这样下来,又有
x = y 0 x=y0 x=y0
y = x 0 − y 0 ∗ ( a / b ) y=x0-y0*(a/b) y=x0−y0∗(a/b)
这样就好了
注意求出一组(x0,y0)后,就可以直接求出任意的其它解,具体就这样(应该知道k为整数吧):
x = x 0 − k ∗ b / g c d ( a , b ) x=x0-k*b/gcd(a,b) x=x0−k∗b/gcd(a,b)
y = y 0 + k ∗ a / g c d ( a , b ) y=y0+k*a/gcd(a,b) y=y0+k∗a/gcd(a,b)
那么对于刚开始说的 a x + b y = c ax+by=c ax+by=c,首先判断是否满足 g c d ( a , b ) ∣ c gcd(a,b)|c gcd(a,b)∣c,如果能,答案乘上个c/gcd(a,b)就行了。
代码如下:
void exgcd(LL a,LL b,LL &d,LL &x,LL &y)
{
if (!b) {d=a; x=1; y=0;}
else {exgcd(b,a%b,d,y,x); y-=x*(a/b);}
}
(PS:据说这种求法可以保证|x|+|y|最小,不过我也不知道,各位神犇如果知道求告知。)
模线性方程
给出整数a,b,p求最小正整数x,满足 a x ≡ b ( m o d p ) ax\equiv b(mod\ p) ax≡b(mod p)
来,转化一下,设存在整数k满足 a ∗ x − b = k ∗ p a*x-b=k*p a∗x−b=k∗p
那么必有 a ∗ x − p ∗ k = b a*x-p*k=b a∗x−p∗k=b
这就又可以用扩展gcd来求了,对应的, a = a , b = p , c = b a=a,b=p,c=b a=a,b=p,c=b,则最后求出的x就是一个可行解
但是大多数时候要求会求出最小正整数解。
设d=gcd(a,b),两个相邻可行解的间距为L,则设当前已求出一个可行解x0,则
a ( x 0 + L ) ≡ b ( m o d p ) a(x0+L)\equiv b(mod\ p) a(x0+L)≡b(mod p)
a x ≡ b ( m o d p ) ax\equiv b(mod\ p) ax≡b(mod p)
两式相减,得
a ∗ L ≡ 0 ( m o d p ) a*L\equiv 0(mod\ p) a∗L≡0(mod p) 即p|a*L
则L的最小值就是p/d。
那么x的最小正整数也有公式了。
x = ( x 0 m o d L + L ) m o d L x=(x0\ mod\ L+L)\ mod\ L x=(x0 mod L+L) mod L
说明,当b=1时,x就是a在模p情况下的逆元,可写成 x − 1 x^{-1} x−1,不过也有O(n)递推求出1到n的,比这里的O(n*logn)好。在此不表。
更新于2017年8月23日19:26:25
加例题一道 POJ 2115
求最小正整数x,满足
(
a
+
c
∗
x
)
m
o
d
2
k
=
b
(a+c*x)\ mod\ 2^k=b
(a+c∗x) mod 2k=b
令
a
=
c
;
b
=
b
−
a
;
c
=
2
k
a=c;\ b=b-a;\ c=2^k
a=c; b=b−a; c=2k,就是经典的模线性方程
a
x
≡
b
(
m
o
d
c
)
ax\equiv b(mod\ c)
ax≡b(mod c)
注意开long long
#include<cstdio>
#define LL long long
using namespace std;
LL x,L,a,b,c,p,k;
inline void readL(LL &x){
x=0; char ch=getchar();
while ('0'>ch||ch>'9') ch=getchar();
while ('0'<=ch&&ch<='9') {x=x*10+ch-'0'; ch=getchar();}
}
void exgcd(LL a,LL b,LL &d,LL &x,LL &y){
if (b) {exgcd(b,a%b,d,y,x); y-=x*(a/b);}
else {x=1; y=0; d=a;}
}
int main()
{
freopen("looop.in","r",stdin);
freopen("looop0.out","w",stdout);
readL(a); readL(b); readL(c); readL(k);
while (a||b||c||k){
b-=a; a=c; c=(LL)1<<k;
exgcd(a,c,p,x,k);
if (b%p) printf("FOREVER\n");
else{L=c/p; x=(x*(b/p))%c; x=(x%L+L)%L; printf("%lld\n",x);
}
readL(a); readL(b); readL(c); readL(k);
}
return 0;
}