歐幾里德算法又稱輾轉(zhuǎn)相除法,用于計(jì)算兩個整數(shù)a,b的最大公約數(shù)。其計(jì)算原理依賴于下面的定理:
定理:gcd(a,b) = gcd(b,a mod b)
其算法用C++語言描述為:
int gcd(int m, int n)
{
if (m == 0)
return n;
if (n == 0)
return m;
if (m < n)
{
int tmp = m;
m = n;
n = tmp;
}
while (n != 0)
{
int tmp = m % n;
m = n;
n = tmp;
}
return m;
}
Stein算法(以下理論請參考http://blog.vckbase.com/arong/archive/2004/06/15/458.html),代碼是我加上的。
歐幾里德算法是計(jì)算兩個數(shù)最大公約數(shù)的傳統(tǒng)算法,他無論從理論還是從效率上都是很好的。但是他有一個致命的缺陷,這個缺陷只有在大素?cái)?shù)時才會顯現(xiàn)出來。
考慮現(xiàn)在的硬件平臺,一般整數(shù)最多也就是64位,對于這樣的整數(shù),計(jì)算兩個數(shù)之間的模是很簡單的。對于字長為32位的平臺,計(jì)算兩個不超過32位的整數(shù)的模,只需要一個指令周期,而計(jì)算64位以下的整數(shù)模,也不過幾個周期而已。但是對于更大的素?cái)?shù),這樣的計(jì)算過程就不得不由用戶來設(shè)計(jì),為了計(jì)算兩個超過64位的整數(shù)的模,用戶也許不得不采用類似于多位數(shù)除法手算過程中的試商法,這個過程不但復(fù)雜,而且消耗了很多CPU時間。對于現(xiàn)代密碼算法,要求計(jì)算128位以上的素?cái)?shù)的情況比比皆是,設(shè)計(jì)這樣的程序迫切希望能夠拋棄除法和取模。
Stein算法由J. Stein 1961年提出,這個方法也是計(jì)算兩個數(shù)的最大公約數(shù)。和歐幾里德算法 算法不同的是,Stein算法只有整數(shù)的移位和加減法,這對于程序設(shè)計(jì)者是一個福音。
為了說明Stein算法的正確性,首先必須注意到以下結(jié)論:
gcd(a,a) = a,也就是一個數(shù)和他自身的公約數(shù)是其自身
gcd(ka,kb) = k gcd(a,b),也就是最大公約數(shù)運(yùn)算和倍乘運(yùn)算可以交換,特殊的,當(dāng)k=2時,說明兩個偶數(shù)的最大公約數(shù)必然能被2整除
有了上述規(guī)律就可以給出Stein算法如下:
1.如果A=0,B是最大公約數(shù),算法結(jié)束
2.如果B=0,A是最大公約數(shù),算法結(jié)束
3.設(shè)置A1 = A、B1=B和C1 = 1
4.如果An和Bn都是偶數(shù),則An+1 =An /2,Bn+1 =Bn /2,Cn+1 =Cn *2(注意,乘2只要把整數(shù)左移一位即可,除2只要把整數(shù)右移一位即可)
5.如果An是偶數(shù),Bn不是偶數(shù),則An+1 =An /2,Bn+1 =Bn ,Cn+1 =Cn (很顯然啦,2不是奇數(shù)的約數(shù))
6.如果Bn是偶數(shù),An不是偶數(shù),則Bn+1 =Bn /2,An+1 =An ,Cn+1 =Cn (很顯然啦,2不是奇數(shù)的約數(shù))
7.如果An和Bn都不是偶數(shù),則An+1 =|An -Bn|,Bn+1 =min(An,Bn),Cn+1 =Cn
8.n++,轉(zhuǎn)4
這個算法的原理很顯然,所以就不再證明了。現(xiàn)在考察一下該算法和歐幾里德方法效率上的差別。
考慮歐幾里德算法,最惡劣的情況是,每次迭代a = 2b -1,這樣,迭代后,r= b-1。如果a小于2N,這樣大約需要 4N次迭代。而考慮Stein算法,每次迭代后,顯然AN+1BN+1≤ ANBN/2,最大迭代次數(shù)也不超過4N次。也就是說,迭代次數(shù)幾乎是相等的。但是,需要注意的是,對于大素?cái)?shù),試商法將使每次迭代都更復(fù)雜,因此對于大素?cái)?shù)Stein將更有優(yōu)勢。
其算法用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庫函數(shù)
n = min(m1, n1);//crt宏
}
}
return c * n;
}