Stein算法

Stein算法是一種計算兩個數最大公約數的算法，是針對歐幾里德算法在對大整數進行運算時，需要試商導致增加運算時間的缺陷而提出的改進算法。

基本介紹

• 中文名：Stein算法
• 性質：算法
• 套用：計算兩個數最大公約數
• 針對：增加運算時間的缺陷

歐幾里德算法缺陷

歐幾里德算法是計算兩個數最大公約數的傳統算法，無論從理論還是從實際效率上都是很好的。但是卻有一個致命的缺陷，這個缺陷在素數比較小的時候一般是感覺不到的，只有在大素數時才會顯現出來。

一般實際套用中的整數很少會超過64位（當然現在已經允許128位了），對於這樣的整數，計算兩個數之間的模是很簡單的。對於字長為32位的平台，計算兩個不超過32位的整數的模，只需要一個指令周期，而計算64位以下的整數模，也不過幾個周期而已。但是對於更大的素數，這樣的計算過程就不得不由用戶來設計，為了計算兩個超過64位的整數的模，用戶也許不得不採用類似於多位數除法手算過程中的試商法，這個過程不但複雜，而且消耗了很多CPU時間。對於現代密碼算法，要求計算128位以上的素數的情況比比皆是，設計這樣的程式迫切希望能夠拋棄除法和取模。

算法思想

由J. Stein 1961年提出的Stein算法很好的解決了歐幾里德算法中的這個缺陷，Stein算法只有整數的移位和加減法，為了說明Stein算法的正確性，首先必須注意到以下結論：

gcd(a,a)=a，也就是一個數和其自身的公約數仍是其自身。

gcd(ka,kb)=k gcd(a,b)，也就是最大公約數運算和倍乘運算可以交換。特殊地，當k=2時，說明兩個偶數的最大公約數必然能被2整除。

當k與b互為質數，gcd(ka,b)=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

以前的算法有錯誤，因為cn根本就沒有用到。我編程的時候才發現。現在我已經修正了這個錯誤。

兩種算法的對比

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

C++/java 實現

#define　CHECK(a)　(!(1&(a)))//判斷是否被2整除#define　CLEAN2(a)　while(CHECK(a))a=a>>=1//移除非公因子的2#define　BIGERA　if(a<b)(t=a,a=b,b=t)//取較大的數為aint　gcd(int　a,int　b){　int　c_2=0,t;　while((CHECK(a))&&(CHECK(b))){　a=a>>=1;b=b>>=1;c_2++;　}　CLEAN2(a);　CLEAN2(b);　BIGERA;　while(a=((a-b)>>1)){　CLEAN2(a);　BIGERA;　}　return　b<<c_2;}//注意，由於中文百科會刪除空格，以上空格皆為全形空格//異或交換比中間量賦值慢//前一版在循環中有重複的判斷，在循環的結尾和絕對值處各判斷了一次ab的大小關係，所以效率受到影響//經過測試，比其他用遞歸的實現稍快，但是非常有限

Php實現

function steincore($a,$b) {    if ($a == 0) {return $b;}    if ($b == 0) {return $a;}    while (($a & 1) == 0) {        $a=$a>>1;    }    if ($a<$b) {        $b=($b-$a)>>1;        return steincore($b,$a);    }    else {        $a=($a-$b)>>1;        return steincore($a,$b);    } }function stein($a,$b) {    $c=0;    while ((($a & 1)==0)&&(( $b & 1 )==0)) {        $a=$a>>1;        $b=$b>>1;        $c=$c+1;    }    if (($a & 1) == 0) {        $a=$a>>1;        return steincore($a,$b)<<$c;    }    else {        return steincore($b,$a)<<$c;    }}//通過改善流程省去了低效的a，b交換，同時用while循環代替了過多的遞歸，以節約空間和時間。//疊代過程化為子過程以減少不必要的判斷//子過程中只有減法運算後的值可能為偶數所以只需要對$a判斷是否為偶數

javascript實現

function stein(a,b) {    if (a==0) {return b;}    if (b==0) {return a;}    var c=0;    while (((a & 1)==0)&&(( b & 1 )==0)) {        a=a>>1;        b=b>>1;        c=c+1;    }    while ((a & 1)==0) {        a=a>>1;    }    while ((b & 1)==0) {        b=b>>1;    }    if (a<b) {        b=(b-a)>>1;        var r=stein(b,a)<<c;        return r;    }    else {        a=(a-b)>>1;        var r=stein(a,b)<<c;        return r;    }}//基於php版修改而成

最佳化的C實現

int gcdcore(int a,int b) {    if (a==0) return b;    if (b==0) return a;    while ((a & 0x1)==0) {        a=a>>1;    }    if (a<b) {        b=(b-a)>>1;        return gcdcore(b,a);    }    else {        a=(a-b)>>1;        return gcdcore(a,b);    }}int gcd(int a,int b) {    int c=0;    while (((a & 0x1)==0)&&(( b & 0x1 )==0)) {        a=a>>1;        b=b>>1;        c++;    }    if ((a & 0x1) == 0) {        a=a>>1;        return gcdcore(a,b)<<c;    }    else {        return gcdcore(b,a)<<c;    }}

Ruby實現

def steincore(a,b)        a === 0 ? (return b) :        b === 0 ? (return a) :        a=a>>1 while ((a&1)===0)        a<b ? (return steincore((b-a)>>1,a)) : (return steincore((a-b)>>1,b))enddef GCD(a,b)        c=0        a,b,c=a>>1,b>>1,c+1 while ((a & 1 === 0)&&( b & 1 === 0))        ((a & 1)===0) ? (a=steincore(a,b)) : (a=steincore(b,a))        return (a << c)end

Pascal實現

function stein(a,b :longint):longint;

Begin

if a < b then

begin

a:=a xor b;

b:=a xor b;

a:=a xor b;

end;

if b = 0 then exit(a);

If (a and 1 = 0) and (b and 1 = 0) then exit(stein(a shr 1,b shr 1) shl 1);

If (a and 1 = 0) then exit(stein(a shr 1,b));

If (b and 1 = 0) then exit(stein(a,b shr 1));

exit(stein((a+b）shr 1,(a-b)shr 1));

End;