1. 程式人生 > >歐幾里德演算法求最大公約數

歐幾里德演算法求最大公約數

歐幾里德演算法又稱輾轉相除法,用於計算兩個整數a,b的最大公約數。

定理:gcd(a,b) = gcd(b,a mod b)。

第一種可以寫成:

int Gcd(int a, int b)
{
    while(b != 0)
    {
        int r = b;
        b = a % b;
        a = r;
    }
    return a;
}

也可寫成下面的形式,本質上一樣的。

int Gcd(int a, int b)
{
    if(b == 0)
        return a;
    return Gcd(b, a % b);
}

 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整除
有了上述規律就可以給出Stein演算法如下:

1.如果A=0,B是最大公約數,演算法結束
2.如果B=0,A是最大公約數,演算法結束
3.設定A1 = A、B1=B和C1 = 1
4.如果An和Bn都是偶數,則An+1 =An /2,Bn+1 =Bn /2,Cn+1 =Cn *2(注意,乘2只要把整數左移一位即可,除2只要把整數右移一位即可)
5.如果An是偶數,Bn不是偶數,則An+1 =An /2,Bn+1 =Bn ,Cn+1 =Cn (很顯然啦,2不是奇數的約數)
6.如果Bn是偶數,An不是偶數,則Bn+1 =Bn /2,An+1 =An ,Cn+1 =Cn (很顯然啦,2不是奇數的約數)
7.如果An和Bn都不是偶數,則An+1 =|An -Bn|,Bn+1 =min(An,Bn),Cn+1 =Cn
8.n++,轉4
這個演算法的原理很顯然,所以就不再證明了。現在考察一下該演算法和歐幾里德方法效率上的差別。

考慮歐幾里德演算法,最惡劣的情況是,每次迭代a = 2b -1,這樣,迭代後,r= b-1。如果a小於2N,這樣大約需要 4N次迭代。而考慮Stein演算法,每次迭代後,顯然AN+1BN+1≤ ANBN/2,最大迭代次數也不超過4N次。也就是說,迭代次數幾乎是相等的。但是,需要注意的是,對於大素數,試商法將使每次迭代都更復雜,因此對於大素數Stein將更有優勢。

判斷n是否為奇數or偶數,奇數返回0,偶數返回1。

bool is_even(int n) 
{
 return !(n & 1);
}

c++演算法程式碼:

bool is_even(int n) 
{
 return !(n & 1);
}
int gcd2(int m, int n)
{
 int c = 1;
 while (m != 0 && n != 0)
 {
  if (is_even(m) && is_even(n))
  {
   m >>= 1;
   n >>= 1;
   c <<= 1;
  }
  else if (is_even(m) && !is_even(n))
  {
   m >>= 1;
  }
  else if (!is_even(m) && is_even(n))
  {
   n >>= 1;
  }
  else if (!is_even(m) && !is_even(n))
  {
   int m1 = m;
   int n1 = n;
   m = abs(m-n);  //crt庫函式
   n = min(m1, n1);//crt巨集
  }
 }

 return c * n;
}

AC code:

#include <iostream>
#include <stdio.h>
using namespace std;
#define min(a,b) ((a)<(b)?(a):(b))  
bool pd(int n)
{
	return !(n&1);
}
int gcd(int n,int m)
{
   int c=1;
   if(n==0 && m)
	   return m;
   if(n && m==0)
	   return n;
   while(m && n && n!=m)
   {
	   if(pd(n) && pd(m))
	   {
		   m>>=1;
		   n>>=1;
		   c<<=1;
	   }
	   else if(pd(n) && !pd(m))
	   {
		   n>>=1;
	   }
	   else if(!pd(n) && pd(m))
	   {
		   m>>=1;
	   }
	   else if(!pd(n) && !pd(m))
	   {
		   int m1=m;
		   int n1=n;
		   m=abs(n-m);
		   n=min(m1,n1);
	   }
   }
       return c*n;
}
int main()
{
	int n,m;
	while(cin>>n>>m)
	{
		cout<<n/gcd(n,m)*m<<endl;
	}
  return 0;
}