Stein.J算法求最大公约数的几个实现及分析

首先给出wikipedia上Stein算法的理论依据：

The algorithm reduces the problem of finding the GCD by repeatedly applying these identities:

1. 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.
2. If u and v are both even, then gcd(u, v) = 2·gcd(u/2, v/2), because 2 is a common divisor.
3. 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).
4. 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]
5. Repeat steps 2–4 until u = v, or (one more step) until u = 0. In either case, the GCD is 2^kv, 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,显然第三个算法更多利用了求公约数的原理和性质。怎麽说呢，跳！

#include <stdio.h>
#include <stdlib.h>
typedef unsigned long long uint64;

uint64 get_greatest_common_divisor(uint64 a,uint64 b)
{
    uint64 x=a;
    uint64 y=b;
    int k=1;

    while(x!=y)
    {
        if(y==0)
        {
            return x*k;
        }
        else if(x==0)
        {
            return y*k;
        }

        if(y%2==0 && x%2==0)
        {
            x=x/2;
            y=y/2;
            k=k*2;
        }
        else if(x%2!=0 && y%2!=0)
        {
            if(x>y)
            {
                x=(x-y)/2;
            }
            else
            {
                y=(y-x)/2;
            }
        }
        else if(x%2==0)
        {
            x=x/2;
        }
        else if(y%2==0)
        {
            y=y/2;
        }
    }

    return x*k;
}

uint64 get_greatest_common_divisor2(uint64 a,uint64 b)
{
    uint64 x=a;
    uint64 y=b;
    int k=1;
    int d=0;
    while(x!=y)
    {
        if(y==0)
        {
            return x*k;
        }
        else if(x==0)
        {
            return y*k;
        }

        if((x&1)==0 && (y&1)==0)
        {
            x=x>>1;
            y=y>>1;
            k=k<<1;
        }
        else if((x&1)!=0 && (y&1)!=0)
        {
            if(x>y)
            {
                x=(x-y)>>1;
            }
            else
            {
                y=(y-x)>>1;
            }
        }
        else if((x&1)==0)
        {
            x=x>>1;
        }
        else if((y&1)==0)
        {
            y=y>>1;
        }
    }

    return x*k;
}

/*
given by wikipedia
http://en.wikipedia.org/wiki/Binary_GCD_algorithm
*/

uint64 gcd(uint64 u, uint64 v)
{
 int shift;

 /* GCD(0,x) := x */
 if (u == 0 || v == 0)
   return u | v;

 /* Let shift := lg K, where K is the greatest power of 2
    dividing both u and v. */
 for (shift = 0; ((u | v) & 1) == 0; ++shift) {
     u >>= 1;
     v >>= 1;
 }

 while ((u & 1) == 0)
   u >>= 1;

 /* From here on, u is always odd. */
 do {
     while ((v & 1) == 0)  /* Loop X */
       v >>= 1;

  /* Now u and v are both odd. Swap if necessary so u <= v,
     then set v = v - u (which is even). For bignums, the
     swapping is just pointer movement, and the subtraction
     can be done in-place. */
  if (u > v) {
     uint64 t = v; v = u; u = t;}  // Swap u and v.
  v = v - u;                       // Here v >= u.
 } while (v != 0);

 return u << shift;
}

 

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的运行步数要比前面两个多一步？

这个留作下回分解吧。