从欧几里得讲起……
欧几里得的辗转相除法计算的是两个自然数a和b的最大公约数g,意思是能够同时整除a和b的自然数中最大的一个。
两个数的最大公约数通常写成GCD(a, b),或者简写成(a, b)。
(Greatest Common Divisor)
计算方法……辗转相除!
作为数论中一个很基础的内容,它的形式很漂亮,代码也很干净。
gcd(a,b)=gcd(b,a
证明过程
令r=a则有a=k×b+r故可知r=a−k×b
现假设d是a,b的一个公约数
则有d∣a , d∣b
又由上式可得d∣(a−k×b)即d∣r⟺d∣(a
由此证明d也是b和(a的公约数
即gcd(a,b)=gcd(b,a%b)
注:d∣a表示d整除a
代码实现
long long gcd(long long a,long long b)
{ return b == 0 ? a : gcd(b,a%b); }
贝祖等式……问题的开端!
贝祖等式(又称裴蜀等式)是一个形如ax+by=m的方程。
而有整数解(x,y)⟺m=k×gcd(a,b),k∈Z。
裴蜀等式有解时必然有无穷多个解。
补充:特别的,ax+by=1⟺gcd(a,b)=1
而利用辗转相除法可以得到这一方程的整数特解。
而且不难发现,大多数情况下x和y总是一正一负。
一段巧妙的证明
从ax+by=gcd(a,b)开始……
由于长得漂亮对称的结构,我们不妨设a>b
- 当b=0时,gcd(a,b)=a这个时候我们的等式就变成了ax=a的形态不难看出这时候的一组解(x,y)=(1,0)
- 当b=0时候,我们可以得到
ax1+by1=gcd(a,b)
bx2+(a%b)y2=gcd(b,a%b)
又由于gcd(a,b)=gcd(b,a%b)
所以有ax1+by1=bx2+(a%b)y2
定义可得a%b=a−⌊ba⌋×b将其带入上式ax1+by1=bx2+(a−⌊ba⌋×b)y2
ax1+by1=bx2+ay2−(⌊ba⌋×b)y2
由a和b的对应系数相等,所以可推出
x1=y2
y1=x2−⌊ba⌋y2
综上,我们推导出了求该等式特解的边界条件和递推过程。
代码实现
注意 gcd , x , y 变量均使用引用传参,这是为了每一步都能改变原始数据。
void exgcd(long long a, long long b, long long &gcd, long long &x, long long &y)
{
if (b == 0)
{
gcd = a;
x = 1;
y = 0;
return;
}
exgcd(b, a % b, gcd, x, y);
long long xt = x, yt = y;
x = yt;
y = xt - a / b * yt;
return;
}
应用举例
求解同余方程的最小正整数解ax≡1(modb),给出a和b。
分析
ax≡1(modb)⟺ax+by=1
其中y是我们引入的一个整数而且十有八九应该是负数。
但是有一个问题,ax+by=1又不是ax+by=gcd(a,b)
这是怎么牵扯上关系的呢?还不是你乱点鸳鸯谱
由最大公约数的定义,总有a∣gcd(a,b)而且b∣gcd(a,b)
如果有x,y∈Z则必定(ax+by)∣gcd(a,b)
即若有ax+by=m则m∣gcd(a,b)
当m=1时,就有1∣gcd(a,b)
故而gcd(a,b)就只能等于1了
但是我们求出的x十有八九是负数也有可能太大,怎么得到最小正整数解呢?
ax+by=1
ax+by+k×ab−k×ab=1
a(x+kb)+b(y−ka)=1
故而我们使用(x′%b+b)%b的方式计算最终结果
这里一开始我不太明白
后来想想,我们要求的是x′在加上k个b之后最小的整数结果
前提:x′<0,b>0
由x′%b可以将其结果限制在[−b,0]
当−b<x′%b<0时候加b便是最终结果了
当x′时候加b再对b取模便是0不会出现错误
这里涉及到负数的取模运算,但这个计算的结果是百家争鸣的,我们先以C++中的结果为准
在C++中,有−7%3=−1
以a=3,b=10为例解得x=−3,y=1
计算得到(−3%10+10)%10=7
代码实现
void exgcd(long long a, long long b, long long &gcd, long long &x, long long &y)
{
if (b == 0)
{
gcd = a;
x = 1;
y = 0;
return;
}
exgcd(b, a % b, gcd, x, y);
long long xt = x, yt = y;
x = yt;
y = xt - a / b * yt;
return;
}
int main()
{
long long a,b;
long long g,x,y;
cin >> a >> b;
exgcd(a,b,g,x,y);
x = (x % b + b)%b;
cout << x << endl;
return 0;
}