看完微軟大神寫的求平均值代碼，我意識到自己還是too young了

行早 2022-02-13 13:22:59 來源：量子位

博雯 發(fā)自 凹非寺

量子位 | 公眾號 QbitAI

取整求個無符號整數的平均值，居然也能整出花兒來？

這不，微軟大神Raymond Chen最近的一篇長文直接引爆外網技術平臺，引發(fā)無數討論：

無數人點進去時無比自信：不就是一個簡單的相加后除二的小學生編程題嗎？

unsigned average(unsigned a, unsigned b)
{
return (a + b) / 2;
}

但跟著大神的一路深挖，卻逐漸目瞪狗呆……

沒那么簡單的求平均值

先從開頭提到的小學生都會的方法看起，這個簡單的方法有個致命的缺陷：

如果無符號整數的長度為32位，那么如果兩個相加的值都為最大長度的一半，那么僅在第一步相加時，就會發(fā)生內存溢出。

也就是average（0x80000000U, 0x80000000U）=0。

不過解決方法也不少，大多數有經驗的開發(fā)者首先能想到的，就是預先限制相加的數字長度，避免溢出。

具體有兩種方法：

1、當知道相加的兩個無符號整數中的較大值時，減去較小值再除二，以提前減少長度：

unsigned average(unsigned low, unsigned high)
{
return low + (high - low) / 2;
}

2、對兩個無符號整數預先進行除法，同時通過按位與修正低位數字，保證在兩個整數都為奇數時，結果仍然正確。

（順帶一提，這是一個被申請了專利的方法，2016年過期）

unsigned average(unsigned a, unsigned b)
{
return (a / 2) + (b / 2) + (a & b & 1);
}

這兩個都是較為常見的思路，不少網友也表示，自己最快想到的就是2016年專利方法。

同樣能被廣大網友快速想到的方法還有SWAR（SIMD within a register）：

unsigned average(unsigned a, unsigned b)
{
return (a & b) + (a ^ b) / 2;// 變體 (a ^ b) + (a & b) * 2

以及C++ 20版本中的std: : midpoint函數。

接下來，作者提出了第二種思路：

如果無符號整數是32位而本機寄存器大小是64位，或者編譯器支持多字運算，就可以將相加值強制轉化為長整型數據。

unsigned average(unsigned a, unsigned b)
{
// Suppose "unsigned" is a 32-bit type and
// "unsigned long long" is a 64-bit type.
return ((unsigned long long)a + b) / 2;
}

不過，這里有一個需要特別注意的點：

必須要保證64位寄存器的前32位都為0，才不會影響剩余的32位值。

像是x86-64和aarch64這些架構會自動將32位值零擴展為64位值：

// x86-64: Assume ecx = a, edx = b, upper 32 bits unknown
mov eax, ecx ; rax = ecx zero-extended to 64-bit value
mov edx, edx ; rdx = edx zero-extended to 64-bit value
add rax, rdx ; 64-bit addition: rax = rax + rdx
shr rax, 1 ; 64-bit shift: rax = rax >> 1
; result is zero-extended
; Answer in eax

// AArch64 (ARM 64-bit): Assume w0 = a, w1 = b, upper 32 bits unknown
uxtw x0, w0 ; x0 = w0 zero-extended to 64-bit value
uxtw x1, w1 ; x1 = w1 zero-extended to 64-bit value
add x0, x1 ; 64-bit addition: x0 = x0 + x1
ubfx x0, x0, 1, 32 ; Extract bits 1 through 32 from result
; (shift + zero-extend in one instruction)
; Answer in x0

而Alpha AXP、mips64等架構則會將32位值符號擴展為64位值。

這種時候，就需要額外增加歸零的指令，比如通過向左進位兩字的刪除指令rldicl：

// Alpha AXP: Assume a0 = a, a1 = b, both in canonical form
insll a0, #0, a0 ; a0 = a0 zero-extended to 64-bit value
insll a1, #0, a1 ; a1 = a1 zero-extended to 64-bit value
addq a0, a1, v0 ; 64-bit addition: v0 = a0 + a1
srl v0, #1, v0 ; 64-bit shift: v0 = v0 >> 1
addl zero, v0, v0 ; Force canonical form
; Answer in v0

// MIPS64: Assume a0 = a, a1 = b, sign-extended
dext a0, a0, 0, 32 ; Zero-extend a0 to 64-bit value
dext a1, a1, 0, 32 ; Zero-extend a1 to 64-bit value
daddu v0, a0, a1 ; 64-bit addition: v0 = a0 + a1
dsrl v0, v0, #1 ; 64-bit shift: v0 = v0 >> 1
sll v0, #0, v0 ; Sign-extend result
; Answer in v0

// Power64: Assume r3 = a, r4 = b, zero-extended
add r3, r3, r4 ; 64-bit addition: r3 = r3 + r4
rldicl r3, r3, 63, 32 ; Extract bits 63 through 32 from result
; (shift + zero-extend in one instruction)
; result in r3

或者直接訪問比本機寄存器更大的SIMD寄存器，當然，從通用寄存器跨越到SIMD寄存器肯定也會增加內存消耗。

如果電腦的處理器支持進位加法，那么還可以采用第三種思路。

這時，如果寄存器大小為n位，那么兩個n位的無符號整數的和就可以理解為n+1位，通過RCR（帶進位循環(huán)右移）指令，就可以得到正確的平均值，且不損失溢出的位。

帶進位循環(huán)右移

// x86-32
mov eax, a
add eax, b ; Add, overflow goes into carry bit
rcr eax, 1 ; Rotate right one place through carry

// x86-64
mov rax, a
add rax, b ; Add, overflow goes into carry bit
rcr rax, 1 ; Rotate right one place through carry

// 32-bit ARM (A32)
mov r0, a
adds r0, b ; Add, overflow goes into carry bit
rrx r0 ; Rotate right one place through carry

// SH-3
clrt ; Clear T flag
mov a, r0
addc b, r0 ; r0 = r0 + b + T, overflow goes into T bit
rotcr r0 ; Rotate right one place through carry

那如果處理器不支持帶進位循環(huán)右移操作呢？

也可以使用內循環(huán)（rotation intrinsic）：

unsigned average(unsigned a, unsigned b)
{
#if defined(_MSC_VER)
unsigned sum;
auto carry = _addcarry_u32(0, a, b, &sum);
sum = (sum & ~1) | carry;
return _rotr(sum, 1);
#elif defined(__clang__)
unsigned carry;
sum = (sum & ~1) | carry;
auto sum = __builtin_addc(a, b, 0, &carry);
return __builtin_rotateright32(sum, 1);
#else
#error Unsupported compiler.
#endif
}

結果是，x86架構下的代碼生成沒有發(fā)生什么變化，MSCver架構下的代碼生成變得更糟，而arm-thumb2的clang 的代碼生成更好了。

// _MSC_VER
mov ecx, a
add ecx, b ; Add, overflow goes into carry bit
setc al ; al = 1 if carry set
and ecx, -2 ; Clear bottom bit
movzx ecx, al ; Zero-extend byte to 32-bit value
or eax, ecx ; Combine
ror ear, 1 ; Rotate right one position
; Result in eax

// __clang__
mov ecx, a
add ecx, b ; Add, overflow goes into carry bit
setc al ; al = 1 if carry set
shld eax, ecx, 31 ; Shift left 64-bit value

// __clang__ with ARM-Thumb2
movs r2, #0 ; Prepare to receive carry
adds r0, r0, r1 ; Calculate sum with flags
adcs r2, r2 ; r2 holds carry
lsrs r0, r0, #1 ; Shift sum right one position
lsls r1, r2, #31 ; Move carry to bit 31
adds r0, r1, r0 ; Combine

微軟大神的思考們

Raymond Chen1992年加入微軟，迄今為止已任職25年，做UEX-Shell，也參與Windows開發(fā)，Windows系統(tǒng)的很多最初UI架構就是他搞起來的。

他在MSDN 上建立的blogThe Old New Thing也是業(yè)內非常出名的純技術向產出網站。

這篇博客的評論區(qū)們也是微軟的各路大神出沒，繼續(xù)深入探討。

有人提出了新方法，在MIPS ASM共有36個循環(huán)：

unsigned avg(unsigned a, unsigned b
{
return (a & b) + (a ^ b) / 2;
}

// lw $3,8($fp) # 5
// lw $2,12($fp) # 5
// and $3,$3,$2 # 4
// lw $4,8($fp) # 5
// lw $2,12($fp) # 5
// xor $2,$4,$2 # 4
// srl $2,$2,1 # 4
// addu $2,$3,$2 # 4

有人針對2016年專利法表示，與其用(a / 2) + (b / 2) + (a & b & 1)的方法，為啥不直接把 (a & 1) & ( b & 1 ) ) 作為進位放入加法器中計算呢？

還有人在評論區(qū)推薦了TopSpeed編譯器，能夠通過指定合適的代碼字節(jié)和調用約定來定義一個內聯(lián)函數，以解決“乘除結果是16位，中間計算值卻不是”的情況。

只能說，學無止境啊。

原文：
https://devblogs.microsoft.com/oldnewthing/20220207-00/?p=106223

參考鏈接：
https://news.ycombinator.com/item?id=30252263