求最大公约数的多种算法比较

一.实验内容

运行最大公约数的常用算法，并进行程序的调式与测试，要求程序设计风格良好，并添加异常处理模块（如输入非法等）。

二.题目分析

1.辗转相除法
辗转相除法（又名欧几里德法）C语言中用于计算两个正整数a,b的最大公约数和最小公倍数，实质它依赖于下面的定理：

2.穷举法（利用数学定义）
穷举法（也叫枚举法）穷举法求两个正整数的最大公约数的解题步骤：从两个数中较小数开始由大到小列举，直到找到公约数立即中断列举，得到的公约数便是最大公约数
3. 更相减损法
更相减损术，是出自《九章算术》的一种求最大公约数的算法，它原本是为约分而设计的，但它适用于任何需要求最大公约数的场合。《九章算术》是中国古代的数学专著，其中的“更相减损术”可以用来求两个数的最大公约数，即“可半者半之，不可半者，副置分母、子之数，以少减多，更相减损，求其等也。以等数约之。”
4.Stein算法
Stein算法由J. Stein 1961年提出，这个方法也是计算两个数的最大公约数。来研究一下最大公约数的性质，发现有 gcd( kx,ky ) = kgcd( x,y ) 这么一个非常好的性质。试取 k=2，则有 gcd( 2x,2y ) = 2 * gcd( x,y )。很快联想到将两个偶数化小的方法。那么一奇一个偶以及两个奇数的情况如何化小呢？
先来看看一奇一偶的情况： 设有2x和y两个数，其中y为奇数。因为y的所有约数都是奇数，所以 a = gcd( 2x,y ) 是奇数。根据2x是个偶数不难联想到，a应该是x的约数。我们来证明一下：(2x)%a=0，设2x=na，因为a是奇数，2x是偶数，则必有n是偶数。又因为 x=(n/2)*a，所以 x%a=0，即a是x的约数。因为a也是y的约数，所以a是x和y的公约数，有 gcd( 2x,y ) <= gcd( x,y )。因为gcd( x,y )明显是2x和y的公约数，又有gcd( x,y ) <= gcd( 2x,y )，所以 gcd( 2x,y ) = gcd( x,y )。至此，我们得出了一奇一偶时化小的方法。
再来看看两个奇数的情况：设有两个奇数x和y，不妨设x>y，注意到x+y和x-y是两个偶数，则有 gcd( x+y,x-y ) = 2 * gcd( (x+y)/2,(x-y)/2 )，那么 gcd( x,y ) 与 gcd( x+y,x-y ) 以及 gcd( (x+y)/2,(x-y)/2 ) 之间是不是有某种联系呢？为了方便设 m=(x+y)/2 ，n=(x-y)/2 ，容易发现 m+n=x ，m-n=y 。设 a = gcd( m,n )，则 m%a=0,n%a=0 ，所以 (m+n)%a=0，(m-n)%a=0 ，即 x%a=0 ，y%a=0 ，所以a是x和y的公约数，有 gcd( m,n )<= gcd(x,y)。再设 b = gcd( x,y )肯定为奇数，则 x%b=0,y%b=0 ，所以 (x+y)%b=0 ，(x-y)%b=0 ，又因为x+y和x-y都是偶数，跟前面一奇一偶时证明a是x的约数的方法相同，有 ((x+y)/2)%b=0,((x-y)/2)%b=0 ，即 m%b=0 ，n%b=0 ，所以b是m和n的公约数，有 gcd( x,y ) <= gcd( m,n )。所以 gcd( x,y ) = gcd( m,n ) = gcd( (x+y)/2,(x-y)/2 )。
整理一下，对两个正整数 x>y ：
1.均为偶数 gcd( x,y ) =2gcd( x/2,y/2 )；
2.均为奇数 gcd( x,y ) = gcd( (x+y)/2,(x-y)/2 )；
2.x奇y偶 gcd( x,y ) = gcd( x,y/2 )；
3.x偶y奇 gcd( x,y ) = gcd( x/2,y ) 或 gcd( x,y )=gcd( y,x/2 )；
现在已经有了递归式，还需要再找出一个退化情况。注意到 gcd( x,x ) = x ，就用这个。

三.算法构造

用流程图的形式表示

1.辗转相除法

2.穷举法

3.更相减损法

4.Stein算法

四.算法实现

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <Windows.h>

int Recurse_Euclid(int a, int b);	//辗转相除法-递归调用
int Eulicd(int a, int b);			//辗转相除法-嵌套调用
int Max_divisor(int a, int b);		//穷举法（定义法）求最大公约数
int GCD(int a, int b);				//更相减损法求最大公约数
int Stein_N(int a, int b);			//Stein 非递归
int Stein_T(int u, int v);			//Stein 递归 

int main() {
	int num1, num2;
	clock_t time1, time2;
	int a[40] = { 2,5,6,9,20,15,46,33,90,10000,19,28,700,3198,47,657,2001,60,77,61875,50,80,666,756,807,687,18459,187,889,45,978,4890,7954,809,55,605,4891,315,614,32 };
	
	time1 = clock();
	for (int i = 0; i < 30; i++)
	{
		Recurse_Euclid(a[i], a[i + 1]);
		Sleep(100);
	}
	printf("辗转相除法-递归调用用时：%.6f\n", (double)(clock() - time1) / CLOCKS_PER_SEC-3.000);
	
	time1 = clock();
	for (int i = 0; i < 30; i++)
	{
		Eulicd(a[i], a[i + 1]);
		Sleep(100);
	}
	printf("辗转相除法-嵌套调用用时：%.6f\n", (double)(clock() - time1) / CLOCKS_PER_SEC - 3.000);

	time1 = clock();
	for (int i = 0; i < 30; i++)
	{
		Max_divisor(a[i], a[i + 1]);
		Sleep(100);
	}
	printf("穷举法（定义法）用时：%.6f\n", (double)(clock() - time1) / CLOCKS_PER_SEC - 3.000);

	time1 = clock();
	for (int i = 0; i < 20; i++)
	{
		GCD(a[i], a[i + 1]);
		Sleep(100);
	}
	printf("更相减损法用时：%.6f\n", (double)(clock() - time1) / CLOCKS_PER_SEC - 2.000);

	time1 = clock();
	for (int i = 0; i < 30; i++)
	{
		Stein_N(a[i], a[i + 1]);
		Sleep(100);
	}
	printf("Stein 非递归用时：%.6f\n", (double)(clock() - time1) / CLOCKS_PER_SEC - 3.000);

	time1 = clock();
	for (int i = 0; i < 30; i++)
	{
		Stein_T(a[i], a[i + 1]);
		Sleep(100);
	}
	printf("Stein 递归用时：%.6f\n", (double)(clock() - time1) / CLOCKS_PER_SEC - 3.000);

	return 0;
}

//辗转相除法-递归调用 
int Recurse_Euclid(int a, int b)
{
	int temp;
	if (a < b)	//先保证a>b 
	{
		temp = a;
		a = b;
		b = temp;
	}
	if (a%b == 0) return b;	//递归结束的出口 
	else return(b, a%b);	//若a%b!=0，则将b作为a，a%b作为b递归调用 
}
//辗转相除法-嵌套调用 
int Eulicd(int a, int b)
{
	int temp;
	if (a < b)
	{
		temp = a;
		a = b;
		b = temp;
	}
	while (b != 0)		//将递归的具体过程用循环实现 
	{
		temp = a % b;
		a = b;
		b = temp;
	}
	return a;
}

//穷举法（定义法） 
int Max_divisor(int a, int b)			//求最大公约数 
{
	int temp;
	temp = a>b ? b : a;					//找出两者较小值 
	for (; temp>0; temp--)				//从0遍历到temp，直到a，b均能被temp整除 
	{
		if (a%temp == 0 && b%temp == 0) break;
	}
	return temp;
}

//更相减损法求最大公约数
int GCD(int a, int b)
{
	while (a != b)
	{
		if (a > b)
			a = a - b;
		else
			b = b - a;
	}
	return a;
}
//Stein 非递归 
int Stein_N(int a, int b)
{
	int factor = 0;
	int t;
	if (a < b)
	{
		t = a;
		a = b;
		b = t;
	}
	while (a != b)
	{
		if (a & 0x1) 		//判断a是否为偶数 
		{
			if (b & 0x1) 	//判断b是否为偶数 
			{
				b = (a - b) >> 1;
				a -= b;
			}
			else
			{
				b >>= 1;
			}
		}
		else
		{
			if (b & 0x1)
			{
				a >>= 1;
				if (a < b)
				{
					t = a;
					a = b;
					b = t;
				}
			}
			else
			{
				a >>= 1;
				b >>= 1;
				++factor;
			}
		}
	}
	return (a << factor);
}

//Stein 递归 
int Stein_T(int u, int v)
{
	if (u == 0) return v;
	if (v == 0) return u;
	if (~u & 1)				//判断u是否为偶数 
	{
		if (v & 1)			//判断v是否为0偶数 
			return Stein_T(u >> 1, v);
		else
			return Stein_T(u >> 1, v >> 1) << 1;
	}
	if (~v & 1)
		return Stein_T(u, v >> 1);
	if (u > v)
		return Stein_T((u - v) >> 1, v);
	return Stein_T((v - u) >> 1, u);
}

五.运行结果

六.经验归纳

这次程序编写发现了诸多问题也学到了很多，首先在使用计时函数时，由于运行速度过快，导致所用时间十分接近0，易被忽略，使用Windows中的Sleep函数可以增加每次循环所用时间，让其很小的时间值不被忽略，最后减去Sleep（）所用时间即可得到算法所用时间，但是由于使用的是clock（）函数，精确度还是不够高，仅能精确到毫秒，预计精确到微秒可以较为准确的比较出各个算法之间的用时差距，此部分也是我在程序后期调试中发现的问题，也花费了较多的时间。这是接下来对代码优化的方向。

在编写程序时常会出现思路较为拘谨或者不严谨的地方，如若有更好的想法或者指出错误之处，欢迎评论，一起学习进步！！