程序人生 - 微軟大神「玩」出新花樣，求平均值代碼還能這樣寫？－鑽石舞台

編譯 | 馬超責編 | 蘇宓

出品 | CSDN（ID：CSDNnews）

近日，微軟神級人物Raymond Chen在個人博客上，發布了一篇關於《如何計算平均值》的博文。這個話題雖然看似平淡無奇，卻意外引爆技術圈，並帶來無數討論。

看完這篇博客之後，也讓我感嘆於國外技術討論氛圍的濃烈，雖然這一話題切入點非常簡單，但是最終能夠升華至編程之道層面的舉輕若重的文章，接下來，我們不妨一起來看看。

有關求平均數算法的最初版本

有關如何求平均數這個問題，Raymond Chen並沒有從一開始就炫技，而是循序漸進先放了一段最普通的實現，如下：

unsignedaverage(unsigneda,unsignedb){return(a+b)/2;}

相信絕大多數程序員都能一眼看出這種方法中可能隱藏的錯誤，那就是無法處理值溢出的問題，在Raymond的原文當中用「if unsigned integers are 32 bits wide, then it says that average(0x80000000U, 0x80000000U) is zero.」一句話來總結。

也就是說一旦（a+b)已經溢出，也就是大於unsign類型所能表示的最大整改，那麼其計算結果將是average（0x80000000U, 0x80000000U）=0。

不過筆者在這裡需要指出0x80000000U是x86平台特有的一個溢出表示方法，即indefinite integer value（不確定數值），不過同樣是溢出ARM等RISC架構處理則非常清晰和簡單，在上溢出或下溢出時，保留整型能表示的最大值或最小值，對照比較如下：

CPU

溢出值轉為long

變量保留值說明

x86

範圍0x8000000000000000

indefinite integer value

x86

範圍0x8000000000000000

indefinite integer value

ARM

範圍0x7FFFFFFFFFFFFFFF

變量賦值最大的正數

ARM

範圍0x8000000000000000

變量賦值最大的正數

因此這段代碼在ARM平台上運行時，如果出現溢出情況也並不會返回0，而會是該類型表示最大整數的一半，當然這個最大整數根據處理器的字長不同可能會有所變化。

return(a+b)/2

低調的改進版本

接下來Raymond又給出了幾種考慮溢出處理，同時又兼顧空間複雜度的方案：

1、變形法：

也就是將（a+b)/2變形，首先找到a和b當中較大的值，設為high，較小的值設為low，然後把（a+b)/2變成high-(high-low)/2或者low+(high-low)/2，如下：

unsignedaverage(unsignedlow,unsignedhigh){returnlow+(high-low)/2;}

這種方法所需要的運算量是先進行一次比較以確定兩個輸入的大小，然後還需要再做兩次加法（在計算機運算中加法和減法其實是基本等效的）和一次除法，最終得到答案。

2、除法前置方案：

也就是先對兩個輸入進行除2操作，即把（a+b)/2轉換為a/2+b/2，當然這種方法需要考慮個位丟失的問題，比如說1/2在整形運算當中的結果會是0，因此1/2+1/2的結果是0而不是1，此時需要把兩個輸入的個位提取出來進行修正，具體如下：

unsignedaverage(unsigneda,unsignedb){return(a/2)+(b/2)+(a&b&1);}

這個算法當中的計算量是兩次除法，兩次加法和一次與運算操作。

3.SWAR法

SWAR法也非常的巧妙，它的本質思路就是把求平均值變成位運算，位操作其實就是二進制的操作，如果我們按位考慮輸入值與輸出結果的對應關係，那麼會有以下的需求要點：

輸入都是0，輸出結果是0

輸入都是1，輸出是1

輸入是一個0一個1，那麼輸出結果就是1/2

而滿足以上條件的位運算，是與運算加上異常運算除2的結果，即(a and b) + (a xor b )/2，如下：

unsignedaverage(unsigneda,unsignedb){return(a&b)+(a^b)/2;//變體(a^b)+(a&b)*2}

至於(a and b) + (a xor b )/2這個等式為什麼能滿足求平均值的要求，大家根據各種輸入的情況都列一下就一目了然了。在這種方案下的計算量是兩次位運算、一次加法運算以及一次除法運算來完成。

空間換時間的改進版本

在算法設計當中有一個最基本的常識，空間複雜度與時間複雜度是對蹺蹺板，上一節的儲多算法當中，基本都是犧牲時間複雜度為代價來換取對於溢出的正確處理，那麼反過來講也完全可以用空間換時間，比如現在我們大多數的終端電腦都是64位機了，沒必要為了32位長的整形溢出問題而煩惱，直接把類型轉換為Long再計算結果就可以了。

unsignedaverage(unsigneda,unsignedb){//Suppose"unsigned"isa32-bittypeand//"unsignedlonglong"isa64-bittype.return((unsignedlonglong)a+b)/2;}

但是只要涉及的轉換就又要針對不同架構的處理器進行特殊處理了，比如x86的64位處理器在進行32位整形轉換為64位長整形時會自動將高32位的值填為0：

//x86-64:Assumeecx=a,edx=b,upper32bitsunknownmoveax,ecx;rax=ecxzero-extendedto64-bitvaluemovedx,edx;rdx=edxzero-extendedto64-bitvalueaddrax,rdx;64-bitaddition:rax=rax+rdxshrrax,1;64-bitshift:rax=rax>>1;resultiszero-extended;Answerineax//AArch64(ARM64-bit):Assumew0=a,w1=b,upper32bitsunknownuxtwx0,w0;x0=w0zero-extendedto64-bitvalueuxtwx1,w1;x1=w1zero-extendedto64-bitvalueaddx0,x1;64-bitaddition:x0=x0+x1ubfxx0,x0,1,32;Extractbits1through32fromresult;(shift+zero-extendinoneinstruction);Answerinx0

Mips64等架構則會將32位的整形轉換為有符號擴展的類型。這時候就需要增加rldicl等刪除符號的指令做特殊處理。

//AlphaAXP:Assumea0=a,a1=b,bothincanonicalforminslla0,#0,a0;a0=a0zero-extendedto64-bitvalueinslla1,#0,a1;a1=a1zero-extendedto64-bitvalueaddqa0,a1,v0;64-bitaddition:v0=a0+a1srlv0,#1,v0;64-bitshift:v0=v0>>1addlzero,v0,v0;Forcecanonicalform;Answerinv0//MIPS64:Assumea0=a,a1=b,sign-extendeddexta0,a0,0,32;Zero-extenda0to64-bitvaluedexta1,a1,0,32;Zero-extenda1to64-bitvaluedadduv0,a0,a1;64-bitaddition:v0=a0+a1dsrlv0,v0,#1;64-bitshift:v0=v0>>1sllv0,#0,v0;Sign-extendresult;Answerinv0//Power64:Assumer3=a,r4=b,zero-extendedaddr3,r3,r4;64-bitaddition:r3=r3+r4rldiclr3,r3,63,32;Extractbits63through32fromresult;(shift+zero-extendinoneinstruction);resultinr3

不過這種向更高位類型轉換的方案也有一定問題，那就是空間的浪費，因為我原本只需要1位去處理溢出就好了，但是做了轉換之後我卻用了白白消費了31位的空間沒有利用。

利用進位處理溢出的改進版本

在現代CPU當中大多都帶有Carry bit(這裡指進位位，不是C位的意思）功能。通過讀取Carry bit的信息，就能達到在不浪費空間的情況下處理溢出的問題。比如在X86-32位處理器的代碼如下：

//x86-32moveax,aaddeax,b;Add,overflowgoesintocarrybitrcreax,1;Rotaterightoneplacethroughcarry//x86-64movrax,aaddrax,b;Add,overflowgoesintocarrybitrcrrax,1;Rotaterightoneplacethroughcarry//32-bitARM(A32)movr0,aaddsr0,b;Add,overflowgoesintocarrybitrrxr0;Rotaterightoneplacethroughcarry//SH-3clrt;ClearTflagmova,r0addcb,r0;r0=r0+b+T,overflowgoesintoTbitrotcrr0;Rotaterightoneplacethroughcarry

而對於那些沒有Carry bit功能的處理器來說，也可以通過自定義carry bit變量的方式來解決這個問題。如下：

unsignedaverage(unsigneda,unsignedb){#ifdefined(_MSC_VER)unsignedsum;autocarry=_addcarry_u32(0,a,b,&sum);return_rotr1_carry(sum,carry);//missingintrinsic!#elifdefined(__clang__)unsignedcarry;autosum=_builtin_adc(a,b,0,&carry);return_builtin_rotateright1throughcarry(sum,carry);//missingintrinsic!#elifdefined(__GNUC__)unsignedsum;autocarry=__builtin_add_overflow(a,b,&sum);return_builtin_rotateright1throughcarry(sum,carry);//missingintrinsic!#else#errorUnsupportedcompiler.#endif}

對應arm-thumb2的clang 匯編代碼如下：

//__clang__withARM-Thumb2movsr2,#0;Preparetoreceivecarryaddsr0,r0,r1;Calculatesumwithflagsadcsr2,r2;r2holdscarrylsrsr0,r0,#1;Shiftsumrightonepositionlslsr1,r2,#31;Movecarrytobit31addsr0,r1,r0;Combine

Quake3中「神」一樣的代碼

可以看到Raymond的博客先從一個簡單問題入手，逐步提出問題並給出解決方案，是一篇闡述編程之道的上乘之作，接下來請允許筆者再推薦一下《Quake3》當中的神級代碼。

《Quake3》這款3D遊戲當年可以在幾十兆內存的環境下跑得飛起，和目前動輒要求幾十G顯存的所謂3A大作形成鮮明對比，而《Quake3》取得這種性價比奇蹟的關鍵在於把代碼寫得像神創造的一樣。

《Quake3》最大的貢獻莫過於提出使用平方根倒數速算法，並引入了0x5f3759df這樣一個魔法數，目前這段代碼的開源地址在：https://github.com/raspberrypi/quake3/blob/8d89a2a3c1707bf0f75b2ea26645b872e97c0b95/code/qcommon/q_math.c

如下：

floatQ_rsqrt(floatnumber){floatint_tt;floatx2,y;constfloatthreehalfs=1.5F;x2=number*0.5F;t.f=number;t.i=0x5f3759df-(t.i>>1);//whatthefuck?y=t.f;y=y*(threehalfs-(x2*y*y));//1stiteration//y=y*(threehalfs-(x2*y*y));//2nditeration,thiscanberemovedreturny;}

這個算法的輸入是一個float類型的浮點數，首先將輸入右移一次（除以2），並用十六進制「魔術數字」0x5f3759df減去右移之後的數字，這樣即可得對輸入的浮點數的平方根倒數的首次近似值；而後重新將其作為原來的浮點數，以牛頓迭代法迭代，目前來看迭代一次即可滿足要求，這個算法避免了大量的浮點計算，比直接使用浮點數除法要快四倍，大幅提升了平方根倒數運算的效率。