最大公约数算法之横向评测

Posted Sun Feb 09 2020. 2897 words. 15 min read.

This article was last updated on days ago, the information described in the article may be outdated.

最大公约数
简称 Gcd
当然还有一个简称 hcf 嘿嘿嘿
(可我不敢用 hcf QAQ)

求法有好多种

  • 穷举法
  • 碾转相除法
    • 欧几里得算法
    • 扩展欧几里得算法
  • Stein 算法
  • 碾转相减法
  • 利用位运算(不知道叫什么,非常快就对了)

算法介绍

1. 穷举法

没什么说的

代码实现

int Gcd1(int a, int b) {
int Gcd = 1, minn = min(a, b);
for (int i=2; i<=minn; i++)
if (!a%i && !b%i)
Gcd = i;
return Gcd;
}

2. 辗转相除法

又称为 欧几里得算法

公式

Gcd(a,  b)=Gcd(b,  a  mod  b)Gcd(a,\;b) = Gcd(b,\;a\;mod\;b)

其中,a mod b是a除以b所得的余数。

证明

  1. 充分性:
    设 g 是 a,b 的公约数,则 a,b 可由 g 来表示:
    a = xg, b = yg (x,y为整数)
    又,a 可由 b 表示(任意一个数都可以由另一个数来表示):
    a = kb + r (k 为整数,r 为 a 除以 b 所得余数)
    => r = a - kb = xg - kyg = (x - ky)g
    即,g 也是 r 的约数。
    这样,g 就是 (b, r) 即 (b, a mod b) 的公约数。

  2. 必要性:
    设g是 (b, a mod b) 的公约数,类似于 1.,同样可以推出 g 也是 (a, b) 的公约数。

综合 1.2.可得 (a, b) 和 (b, a mod b) 的公约数是一样的,当然最大公约数也是一样的。

代码如下:

// 碾转相除(递归写法)
int Gcd2(int a, int b) {
if (!b) return a;
return Gcd2(b, a % b);
}

//碾转相除(迭代写法)
int Gcd3(int a, int b) {
while(b) {
int t = b;
b = a % b;
a = t;
}
return a;
}

补充:
(看大佬博客发现的,以下内容均来自 百度百科大佬博客
扩展欧几里得算法(英语:Extended Euclidean algorithm)是欧几里得算法(又叫辗转相除法)的扩展。
已知整数 a、b,扩展欧几里得算法可以在求得a、b的最大公约数的同时,能找到整数 x、y (其中一个很可能是负数),使它们满足贝祖等式,即

ax+by=Gcd(a,  b)ax + by = Gcd(a,\;b)

如果 a 是负数,可以把问题转化成

a(x)+by=Gcd(a,  b)| a | ( -x ) + by = Gcd(| a |,\;b)

然后令x’ = (-x)。

通常谈到最大公约数时,我们都会提到一个非常基本的事实:给予二个整数 a、b,必存在整数 x、y 使得 ax + by = Gcd(a,b)。
有两个数 a,b,对它们进行辗转相除法,可得它们的最大公约数——这是众所周知的。然后,收集辗转相除法中产生的式子,倒回去,可以得到 ax+by=Gcd (a,b) 的整数解。
扩展欧几里得算法可以用来计算模反元素(也叫模逆元),而模反元素在 RSA 加密算法中有举足轻重的地位。

代码如下:

int exGcd(int a, int b, int &x, int &y) {
if(b == 0) {
x = 1;
y = 0;
return a;
}
int r = exGcd(b, a % b, x, y);
int t = x;
x = y;
y = t - a / b * y;
return r;
}

把这个实现和 Gcd 的递归实现相比,发现多了下面的 x,y 赋值过程,这就是扩展欧几里德算法的精髓。
可以这样思考:
对于

a=b,b=a%ba' = b, b' = a \% b

而言,我们求得 x, y使得

ax+by=Gcd(a,b)a'x + b'y = Gcd(a', b')

由于

b=a%b=aa/bbb' = a \% b = a - a / b * b

(注:这里的/是程序设计语言中的除法)
那么可以得到:

ax+by=Gcd(a,b)bx+(aa/bb)y=Gcd(a,b)=Gcd(a,b)ay+b(xa/by)=Gcd(a,b)a'x + b'y = Gcd(a', b')\\⟹ bx + (a - a / b * b)y = Gcd(a', b') = Gcd(a, b)\\⟹ ay +b(x - a / b*y) = Gcd(a, b)

因此对于 a 和 b 而言,他们的相对应的 p,q 分别是 y 和 (x-a/b*y)

3. Stein 算法

(以下内容参考 大佬博客

欧几里得算法缺陷

欧几里德算法是计算两个数最大公约数的传统算法,他无论从理论还是从效率上都是很好的。但是他有一个致命的缺陷,这个缺陷只有在大素数时才会显现出来。

考虑现在的硬件平台,一般整数最多也就是 64 位,对于这样的整数,计算两个数之间的模是很简单的。对于字长为 32 位的平台,计算两个不超过32位的整数的模,只需要一个指令周期,而计算 64 位以下的整数模,也不过几个周期而已。但是对于更大的素数,这样的计算过程就不得不由用户来设计,为了计算两个超过 64 位的整数的模,用户也许不得不采用类似于多位数除法手算过程中的试商法,这个过程不但复杂,而且消耗了很多 CPU 时间。对于现代密码算法,要求计算 128 位以上的素数的情况比比皆是,设计这样的程序迫切希望能够抛弃除法和取模。 (注:说到抛弃除法和取模,其实辗转相除法可以写成减法的形式)

算法思想

Stein 算法由 J. Stein 1961 年提出,这个方法也是计算两个数的最大公约数。和欧几里德算法 算法不同的是,Stein 算法只有整数的移位和加减法,这对于程序设计者是一个福音。

为了说明 Stein 算法的正确性,首先必须注意到以下结论:

Gcd(a,a) = a,也就是一个数和他自身的公约数是其自身
Gcd(ka,kb) = k Gcd(a,b),也就是最大公约数运算和倍乘运算可以交换,特殊的,当 k=2 时,说明两个偶数的最大公约数必然能被 2 整除。

算法步骤

1、如果 An=Bn,那么 An(或Bn)Cn 是最大公约数,算法结束
2、如果 An=0,Bn 是最大公约数,算法结束
3、如果 Bn=0,An 是最大公约数,算法结束
4、设置 A1=A、B1=B 和 C1=1
5、如果 An 和 Bn 都是偶数,则 An+1=An/2,Bn+1=Bn/2,Cn+1=Cn
2 (注意,乘 2 只要把整数左移一位即可,除 2 只要把整数右移一位即可)
6、如果 An 是偶数,Bn 不是偶数,则 An+1=An/2,Bn+1=Bn,Cn+1=Cn(很显然啦,2 不是奇数的约数)
7、如果 Bn 是偶数,An 不是偶数,则 Bn+1=Bn/2,An+1=An,Cn+1=Cn(很显然啦,2 不是奇数的约数)
8、如果 An 和 Bn 都不是偶数,则 An+1=|An-Bn|/2,Bn+1=min(An,Bn),Cn+1=Cn
9、n 加 1,转 1

C++的实现:

int Gcd4(int a, int b) {
if(a == 0) return b;
if(b == 0) return a;
if(a % 2 == 0 && b % 2 == 0) return 2 * Gcd4(a >> 1, b >> 1);
else if(a % 2 == 0) return Gcd4(a >> 1, b);
else if(b % 2 == 0) return Gcd4(a, b >> 1);
else return Gcd4(abs(a - b), min(a, b));
}

考虑欧几里德算法,最恶劣的情况是,每次迭代 a = 2b -1,这样,迭代后,r= b-1。如果 a 小于 2N,这样大约需要 4N 次迭代。而考虑 Stein 算法,每次迭代后,显然 AN+1BN+1≤ ANBN/2,最大迭代次数也不超过 4N 次。也就是说,迭代次数几乎是相等的。但是,需要注意的是,对于大素数,试商法将使每次迭代都更复杂,因此对于大素数 Stein 将更有优势

4. 辗转相减法

即尼考曼彻斯法

若 a > b,则 a = a - b
若 b > a,则 b = b - a
若 a == b,则 a(或b) 即为最大公约数
若 a != b,则回到 1

证明


a = bq1 + r1 (0 < r1 < b)
b = r1q2 + r2 (0 < r2 < r1)
r1 = r2q3 + r3 (0 < r3 < r2)
……
只要 r1, r2, r3…… 不是 0 就可以继续写下去
我们看到:
b > r1 > r2 > r3 > …… > 0
b 是有限的 r1, r2, r3 是整数
所以至多 b 步后,必有 rn = 0
rn - 2 = rn - 1qn + rn
rn - 1 = rn * qn + 1 + 0
由此可以得到 (a, b) = rn
(来自百度百科)

代码实现

// 辗转相减法(递归写法)
int Gcd5(int a, int b) {
if(a == b) return a;
return a > b ? Gcd5(a-b, b) : Gcd5(a, b-a);
}

// 辗转相减法(迭代写法)
int Gcd6(int a, int b) {
while(a != b) {
if(a > b)
a = a - b;
else
b = b - a;
}
return a;
}

5.利用位运算

不会证明,也没找到证明,直接上代码

代码实现

int Gcd7(int a, int b) {
  while(b^=a^=b^=a%=b);
  return a;
}

是不是超短,但是非常快

运行时间对比

// 测试代码
// 对于同种算法仅测试递归写法
#include<iostream>
#include<windows.h>

// 穷举法
void Gcd1(long long a, long long b) {
long long Gcd = 1, minn = std::min(a, b);
for (long long i=2; i<=minn; i++)
if (!a%i && !b%i)
Gcd = i;
}

// 碾转相除(递归写法)
long long Gcd2(long long a, long long b) {
if (!b) return a;
return Gcd2(b, a % b);
}


long long Gcd4(long long a, long long b) {
if(a == 0) return b;
if(b == 0) return a;
if(a % 2 == 0 && b % 2 == 0) return 2 * Gcd4(a >> 1, b >> 1);
else if(a % 2 == 0) return Gcd4(a >> 1, b);
else if(b % 2 == 0) return Gcd4(a, b >> 1);
else return Gcd4(abs(a - b), std::min(a, b));
}

// 辗转相减法(递归写法)
long long Gcd5(long long a, long long b) {
if(a == b) return a;
return a > b ? Gcd5(a-b, b) : Gcd5(a, b-a);
}

long long Gcd7(long long a, long long b) {
while(b^=a^=b^=a%=b);
return a;
}

int main() {
double time=0;
LARGE_INTEGER nFreq;
LARGE_INTEGER nBeginTime;
LARGE_INTEGER nEndTime;
QueryPerformanceFrequency(&nFreq);
QueryPerformanceCounter(&nBeginTime);

long long a = 77353649739;
long long b = 2147483647;
Gcd1(a, b);

QueryPerformanceCounter(&nEndTime);
time=(double)(nEndTime.QuadPart-nBeginTime.QuadPart)/(double)nFreq.QuadPart;
std::cout<<"运行时间:"<<time*1000<<"ms\n";
return 0;
}

结果:

Gcd1(a, b) 运行时间:14641.3 ms
Gcd2(a, b) 运行时间:0.0007 ms
Gcd4(a, b) 运行时间:0.0051 ms
Gcd5(a, b) 运行时间:0.0057 ms
Gcd6(a, b) 运行时间:0.0017 ms
Gcd7(a, b) 运行时间:0.0007 ms

可见辗转相除法利用位运算的两个算法最快。

那么加大数据,再次比较这两种算法

unsigned long long a = 18446344013706551617;
unsigned long long b = 12416744005372551613;

多次运行结果:

Gcd2(a, b) 运行时间:0.0007~0.0008 ms

Gcd7(a, b) 运行时间:0.0002~0.0005 ms

可见利用位运算的算法是最快的,又短又快 嘿嘿嘿

2.20更新

加大数据,并用 Linux 内置 time 测试时间

数据越加越大

//辗转相除
real 0m0.014s
user 0m0.000s
sys 0m0.016s

//位运算
real 0m0.013s
user 0m0.000s
sys 0m0.016s

好像也没差多少啊?

还是我测的方法不对

测试方法如下:

两个方法分两个文件写

//test1
#include<iostream>
long long Gcd(long long a, long long b) {
if (!b) return a;
return Gcd(b, a % b);
}
int main() {

unsigned long long a = 1844674407370955143;
unsigned long long b = 1844674407370955137;
Gcd(a, b);
return 0;
}
//test2
#include<iostream>
unsigned long long Gcd(unsigned long long a, unsigned long long b) {
while(b^=a^=b^=a%=b);
return a;
}
int main() {
unsigned long long a = 1844674407370955143;
unsigned long long b = 1844674407370955137;
Gcd(a, b);
return 0;
}

然后命令行 g++ test1.cpp -o test1g++ test2.cpp -o test2

然后 time ./test1time ./test2

(应该没错啊)


  1. 部分介绍内容来源百度百科

Comments

Unable to load Disqus, please make sure your network can access.