如果有一个自然数a能被自然数b整除，则称a为b的倍数，b为a的约数。几个自然数公有的约数，叫做这几个自然数的公约数。公约数中最大的一个公约数，称为这几个自然数的最大公约数（Greatest Common Divisor，简写为GCD）。例如，自然数12和30的公约数有1、2、3、6，其中6就是12和30的最大公约数。
求最大公约数有多种方法，常见的有辗转相除法、短除法、更相减损法、质因数分解法等。
1．辗转相除法。
辗转相除法的基本做法是用较大的数除以较小的数，若余数为0，则其中较小的数（即除数）就是最大公约数；若余数不为0，则用前面较小的除数和得出的余数构成新的一对数，继续做上面的除法，直到出现能够整除的两个数（余数为0）。
下面以求48和27的最大公约数为例，求解过程如下：
48 ÷ 27 = 1余 21
27 ÷ 21 = 1余6
21 ÷ 6 =3余3
6 ÷ 3=2
所以3就是48和27的最大公约数。
辗转相除法使用到的原理也很简单，假设用f（x, y）表示x，y的最大公约数，取k = x/y，b = x%y，则x = ky + b，如果一个数能够同时整除x和y，则必能同时整除b和y；而能够同时整除b和y的数也必能同时整除x和y，即x和y的公约数与b和y的公约数是相同的，其最大公约数也是相同的，则有f（x, y）= f（y, x%y）（y > 0），如此便可把原问题转化为求两个更小数的最大公约数，直到其中一个数为0，剩下的另外一个数就是两者最大的公约数。
采用辗转相除法编写的求两个整数m和n的最大公约数的函数及测试主程序如下。
#includeint gcd(int m,int n){ int r=m%n; while (r!=0) { m = n; n = r; r = m%n; } return n;}int main(){ int m,n,t; while (scanf("%d %d",&m,&n)!=EOF) { t = gcd(m,n); printf("%d\n",t); } return 0;}
这个函数也可以写成递归的形式，如下。
int gcd(int m,int n){ if (m%n==0) return n; return gcd(n,m%n);}
2．短除法。
短除法是不断地找两个数的最小公约数，直到找不到为止，然后把所有的最小公约数乘起来，就是最大公约数。
下面以求48和36的最大公约数为例，操作如下：
从除1外的最小的约数2开始试探起
48 ÷ 2 = 24 36 ÷ 2 =18 ，2是约数
24 ÷ 2 = 12 18 ÷ 2 =9 ， 2还是约数
12 ÷ 3 = 4 9 ÷ 3 =3 ， 3是约数
4和3不存在公约数，结束。
所以，2×2×3=12就是48和36的最大公约数。
采用短除法编写的求两个整数m和n的最大公约数的函数如下。
int gcd(int m,int n){ int i,r = 1; for (i=2;i<=m && i<=n;i++) { while (m%i==0&&n%i==0) { r = r*i; m = m/i; n = n/i; } } return r;}
3．更相减损法。
更相减损术出自于中国古代的《九章算术》，也是一种求最大公约数的算法，其基本方法是：
① 先判断两个数的大小，如果两数相等，则这个数本身就是就是它的最大公约数。
② 如果不相等，则用大数减去小数，然后用这个较小数与它们相减的结果相比较，如果相等，则这个差就是它们的最大公约数；如果不相等，则继续执行②操作。
例如，求24与15的最大公约数的逐步减损过程如下：
(24,15) -> (9,15) -> (9,6) -> (3,6) -> (3,3)
所以，3就是24和15的最大公约数。
采用更相减损法编写的求两个整数m和n的最大公约数的函数如下。
int gcd(int m,int n){ int t; while (m!=n) { if (m<n) {t=m; m=n; n=t; } m=m-n; } return n;}
4．Stein算法。
Stein算法由J. Stein 1961年提出，这个方法在计算两个数的最大公约数时，只有整数的移位和加减法，对于程序来说，运算方便效率高。
Stein算法的基本原理是：
1）若a和b都是偶数，则记录下公约数2，然后都除以2（右移1位）；
2）若其中一个数是偶数，则偶数除以2，因为此时2不可能是这两个数的公约数了；
3）若两个都是奇数，则a = |a-b|，b = min(a,b)，因为若d是a和b的公约数，那么d也是|a-b|和min(a,b)的公约数。
设gcb(a,b)表示求整数a、b的最大公约数，则这个算法可描述为
当a和b均为偶数时，gcb(a,b) = 2gcb(a/2, b/2) = 2 * gcb(a>>1, b>>1)；
当a为偶数，b为奇数时，gcb(a,b) = gcb(a/2, b) = gcb(a>>1, b)
当a为奇数，b为偶数时，gcb(a,b) = gcb(a, b/2) = gcb(a, b>>1)
当a和b均为奇数，利用更相减损术运算一次，gcb(a,b) = gcb(b, a-b)，此时a-b的结果必然是偶数，又可以继续进行移位运算。
采用Stein算法编写的求两个整数m和n的最大公约数的函数如下。
int gcd(int m,int n){ int p=1,t; while (m%n!=0) { if (m<n) { t=m; m=n; n=t; } if (m%2==0) { if (n%2==0) { p=p*2; m=m>>1; n=n>>1; } else { m=m>>1; } } else { if (n%2==0) { n=n>>1; } else { m=m-n; m=m>>1; } } } return p*n;}
这一算法可以写成递归的形式，如下。
int gcd(int m,int n){ if (m < n) { int tmp = m; m = n; n = tmp; } if (m%n==0) return n; if (m % 2 == 0 && n % 2 == 0) return 2*gcd(m >> 1, n >> 1); else if (m%2 == 0 && n%2 != 0) return gcd(m >> 1, n); else if (m%2 != 0 && n % 2 == 0) return gcd(m, n >> 1); else if (m % 2 != 0 && n % 2 != 0) return gcd(n, (m - n) >> 1);}
【例1】最大最小公倍数。
问题描述
已知一个正整数N，问从1~N中任选出三个数，他们的最小公倍数最大可以为多少。
输入
输入一个正整数N（1 <= N <= 106）。
输出
输出一个整数，表示你找到的最小公倍数。
输入样例
9
输出样例
504
（1）编程思路。
设有整数x和y，若x和y的最大公约数为k，则它们的最小公倍数为x*y/k。显然，若k=1，则x和y的最小公倍数才尽可能大。因此，要在1~N中选尽可能大的3个互质（最大公约数为1）的数，它们的最小公倍数才最大。
显然，x与x-1是互质的，x-1与x-2也是互质的。当然 x 与 x-2就不一定互质了。若 x是偶数，它们就存在约数2。若x是奇数，x与x-2也是互质的。
因此，在1~N中要选3个数的话，若N是奇数，就直接选 N、N-1和N-2这3个数，它们彼此互质，最小公倍数就是N*(N-1)*(N-2)，为最大。
若N是偶数，则先选N、N-1和N-3这3个数，若这3个数互质，则满足要求；若不互质，则选N-1、N-2和N-3这3个数即可。
（2）源程序。
#include int gcd(int a,int b){ if (a%b==0) return b; else return gcd(b,a%b);}int main(){ int n,n1,n2,n3; scanf("%d",&n); if (n%2==1) { n1=n; n2=n-1; n3=n-2; } else { n1=n; n2=n-1; n3=n-3; while (gcd(n1,n3)!=1 || gcd(n2,n3)!=1) n3--; if (n3!=n-3) { n1=n-1; n2=n-2; n3=n-3; } } long long ans; ans=1ll*n1*n2*n3; printf("%lld\n",ans); return 0;}
【例2】最大公约数和最小公倍数问题。
问题描述
输入二个正整数x0,y0(2<=x0<=100000，2<=y0<=100000),求出满足下列条件的P，Q的个数
条件: 1.P,Q是正整数
2.要求P,Q以x0为最大公约数,以y0为最小公倍数.
试求:满足条件的所有可能的两个正整数的个数。
输入
两个数x0,y0
输出
一个数表示答案
输入样例
3 60
输出样例
4
（1）编程思路。
显然P和Q的最小值为x0，最大值为y0，直接对P进行穷举即可。并且穷举时，P的值的变化是 x0、2*x0、3*x0、…直到y0。
（2）源程序。
#include int gcd(int m, int n){ int r; while(m%n!=0) { r=m%n; m = n; n = r; } return n;}int main(){ int x,y,i,cnt=0; scanf("%d%d",&x,&y); for(i=x;i<=y;i+=x) // 因为所求的两个数都是x的整数倍，所以每次加x if ((x*y)%i==0&& gcd(x*y/i,i)==x) // ∵i*j=x*y ∴j=x*y/i ∴只要枚举i既可 cnt++; printf("%d\n",cnt); return 0;}
将上面的源程序提交给洛谷题库P1029 [NOIP2001 普及组] 最大公约数和最小公倍数问题 （https://www.luogu.com.cn/problem/P1029），可以Accepted。
【例3】第k个与M互质的数。
问题描述
如果两个整数的最大公因数（GCD）为1，则两个正整数被称为彼此互质的。例如，1、3、5、7、9……都与2006是互质的。
对于给定的整数m，当把与m互质的数按升序排列后，找到第K个数，它与m是互质的。
输入
输入包含多个测试用例。对于每个测试用例，它包含两个整数m（1<=m<=1000000）、K（1<=K<=100000000）。
输出
在单行中输出第K个元素。
输入样例
2006 1
2006 2
2006 3
输出样例
1
3
5
（1）编程思路。
因为 gcd(a+b*k,b)=gcd(b,a%b)， gcd(a,b)=gcd(b,a%b)，
所以 gcd(a+b*k,b)=gcd(a, b)， k为常数。
这表明了：对于与b互质的数，它们对b取模的余数会周期性出现。
故只需要计算出在b的范围内， 与b互质的数有哪些。然后看第k个与b互质的数是在第几个周期的第几个就可以了。
（2）源程序。
#includeint s[1000005];int gcd(int a,int b){ int r; while (a%b!=0) { r=a%b; a=b; b=r; } return b;}int main(){ int m,k,i,num; while (scanf("%d%d",&m,&k)!=EOF) { num=0; for(i=1;i<=m;i++) { if(gcd(m,i)==1) { s[num++]=i; } } if (k%num==0) { printf("%d\n",(k/num-1)*m + s[num-1]); } else { printf("%d\n",k/num*m + s[k%num-1]); } } return 0;}
将上面的源程序提交给北大POJ题库 POJ 2773 Happy 2006 （http://poj.org/problem?id=2773），可以Accepted。
【例4】可见的点。
问题描述
对于除原点（0，0）外的第一象限的任意一个格点（x，y）（x和y是大于或等于0的整数），如果从（0，0）到（x，y）的连线没有穿过任何其它的格点，则称格点（x，y）是从原点可见的。
例如，点（4，2）不可见，因为来自原点的线穿过了（2，1），如下图所示。也就是说前面的点（2，1）可以挡住后面的点（4，2）。
编写一个程序，在给定大小N的值后，计算可见点（x，y）的数量，其中0≤ x、y≤N。
输入
第一行输入包含单个整数C（1≤C≤1000），表示测试用例的数量。
每个测试用例由包含单个整数N（1≤ N≤ 1000）的单行输入组成。
输出
对于每个测试用例，都有一行输出，包括：从1开始的测试用例编号、N的大小以及该大小下可见点的数量，数据间用空格分隔。
输入样例
4
2
4
5
231
输出样例
1 2 5
2 4 13
3 5 21
4 231 32549
（1）编程思路。
本题实际上就是要求在这个N*N的网格上从原点到各格点的直线有多少种不同的斜率。
边上的两点先不管，然后将整个正方形分成上三角和下三角两部分，由于对称性，两边可以看到的点的数目肯定一样多。以下三角为例进行说明。
1×1只有一个斜率为0的。
2×2 斜率有0，1/2（0已经算过了，以后不再算了），其实就多了一个斜率为1/2的。
3×3 斜率有1/3，2/3两个，比以前多了2个。
4×4 斜率有1/4，2/4（1/2已经有过了），3/4，所以也是2个。
5×5 斜率有1/5，2/5，3/5，4/5，之前都没有，所以多了4个
6×6 斜率有1/6，2/6（1/3已经有了），3/6(1/2已经有了)，4/6(2/3已经有了)，5/6，所以只多2个。
更一般地 n×n，可以从0,0连接到（n,0）到(n,n)上，斜率将会是1/n，2/n，……，(n-1)/n。
由此发现，对于所有能看到的点，它们有着一个共同的特征，
那就是gcd(x,y）=1，若不为1，则它前面肯定有一个点挡住了这个点。
由此问题就转换为：对于一个数n，求小于n的与n互质的数的个数，这就是欧拉函数。
欧拉函数的定义: E(k)=([1,n-1]中与n互质的整数个数)。
因为任意正整数k都可以唯一表示成如下形式：
k=p1^a1*p2^a2*……*pi^ai;(即分解质因数形式)
可以推出： E(k)=(p1-1)(p2-1)……(pi-1)*(p1^(a1-1))(p2^(a2-1))……(pi^(ai-1))
=k*(p1-1)(p2-1)……(pi-1)/(p1*p2*……pi);
=k*(1-1/p1)*(1-1/p2)….(1-1/pk)
利用欧拉函数如下性质，可以快速求出欧拉函数的值(a为N的质因子)。
若 (N%a==0 && (N/a)%a==0) 则有 E(N)=E(N/a)*a;
若 (N%a==0 && (N/a)%a!=0) 则有 E(N)=E(N/a)*(a-1);
计算好下三角后，最终结果=下三角的个数×2 再加1（斜率为1的点）。
（2）源程序。
#include #include <string.h>#define MAXN 1001int main(){ int c,n,i,j,sum; int prime[MAXN], pNum, phi[MAXN]; bool vis[MAXN]; memset(vis,true,sizeof(vis)); // 下面程序段既求MAXN以内的质数又求欧拉数。 pNum=0; phi[1] = 1; for (i = 2; i < MAXN; i++) { if (vis[i]) // i是质数 { prime[pNum++] = i; phi[i] = i-1; } for (j = 0; j < pNum && prime[j]*i <MAXN; j++ ) { vis[prime[j]*i] = false; if (i % prime[j] == 0) { phi[i*prime[j]] = phi[i] * prime[j]; break; } else phi[i*prime[j]] = phi[i] *(prime[j] - 1); } } scanf("%d", &c); for (i =1; i <=c; i++) { scanf("%d", &n); sum = 0; for (j = 1; j <= n ; j++) { sum += phi[j]; } printf("%d %d %d\n", i, n, 2 * sum + 1); } return 0;}
将上面的源程序提交给北大POJ题库 POJ 3090 Visible Lattice Points （http://poj.org/problem?id=3090），可以Accepted。
【例5】最大公约数之和。
问题描述
给定一个整数N（1<N<2^31），计算∑gcd（i，N） 1<=i<=N。
输入
输入包含多个测试用例。
每行一个整数N。
输出
对于每个N个输出1-N中所有的数与N的最大公约数的和。
输入样例
2
6
输出样例
3
15
（1）编程思路。
设x与n的最大公约数是a，也就是gcd(x，n)=a，那么我们就可以知道gcd（x/a，n/a）=1，也就是x/a与n/a互质。
这说明在1-n中与 n 的最大公约数为a的数都与n/a互质，所以我们可以用欧拉函数求出1-n/a中，与n/a互质的数有多少个，其中，a是n的约数，所以遍历所有n的约数。
本题要求的是1~n中各数与n的最大公约数的和sum，所以可以枚举N的所有约数，当a是n的约数时，和值sum就加上a*euler（n/a）。
其中， euler(x) 是欧拉函数，其返回值是 1 到 x 中有多少个数字与 x 互质。
（2）源程序。
#includelong long euler(long long x){ long long i,ret=x; for (i=2;i*i<=x;i++) { if (x%i==0) ret=ret/i*(i-1); while (x%i==0) x/=i; } if (x>1) ret=ret/x*(x-1); return ret;}int main(){ long long n; while (scanf("%lld",&n)!=EOF) { long long i,ans=0; for (i=1;i*i<=n;i++) { if (n%i==0) { ans+=euler(n/i)*i; if (i*i!=n) { ans+=euler(i)*(n/i); } } } printf("%lld\n",ans); } return 0;}
将上面的源程序提交给北大POJ题库 POJ 2480 Longge’s problem （http://poj.org/problem?id=2480），可以Accepted。
【例6】行星共线。
问题描述
X星的行星系统中有n颗行星，它们以位于同一平面的圆形轨道围绕X星运行。每个行星的角速度是恒定的。所有行星的轨道方向都是相同的。
编写程序求连续两次所有行星和恒星X位于同一直线上的时间间隔长度。
输入
输入文件的第一行包含n-行星数（2≤ n≤ 1000)。
第二行包含n个整数ti-行星的轨道周期（1≤ti≤10000)，并非所有的ti都是相同的。
输出
将答案输出为不可约分数，用空格分隔分子和分母。
输入样例
3
6 2 3
输出样例
3 1
（1）编程思路。
任意两颗行星共线，其运行距离差为半个周长的整数倍即可。
已知每个行星的角速度为vi = 2*π/ti，选择一个行星T0作为坐标系，则其他行星的相对速度为vi’ = (2*π/t0-2*π/ti) = (t0 – ti)*2π/(t0*ti)，则角度绕过半个圆周的时间为Ti’ = π/vi’ = (t0*ti)/((t0 – ti)*2)。
所求的答案就是求所有Ti‘的分子的最小公倍数LCM和所有Ti’分母的最大公约数GCD。
由于题目给定2≤n≤1000，1≤ti≤10000，因此最小公倍数会超出长整数表数范围，需要采用高精度计算。在求最小公倍数时，采用将最小公倍数分解质因数的方法，即将LCM分解为Π(pi^bi)的形式，最后根据这个分解式计算出所有Ti的最小公倍数即可。
另外，在求Ti’=(t0*ti)/((t0 – ti)*2)的时候先对分子和分母进行约分，这样就能够保证最后分子分母互质。
（2）源程序。
#include #include <string.h>int prime[5000],pcnt=0;void init() // 生成10000以内的质数表{ int i,j; int vis[10005]={0}; for (i=2;i<10000;i++) { if (!vis[i]) { prime[pcnt++]=i; for (j=2*i;j<10000;j+=i) vis[j]=1; } }}int gcd(int a,int b){ if (b==0) return a; return gcd(b,a%b);}int lcm(int a,int b){ return a*b/gcd(a,b);}int abs(int x){ return x>0?x:-x;}int main(){ init(); int n; while (scanf("%d",&n)!=EOF) { int fm=0; int fz[1001]; // 用于保存分子（各数的最大公约数），是一个大整数，每个元素保存4位数字 int p[5001]; // 用于保存各数最小公倍数的质因子的指数 memset(p,0,sizeof(p)); memset(fz,0,sizeof(fz)); int i,j,t[1001]; for (i=0;i<n;i++) scanf("%d",&t[i]); for (i=1;i<n;i++) { if (t[i]==t[0]) continue; int tmp1=lcm(t[i],t[0]); int tmp2=2*abs(tmp1/t[i]-tmp1/t[0]); int temp=gcd(tmp1,tmp2); tmp1/=temp; // 分子与分母约分，保证两两互质 tmp2/=temp; fm=gcd(fm,tmp2); // 求分母的最大公约数 for (j=0;j<pcnt;j++) // 将分子的最小公倍数分解质因数 { int cnt=0; while (tmp1%prime[j]==0) { cnt++; tmp1/=prime[j]; } if (p[j]<cnt) p[j]=cnt; // 保存质因数prime[j]较大的指数 if (tmp1==1)break; } } fz[0]=1; // 按各质因子及指数求得分子 int len=1; // 分子大整数的位数（按万进制） for (i=0;i<pcnt;i++) { for (j=0;j<p[i];j++) { int k,cf=0; for (k=0;k<len;k++) { fz[k]=fz[k]*prime[i]+cf; cf=fz[k]/10000; fz[k]%=10000; } if (cf!=0) fz[len++]=cf; } } printf("%d",fz[len-1]); for (i=len-2;i>=0;i--) printf("%04d",fz[i]); printf(" %d\n",fm); } return 0;}
将上面的源程序提交给北大POJ题库POJ 3101Astronomy （http://poj.org/problem?id=3101） ，可以Accepted。
【练习1】POJ1284 Primitive Roots （http://poj.org/problem?id=1284）。
#include int getPhiValue(int n) // 求n的欧拉函数值phi(n){ int i,sum,t; sum = n; t=n; for (i=2;i*i<=t;i++) if (t % i == 0) // 找到一个质因数i { sum -= sum/i; while (t % i == 0) // 将质因数i全去掉 t /= i; } if (t!=1) sum -= sum/t; return sum;}int main(){ int p; while (scanf("%d", &p) !=EOF) { printf("%d\n",getPhiValue(p-1)); } return 0;}
参考程序1
#include #define MAXN 65536int phi[MAXN];int prime[MAXN], pNum;int vis[MAXN]={0};void getPhiTable(int n) // 求2~n的欧拉函数值phi[i] （1<i<=n){ int i,j; pNum=0; phi[1] = 1; for (i = 2; i < MAXN; i++) { if (vis[i]==0) // i是素数 { prime[pNum++] = i; phi[i] = i-1; } for (j = 0; j < pNum && prime[j]*i <MAXN; j++ ) { vis[prime[j]*i] = 1; if (i % prime[j] == 0) { phi[i*prime[j]] = phi[i] * prime[j]; break; } else phi[i*prime[j]] = phi[i] *(prime[j] - 1); } }}int main(){ int p; getPhiTable(MAXN); while (scanf("%d", &p) !=EOF) { printf("%d\n",phi[p-1]); } return 0;}
参考程序2
【练习2】POJ 1580 String Matching （http://poj.org/problem?id=1580）。
// 采用穷举法搜索两串匹配的首字母，不断更新最大匹配值即可，// 最后的答案需要对最大匹配值和两串总长进行约分 #include #include <string.h> int gcd(int x,int y){ int r; while (x%y!=0) { r=x%y; x=y; y=r; } return y;}int min(int x,int y){ return x<y?x:y;}int main(){ int i,j,k,len1,len2,a,b,res,max; char c1[1000],c2[1000]; while(1) { scanf("%s",c1); if(strcmp(c1,"-1")==0) break; scanf("%s",c2); len1=strlen(c1); len2=strlen(c2); max=0; for(i=0;i<len1;i++) // 穷举两字符串首字母 for(j=0;j<len2;j++) { res=0; for(k=0;k<min(len1-i,len2-j);k++) if(c1[k+i]==c2[k+j]) res++; if (max<res) max=res; } if(max==0||2*max-len1-len2==0) printf("appx(%s,%s) = %d\n",c1,c2,2*max/(len1+len2)); else { a=2*max; b=len1+len2; i=gcd(a,b); a/=i; b/=i; printf("appx(%s,%s) = %d/%d\n",c1,c2,a,b); } } return 0;}
参考程序
【练习3】POJ 2478 Farey Sequence （http://poj.org/problem?id=2478）。
#include #include <string.h>#define MAXN 1000005int prime[MAXN], pNum, phi[MAXN];long long num[MAXN]={0};int vis[MAXN];int main(){ int n,i,j; memset(vis,1,sizeof(vis)); pNum=0; phi[1] = 1; for (i = 2; i < MAXN; i++) { if (vis[i]) // i是素数 { prime[pNum++] = i; phi[i] = i-1; } for (j = 0; j < pNum && prime[j]*i <MAXN; j++ ) { vis[prime[j]*i] = 0; if (i % prime[j] == 0) { phi[i*prime[j]] = phi[i] * prime[j]; break; } else phi[i*prime[j]] = phi[i] *(prime[j] - 1); } } for (i=2;i<MAXN;i++) { num[i]=num[i-1]+phi[i]; } while (scanf("%d", &n) && n!=0) { printf("%lld\n",num[n]); } return 0;}
参考程序
【练习4】POJ 3002 Jackpot （http://poj.org/problem?id=3002）。
#includelong long gcd(long long a,long long b){ if (a%b==0) return b; return gcd(b,a%b);}long long lcm(long long a,long long b){ return a/gcd(a,b)*b;}int main(){ int t; scanf("%d",&t); while(t--) { int i,n; scanf("%d",&n); long long a,b; scanf("%lld",&a); for (i=1;i<n;i++) { scanf("%lld",&b); a=lcm(a,b); } if (a<=1000000000) printf("%lld\n",a); else printf("More than a billion.\n"); } return 0;}
参考程序