GCD
Greatest Common Divisor，常缩写为 gcd，即我们在数学中熟悉的最大公约数。所谓约数，即指一个数的因子，比如 6 能被 3 整除，所以 3 是 6 的约数，一个数的约数的个数是有限的。一组数的公约数，是指同时是这组数中每一个数的约数的数，而最大公约数，则是指所有公约数里面最大的一个。
一、欧几里得算法
求两个数的最大公约数有很多方法，如质因数分解法、短除法、辗转相除法、更相减损法等，本文介绍的就是我们在算法竞赛中最常用的辗转相除法，即欧几里得算法，它的代码很简单：
// 欧几里得(递归版)
inline int gcd(int a, int b)
{
return b == 0 ? a : gcd(b,a%b);
}
// 欧几里得(非递归版)
inline int gcd2(int a,int b)
{
while(1)
{
if(a<b) {int t=a;a=b;b=t;}
if(b==0) return a;
int x=a%b; a=b; b=x;
}
}算法递归至 b==0（即上一步的 a%b==0) 的情况再返回值即可，上述算法也被称作欧几里得算法（Euclidean algorithm），特殊的，如果 gcd(a,b) = 1 ，我们称 a 与 b 互质。欧几里得算法的整个流程图大致如下：
为什么辗转相除法是成立的？我们不妨设 a > b ，可以发现如果 b 是 a 的约数，那么 b 就是两者最大的公约数。而如果 b 并无法整除 a ，我们可以得到 a = b × k + r ， k 为商(即 a/b )， r 为余数(即 a%b )，可知 r < b 。通过证明可知 gcd(a,b) = gcd(b, a%b)，过程如下：
设 a = b × k + r ， 显然有 r = a % b ，设 d 为 a、b 的公约数，则有 r = a - bk ， $\frac{r}{d}$ = $\frac{a}{d}$ - $\frac{b}{d}$ k 。 由右边的式子可知 $\frac{r}{d}$ 为整数，即 d 也能整除 r，所以对于 a、b 的公约数，它也会是 a % b 的公约数。
反过来也需要证明：设 a = b × k + r ，显然有 r = a % b ， 设 d 是 b、a % b 的公约数，我们依然可以得到等式 $\frac{a}{d}$ = $\frac{b}{d}$ k + $\frac{r}{d}$ ，因为右边的式子显然为整数，所以 $\frac{a}{d}$ 也为整数，即 d 也能整除 a ，所以对于 b、a % b 的公约数，它也会是 a 的公约数。
既然两式的公约数是相同的，那么最大公约数也是相同的，所以得到式子 gcd(a,b) = gcd(b,r) ，即 gcd(a,b) = gcd(b,a%b) 。我们也就得到了代码里面的那个递归求法。
虽然我们在一开始假设了 a > b ，但是在代码里并不需要去对传入参数做判断，是因为如果 a < b ， a % b 此时就等于 a ，所以函数会返回 gcd(b,a) ，其实就完成了一次交换。
欧几里得算法的时间复杂度为 $O(log{n})$ ，证明如下，当我们求 gcd(a,b) 时，会遇到两种情况：
- a < b ，这时 gcd(a,b) = gcd(b,a)
- a >= b ，这时 gcd(a,b) = gcd(b,a%b) ，而对于 a 取模会让 a 至少折半，这意味着这一过程最多发生 $O(log{n})$ 次
第一种情况发生后一定会发生第二种情况，因此第一种情况发生的次数一定不多于第二种情况，从而我们最多递归 $O(log{n})$ 次就可以得出结果。事实上假设我们用欧几里得算法去求斐波那契数列相邻两项的最大公约数，会让该算法达到最坏复杂度。
二、多个数的gcd
那怎么求多个数的最大公约数呢？显然答案一定是每个数的约数，那么也一定是每相邻两个数的约数。我们采用归纳法，可以证明，每次取出两个数求出答案后再放回去，不会对所需要的答案造成影响。
LCM
Least Common Multiple，常缩写为 lcm，即我们在数学中熟悉的最小公倍数。一个整数能够被另一个整数整除，这个整数就是另一整数的倍数，比如 6 能被 3 整除，所以 6 是 3 的倍数，一个数的倍数是无限的。一组数的公倍数，是指同时是这组数中每一个数的倍数的数，而最小公倍数，则是指所有公倍数里面最小的一个。
一、两个数的lcm
在知道两个数的gcd得前提下，lcm是很容易求的，公式是 d = $\frac{a×b}{gcd(a,b)}$ ，代码如下：
// 欧几里得算法
inline int gcd(int a, int b)
{
return b == 0 ? a : gcd(b,a%b);
}
// 最小公倍数
inline int lcm(int a, int b)
{
// 公式是a*b/gcd(a,b),因为直接a*b容易爆范围,所以先算除法
return a / gcd(a,b) * b;
}接下来给出证明，首先我们求出 gcd(a,b) ，显然我们能知道 $\frac{a}{gcd(a,b)}$ 和 $\frac{b}{gcd(a,b)}$ 互质，此时 $\frac{a}{gcd(a,b)}$ 和 $\frac{b}{gcd(a,b)}$ 的最小公倍数是它们两的乘积，即 $\frac{a}{gcd(a,b)}$ × $\frac{b}{gcd(a,b)}$ 。现在如果要让它变成 a 和 b 的共同倍数，则最少乘上 gcd(a,b) 即可，所以 a 和 b 的最小公倍数为 d = $\frac{a×b}{gcd(a,b)}$ 。
二、多个数的lcm
可以发现，当我们求出两个数的 gcd 时，求最小公倍数是 O(1) 的复杂度，那么对于多个数，我们其实没有必要求一个共同的最大公约数再去处理。前面提到求多个数 gcd 时，当我们算出两个数的 gcd 后，我们会将它放回序列对后面的数继续求解，那么求最大公倍数则可以转换一下，我们每两个数求 lcm ，求完后把 lcm 放回序列序列即可。
HDU 1019 Least Common Multiple
题目描述
给出 n 个数字，求这 n 个数字的最小公倍数
思路
每两个数求一次 lcm ，求完后的结果放回序列中，直到序列只剩一个数
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
inline ll gcd(ll a,ll b)
{
return b==0?a:gcd(b,a%b);
}
inline ll lcm(ll a,ll b)
{
return a/gcd(a,b)*b;
}
ll a[100010];
int main()
{
ll n;
while(scanf("%lld",&n)!=EOF)
{
while(n--)
{
ll m;
scanf("%lld",&m);
for(ll i=1;i<=m;i++)
scanf("%lld",&a[i]);
for(ll i=1;i<=m-1;i++)
a[i+1]=lcm(a[i],a[i+1]);
printf("%lld\n",a[m]);
}
}
return 0;
}拓展GCD
一、拓展gcd原理
那么什么是拓展欧几里得呢？拓展gcd是一种算法，它可以在辗转相除求最大公约数的过程中求出不定方程 ax + by = c 的一组解，它的流程跟普通gcd很相似
可以注意到上图倒数第二行的 3 + 6 × 3 = 21 可以改写成 3 = 6 × (-3) + 21 ，也就是说 3 可以表示为 6 和 21 的线性组合。接着看倒数第三行的 6 + 21 × 1 = 27 ，说明 6 可以被表示为 21 和 27 的线性组合，即 6 = 21 × (-1) + 27 ，具体的我们可以得到：3 = 6 × (-3) + 21 = (21 × (-1) + 27) × (-3) + 21 = 27 × (-3) + 21 × 4 。这样一路推下去，我们最终可知 3 可以表示为75 和 48 的线性组合，那么方程 75x + 48y = 3 就能找到解了。
我们知道 3 = gcd(75,48) ，经过上面的流程，我们知道了求 ax + by = gcd(a,b) 的一种方法，那如果 c 是其他数，可否找到解呢？事实上 c 必须是 gcd(a,b) 的整数倍，我们把方程 ax + by = c 的两边除以 gcd(a,b) 可得 $\frac{a}{gcd(a,b)}$x + $\frac{b}{gcd(a,b)}$y = $\frac{c}{gcd(a,b)}$ ，方程左边显然是整数，那么方程右边也必须是整数，所以可知 c 必须是 gcd(a,b) 的整数倍
实际上，上述过程并不是很严谨的证明了一个数论定理：裴蜀定理(贝祖定理)
裴蜀定理：
设 a、b 为正整数，则对于任意的整数 x、y 而言，ax + by 一定是 gcd(a,b) 的倍数，
特别的，当且仅当 c 是 gcd(a,b) 的倍数时，方程 ax + by = c 一定有整数解重要推论：
a、b 互质的充分必要条件是存在整数 x、y 使 ax + by = 1
把刚刚提到过的流程画下来是这样的（右下角按规律添了两行）：
可以看到，通过求 $bx_0 + (a%b) y_0 = c$ 的解，可以得出 ax + by = c 的解。实际上，前者等价于 $bx_0 + (a - b \lfloor a/b \rfloor) y_0 = c$ ，也就是 $ay_0 + b(x_0 - \lfloor a/b \rfloor y_0) = c$ ，根据待定系数法可得到下列通项公式：
和普通辗转相除法相同，当 b = 0 时递归结束，由上图右侧底部的式子可知，此时应该令 x = 1 ， y = 0 。拓展gcd的代码如下：
// 拓展gcd
int exgcd(int a, int b, int &x, int &y)
{
if (b == 0)
{
x = 1;
y = 0;
return a;
}
int d = exgcd(b, a % b, x, y), x0 = x, y0 = y;
x = y0;
y = x0 - (a / b) * y0;
return d;
}可以看到，确实就是在欧几里得算法的基础上进行了一些拓展。利用C++的引用特性，我们得以简洁地实现这一算法。其实也可以有简单一点点的写法
int exgcd(int a, int b, int &x, int &y)
{
if (b == 0)
{
x = 1;
y = 0;
return a;
}
int d = exgcd(b, a % b, y, x); //这里交换了x和y
y -= (a / b) * x;
return d;
}有注释的那一行，交换了 x 和 y ，相当于令 x = $y_0$，y = $x_0$ 。其实本质上没什么区别，只是书写和记忆起来更方便一点。
这样我们就可以求得 ax + by = gcd(a,b) 的一组特解，那么通解该怎么求呢？过程大致如下：
二、拓展gcd应用
一般的题目不会让我们求通解，而会让我们求符合某些特殊条件的解，例如这道：
题目描述
求关于 x 的同余方程 ax ≡ 1 (mod b) 的最小整数解
输入格式
一行，包括两个正整数 a、b ，用一个空格隔开
输出格式
一个正整数 $x_0$ ，即最小正整数解。输入的数据保证一定有解
数据范围
对于 100% 的数据， 2 <= a,b <= 2000000000
所谓 ax ≡ 1 (mod b) ，其实相当于 ax = by + 1 ，于是我们可以用拓展gcd算法求解（要求的是 x ，显然 by 还是 -by 并不重要）。根据上面裴蜀定理的重要推论，ax - by = 1 有解的条件是 gcd(a,b) = 1 ，即 a、b 互质，那么通解就是 $x_k = x + kb$ 。那么已知 x 要求最小整数解，则只需要求 x mod b 即可，但如果 x 是负数，则取模会求到最大负整数解，所以需要加上一个 b 再取模，出于简化书写，可以统一加上 b 然后再模 b 。
#include <bits/stdc++.h>
using namespace std;
// 拓展gcd
int exgcd(int a, int b, int &x, int &y)
{
if (b == 0)
{
x = 1;
y = 0;
return a;
}
int d = exgcd(b, a % b, x, y), x0 = x, y0 = y;
x = y0;
y = x0 - (a / b) * y0;
return d;
}
int main()
{
int a, b, x, y;
scanf("%d%d", &a, &b);
exgcd(a,b,x,y);
printf("%d\n", ((x%b)+b)%b);
return 0;
}上面那题给的数据比较友善，保证了方程一定有解，但是大多数情况下给的不定方程的 a、b 其实是任意的，并不一定保证有解，比如下面这题：
通过题目我们可以看出，如果两只青蛙跳 k 次后会相遇，则存在同余方程 x + km ≡ y + kn (mod L) ，移项可得 x - y + k(m-n) ≡ 0 (mod L) ，即 x - y + k(m-n) = pL ，最后获得等式 k(n-m) +pL = x-y 。
观察等式，我们可以发现跟拓展gcd中的不定方程 ax + by = gcd(a,b) 很类似，k 和 p 即我们需要求的 x、y 。我们可知，对于 ax + by = c ，有解的条件是 c 可以整除 gcd(a,b) ，若式子有解，则通解如下图所示：
前置问题：求解 ax + by = c 的一个解，c 不一定等于 gcd(a,b)
设 g = gcd(a,b) ，则我们可以利用 exgcd 求出方程 ax + by = g 的其中一个解 $x_0$ ，因为 $\frac{ax+by}{g}$ 一定是整数，所以如果 $\frac{c}{g}$ 不是整数，则无解。代入 $x_0$ ，并在在等式两边乘 $\frac{c}{g}$ ，得 $a \frac{cx_0}{g} + b \frac{cy_0}{g} = c$ ，所以我们可知原方程 ax + by = c 的一个解 x = $\frac{cx_0}{g}$ 。由于 gcd 只对非负整数有意义，如果 a < 0 ，等式两边要同时取反， b 本来就是正数，不能变也不用变。
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
// 拓展gcd
ll exgcd(ll a, ll b, ll &x, ll &y)
{
if (b == 0)
{
x = 1;
y = 0;
return a;
}
ll d = exgcd(b, a % b, x, y), x0 = x, y0 = y;
x = y0;
y = x0 - (a / b) * y0;
return d;
}
int main()
{
ll x,y,m,n,L;
scanf("%lld%lld%lld%lld%lld", &x, &y, &m, &n, &L);
ll x0, y0;
ll A = (n-m), c = (x-y);
if(A<0) A = -A, c = -c;
ll gc = exgcd(A,L,x0,y0); // (n-m)k + pL = x-y
if(c%gc != 0) printf("Impossible"); // 无解
else printf("%lld", ((x0*c/gc)%(L/gc)+L/gc)%(L/gc));
return 0;
}除此之外，拓展gcd的典型题型还有求不定方程在一定范围内解的个数、求解模线性方程、求解模的逆元，这些留到下次再展开。
快速GCD
在我学习 gcd 算法的过程中，接触到了一种用位运算来优化普通gcd的方式，这种优化的思路主要是在于尽量消除 % 运算，因为 % 运算会加大时间的花销。此方法分了六种情况讨论：
- 若 x = 0 或 y = 0，返回 y 或 x
- 若x < y ，交换 x , y
- 若 x 为偶数且 y 为偶数，返回 gcd(x/2,x/2) * 2
- 若 x 为偶数且 y 为奇数，返回 gcd(x/2,y)
- 若 x 为奇数且 y 为偶数，返回 gcd(x,y/2)
- 若 x 为奇数且 y 为偶数，返回 gcd(x-y,y)
// 快速gcd
inline int kgcd(int a, int b)
{
if (a == 0) return b;
if (b == 0) return a;
if (!(a & 1) && !(b & 1)) return kgcd(a>>1, b>>1)<<1;
else if (!(a & 1)) return kgcd(a>>1, b);
else if (!(b & 1)) return kgcd(a, b>>1);
else return kgcd(abs(a - b), min(a, b)); // 此处使用abs()和min()实际上已经做了一次交换
}算出来的结果与普通 gcd 并无差别，不过既然说快速，那么还是想看看是不是真的快速，因此我本人不严谨的做了以下的尝试。
第一场
随机生成了9组数据，每组数据里有1000,000个大数对，gcd大多为1，数据生成器如下：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll _rand()
{
return abs((ll)rand()*134982003ll+rand()*32498ll+rand()+1);
}
ll gcd(ll a,ll b)
{
return b == 0 ? a : gcd(b,a%b);
}
int main()
{
int m=1000000;
srand((unsigned)time(NULL));
string fil="gcd";
for(int i=1; i<=9; ++i)
{
printf("Finishing %d...\n",i);
char str[5];
sprintf(str,"%d",i);
FILE* in=fopen((fil+str+".in").c_str(),"w");
FILE* out=fopen((fil+str+".out").c_str(),"w");
fprintf(in,"%d\n",m);
for(int i=1; i<=m; ++i)
{
ll a=_rand(),b=_rand();
fprintf(in,"%lld %lld\n",a,b);
fprintf(out,"%lld\n",gcd(a,b));
}
printf("Complete.\n");
fclose(in);
fclose(out);
}
return 0;
}运行结果如下：普通gcd平均用时 2.006s ，快速gcd平均用时 2.656s，第一场， 普通gcd获胜！！！
第二场
随机生成了9组数据，每组数据里有1000,000个大数对，gcd都较大，数据生成器如下：
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll __rand()
{
return rand()%20+1;
}
ll ___rand()
{
return __rand()*__rand()*__rand()*__rand()*__rand()*__rand();
}
ll _rand()
{
return abs((ll)___rand()*___rand()*___rand());
}
ll gcd(ll a,ll b)
{
return a%b==0?b:gcd(b,a%b);
}
int main()
{
int m=1000000;
srand(time(0));
string fil="gcd";
for(int i=1; i<=9; ++i)
{
printf("Finishing %d...\n",i);
char str[5];
sprintf(str,"%d",i);
FILE* in=fopen((fil+str+".in").c_str(),"w");
FILE* out=fopen((fil+str+".out").c_str(),"w");
fprintf(in,"%d\n",m);
for(int i=1; i<=m; ++i)
{
ll a=_rand(),b=_rand();
fprintf(in,"%lld %lld\n",a,b);
fprintf(out,"%lld\n",gcd(a,b));
}
printf("Complete.\n");
fclose(in);
fclose(out);
}
return 0;
}运行结果如下：普通gcd平均用时 2.185s ，快速gcd平均用时 2.553s，第二场， 普通gcd获胜！！！
可以看到，虽然是叫做快速gcd，但是测试的结果却并不像预料之中那样，我个人猜测原因是，虽然位运算比模运算要快得多，但是在随机数据面前可能并没有太大的优势，所以我们平常使用普通的gcd就可以了
/* 测试代码 */
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
// 普通gcd
inline ll gcd(ll a, ll b)
{
return b == 0 ? a : gcd(b,a%b);
}
// 快速gcd
inline ll kgcd(ll a, ll b)
{
if (a == 0) return b;
if (b == 0) return a;
if (!(a & 1) && !(b & 1)) return kgcd(a>>1, b>>1)<<1;
else if (!(a & 1)) return kgcd(a>>1, b);
else if (!(b & 1)) return kgcd(a, b>>1);
else return kgcd(abs(a - b), min(a, b)); // 此处使用abs()和min()实际上已经做了一次交换
}
ll m,a,b;
int main()
{
string fil="gcd";
double time1 = 0, time2 = 0;
for(int i = 1; i <= 9; ++i)
{
printf((fil+char(i+'0')+".in\n").c_str());
freopen((fil+char(i+'0')+".in").c_str(),"r",stdin);
scanf("%lld", &m);
double s1=clock();
for(int i = 1; i <= m; ++i)
{
scanf("%lld%lld", &a, &b);
gcd(a,b);
}
double t1 = double(clock()-s1)/CLOCKS_PER_SEC;
printf("普通gcd：time used=%.3fs ", t1);
time1 += (t1);
freopen((fil+char(i+'0')+".in").c_str(),"r",stdin);
scanf("%lld", &m);
double s2=clock();
for(int i = 1; i <= m; ++i)
{
scanf("%lld%lld", &a, &b);
kgcd(a,b);
}
double t2 = double(clock()-s2)/CLOCKS_PER_SEC;
printf("快速gcd：time used=%.3fs\n", t2);
time2 += t2;
printf("\n");
}
printf("平均用时\n普通gcd：time used=%.3fs 快速gcd：time used=%.3fs\n", time1/9.0, time2/9.0);
return 0;
} 参考资料
版权属于：PCsky
本文链接：http://www.hyouka.club/index.php/archives/151/
转载时须注明出处及本声明