来学数学!
拓展欧几里得算法 (EXGCD)
求 ax+by=gcd(a,b) 的整数解。
例:求 6x+15y=gcd(6,15)=3 的整数解。
有一组整数解 x=3,y=−1。(解不唯一)
求一组可行解
设存在 x1,y1,x2,y2,满足
{ax1+by1=gcd(a,b)bx2+(amodb)y2=gcd(b,amodb)
因为 gcd(a,b)=gcd(b,amodb),所以 ax1+by1=bx2+(amodb)y2。
又因为 amodb=a−b⌊ba⌋,所以 ax1+by1=bx2+(a−b⌊ba⌋)y2。
整理得 ax1+by1=ay2+b(x2−⌊ba⌋y2)。
即 x1=y2,y1=x2−⌊ba⌋y2。
则我们可以用 a′=b,b′=amodb 的方程的解来反向推出系数为 a,b 方程的解。
边界条件: b=0 时, 解为 x=1,y=0。
由此,我们得到重要结论 1:
重要结论 1
若 bx+(amodb)y=gcd(b,amodb) 的一组整数解为 x=x0,y=y0,则 ax+by=gcd(a,b) 的一组整数解为 x=y0,y=x0−y0⌊ba⌋。
模板代码
由此,可以写出代码。
1 2 3 4 5 6
| void exgcd(int a,int b,int& x,int& y){ if(!b){x=1,y=0;return; } exgcd(b,a%b,x,y); swap(x,y); y -= (a/b)*x; }
|
使用引用来递归计算。x
和 y
,在回溯时一层一层更新。注意时间复杂度和欧几里得算法相同,为 O(logn)。
对其稍加改造还可以一并求出 gcd(a,b),并且可以减少掉这个 swap 操作。
1 2 3 4 5 6
| int exgcd(int a,int b,int& x,int& y){ if(!b){x=1,y=0;return a;} int g = exgcd(b,a%b,y,x); y -= (a/b)*x; return g; }
|
返回值即为 gcd(a,b)。
求出所有解
我们已经求出一组解,接下来尝试求出所有整数解。
设两组整数解为 x1,y1 和 x2,y2,则 ax1+by1=ax2+by2。
整理,a(x1−x2)=b(y2−y1)。
设 gcd(a,b)=g,左右两边同时除 g,a′(x1−x2)=b′(y2−y1)。
注意到此时 a′,b′ 互素,则 x1−x2=kb′(k∈Z),带回得 y1−y2=−ka′。
由此,我们得到重要结论 2:
重要结论 2
若 ax+by=gcd(a,b) 的一组整数解为 x=x0,y=y0,则其所有整数解为 x=x0+kb′,y=y0−ka′,式中 a′=a/gcd(a,b),b′=b/gcd(a,b),k∈Z 。
注意到在证明过程中我们没有用到 ax+by 等于什么,因此得到重要结论 2 一般形式:
重要结论 2 一般形式
若 ax+by=c (c 为任意整数) 的一组整数解为 x=x0,y=y0,则其所有整数解为 x=x0+kb′,y=y0−ka′,式中 a′=a/gcd(a,b),b′=b/gcd(a,b),k∈Z 。
求解线性丢番图方程
求方程 ax+by=c 的所有整数解。
事实上 c 不是 gcd(a,b) 的倍数时,此方程没有整数解。
下证:设 gcd(a,b)=g,方程两边同除 g,得 a′x+b′y=c/g。此时等号左边时整数,故当且仅当 c 为 g 倍数时有解。(若 ab=0,可以特殊判断)
在其他情况可以用上面的进行推导,得到重要结论 3:
重要结论 3
对于方程 ax+by=c,当且仅当 c 为 gcd(a,b) 倍数时有无数组整数解。
且若 ax+by=gcd(a,b) 的一组解为 x=x0,y=y0,则此方程的一组解为 x=x0gcd(a,b)c,y=y0gcd(a,b)c。
要求出这个方程的所有解用重要结论 2 一般形式即可。
例题
直线上的点:对于直线 l:ax+by+c=0,求有多少个整点 (x,y) 满足 x∈[x1,x2],y∈[y1,y2]?保证给定直线不与坐标轴垂直。
例:l:2x+3y+5=0,x1=y1=−5,x2=y2=5,共有 4 个满足要求的点,分别为 (2,−3),(5,−5),(−1,−1),(−4,1)。
首先当 c 不为 gcd(a,b) 倍数时,直线不经过任何整点,直接输出 0
。
否则就可以求出 ax+by=gcd(a,b) 的一组整数解 (x,y)。则一个整点的坐标可以表示为 (x+gcd(a,b)−c,ygcd(a,b)−c) 记作 (x0,y0)。
则有
x1≤x0+kb′≤x2
y1≤y0−ka′≤y2
求解不等式即可。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
| #include<bits/stdc++.h> using namespace std;
int exgcd(int a,int b,int& x,int& y){ if(!b){x=1,y=0;return a;} int g=exgcd(b,a%b,y,x); y-=x*(a/b); return g; }
int main(){ int a,b,c; int x1,x2,y1,y2; scanf("%d%d%d",&a,&b,&c); scanf("%d%d%d%d",&x1,&x2,&y1,&y2);
int x,y; int g = exgcd(a,b,x,y);
if(c%g !=0){ puts("0");return 0; } int x0 = x*-c/g, y0=y*-c/g; int ap = a/g,bp=b/g;
double l1 = (x1-x0)*1.0/bp; double r1 = (x2-x0)*1.0/bp; if(l1>r1)swap(l1,r1);
double l2 = -(y1-y0)*1.0/ap; double r2 = -(y2-y0)*1.0/ap; if(l2>r2)swap(l2,r2);
int l = max(ceil(l1),ceil(l2)); int r = min(floor(r1),floor(r2));
if(l>r)puts("0"); else printf("%d\n",r-l+1);
return 0; }
|
线性同余方程
求解方程 ax≡b(modn)。 a,b,n≤109。
注意到方程等价于 ax+ny=b,即转化为一个不定方程。
方程有解的条件为 b 为 gcd(n,a) 倍数。
此时若一个解为 x0,由上文的重要结论 2 可知,方程的所有解为 x0+kgcd(a,n)b。