歐幾里德算法缺陷
歐幾里德算法是計算兩個數最大公約數的傳統算法，無論從理論還是從實際效率上都是很好的。但是卻有一個致命的缺陷，這個缺陷在素數比較小的時候一般是感覺不到的，只有在大素數時才會顯現出來。
一般實際套用中的整數很少會超過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 實現
Php實現
javascript實現
最佳化的C實現
Ruby實現
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;
