首先给出wikipedia上Stein算法的理论依据:
The algorithm reduces the problem of finding the GCD by repeatedly applying these identities:
- gcd(0, v) = v, because everything divides zero, and v is the largest number that divides v. Similarly, gcd(u, 0) = u. gcd(0, 0) is not typically defined, but it is convenient to set gcd(0, 0) = 0.
- If u and v are both even, then gcd(u, v) = 2·gcd(u/2, v/2), because 2 is a common divisor.
- If u is even and v is odd, then gcd(u, v) = gcd(u/2, v), because 2 is not a common divisor. Similarly, if u is odd and v is even, then gcd(u, v) = gcd(u, v/2).
- If u and v are both odd, and u ≥ v, then gcd(u, v) = gcd((u − v)/2, v). If both are odd and u < v, then gcd(u, v) = gcd((v − u)/2, u). These are combinations of one step of the simple Euclidean algorithm, which uses subtraction at each step, and an application of step 3 above. The division by 2 results in an integer because the difference of two odd numbers is even.[3]
- Repeat steps 2–4 until u = v, or (one more step) until u = 0. In either case, the GCD is 2kv, where k is the number of common factors of 2 found in step 2.
所以相比较欧几里德算法,该算法在二进制的计算机上可以更多发挥位操作的威力,而避免了模运算带来的处理器性能损失。wikipedia上说This is particularly critical on embedded platforms that have no direct processor support for division。不过我想,今天的智能手机估计不存在这麽弱的处理器,但如果还是比较低端的工控芯片...but but,谁会在那种设备上运行大数的公约数计算呢。
依据Stein算法的理论描述,这里给出三个算法。
第一个函数并没有利用位运算操作符,但是思想还是Stein算法的思想,第二个函数在第一个函数的基础上把运算改成位运算。这两个函数是自己实现的,再来看看第三个函数摘自wikipedia,显然第三个算法更多利用了求公约数的原理和性质。怎麽说呢,跳!
1: #include <stdio.h>
2: #include <stdlib.h>
3: typedef unsigned long long uint64;
4:
5: uint64 get_greatest_common_divisor(uint64 a,uint64 b)
6: {
7: uint64 x=a;
8: uint64 y=b;
9: int k=1;
10:
11: while(x!=y)
12: {
13: if(y==0)
14: {
15: return x*k;
16: }
17: else if(x==0)
18: {
19: return y*k;
20: }
21:
22: if(y%2==0 && x%2==0)
23: {
24: x=x/2;
25: y=y/2;
26: k=k*2;
27: }
28: else if(x%2!=0 && y%2!=0)
29: {
30: if(x>y)
31: {
32: x=(x-y)/2;
33: }
34: else
35: {
36: y=(y-x)/2;
37: }
38: }
39: else if(x%2==0)
40: {
41: x=x/2;
42: }
43: else if(y%2==0)
44: {
45: y=y/2;
46: }
47: }
48:
49: return x*k;
50: }
51:
52: uint64 get_greatest_common_divisor2(uint64 a,uint64 b)
53: {
54: uint64 x=a;
55: uint64 y=b;
56: int k=1;
57: int d=0;
58: while(x!=y)
59: {
60: if(y==0)
61: {
62: return x*k;
63: }
64: else if(x==0)
65: {
66: return y*k;
67: }
68:
69: if((x&1)==0 && (y&1)==0)
70: {
71: x=x>>1;
72: y=y>>1;
73: k=k<<1;
74: }
75: else if((x&1)!=0 && (y&1)!=0)
76: {
77: if(x>y)
78: {
79: x=(x-y)>>1;
80: }
81: else
82: {
83: y=(y-x)>>1;
84: }
85: }
86: else if((x&1)==0)
87: {
88: x=x>>1;
89: }
90: else if((y&1)==0)
91: {
92: y=y>>1;
93: }
94: }
95:
96: return x*k;
97: }
98:
99: /*
100: given by wikipedia
101: http://en.wikipedia.org/wiki/Binary_GCD_algorithm
102: */
103:
104: uint64 gcd(uint64 u, uint64 v)
105: {
106: int shift;
107:
108: /* GCD(0,x) := x */
109: if (u == 0 || v == 0)
110: return u | v;
111:
112: /* Let shift := lg K, where K is the greatest power of 2
113: dividing both u and v. */
114: for (shift = 0; ((u | v) & 1) == 0; ++shift) {
115: u >>= 1;
116: v >>= 1;
117: }
118:
119: while ((u & 1) == 0)
120: u >>= 1;
121:
122: /* From here on, u is always odd. */
123: do {
124: while ((v & 1) == 0) /* Loop X */
125: v >>= 1;
126:
127: /* Now u and v are both odd. Swap if necessary so u <= v,
128: then set v = v - u (which is even). For bignums, the
129: swapping is just pointer movement, and the subtraction
130: can be done in-place. */
131: if (u > v) {
132: uint64 t = v; v = u; u = t;} // Swap u and v.
133: v = v - u; // Here v >= u.
134: } while (v != 0);
135:
136: return u << shift;
137: }
114~117算出u,v的因数2的数量,也即是所要求的2的最大幂。该循环将停止在u或者v不能被2整除也即不在含有2的因数或不再是偶数为止,很容易推理通过这一次简单的循环就可以求出u,v含有公因数2的个数。
经过114~117,u,v将至多存一个偶数项或者都为奇数。
如果u为偶数,119到120行将u减半直至为奇数。然后进入123的循环。
123~134循环始终保证u<=v,并不断进行算法第4步的操作也即是v=v-u,v=v/2直至奇数(124~125行)。最终循环停止在v==0上。
可以看出最后一个算法的实现中,循环停止的条件并没有用到u==v这个条件,而是让v减到0,所以实现3的运行步数要比前面两个多一步?
这个留作下回分解吧。