關於 C/C++ 負整數的 shift 以及 arithmetic right shift 的小技巧

在看 《Tonc: Whirlwind Tour of ARM Assembly》 的時候學到一個小技巧: 若 n 是有號 32bit 整數, 則 n>>31 相當於 n >= 0 ? 0 : -1。

看完 arithmetic shift 的定義後, 才發現我沒弄懂 arithmetic right shift, 原來是一直補進 signed bit (most significant bit), 當它是負數時, 等於一直補 1 進來。所以, -1 不管右移幾次, 結果仍是 -1, 不是 0。

將這訊息發到 G+ 和 plurk 後, 看到許多人提到 C、C++ 沒定義 signed shift 的行為, 所以用的時候要小心。jserv 提到 compiler 可以加參數改變位元運算的結果。

查看 gcc 的說明, 得知 gcc 有實作 arithmetic right shift, 但看不懂 left shift 的部份。Scott 提到這篇有解釋, 然後我明白了一件事: 關於 C99、C++0x draft 的定義, 我都看不太懂啊 ...。孩子的教育不能等, 要好好學英文。

btw, Louis 提到一個有趣的應用, ((n>>31)^n) - (n>>31) 等同於 abs(n)。拆解這個式子可得知, n 是正數時:

• n>>31 是 0
• n^0 仍是 n
• n - 0 仍是 n

n 是負數時:

• n>>31 是 -1
• -1 以 2 補數表示為 0xFFFFFFFF
• n^(-1) 等同於 1 補數運算
• 最後再減 -1, 得到 2 補數運算的值

順便試著寫了以下的 C code, 再編成 assembly code, 想說看看會差多少:

int abs1(int n) {
    return n >= 0 ? n : -n;
}

int abs2(int n) {
    if (n >= 0)
        return n;
    return -n;
}

int abs3(int n) {
    return ((n>>31)^n) - (n>>31);
}

結果 abs1 和 abs3 差不多, 以下是 abs1 用 gcc 4.4.3 -O3 編後的結果:

abs1:
 pushl %ebp
 movl %esp, %ebp
 movl 8(%ebp), %eax
 popl %ebp
 movl %eax, %edx
 sarl $31, %edx
 xorl %edx, %eax
 subl %edx, %eax
 ret

gcc 偷吃步用 arithmetic right shift 避免使用 branch。不知這是針對特例的最佳化, 還是用某種神祕的機制算出來的。附帶一提, 我用不同版本的 gcc 編出來的結果不同, 用 4.5.2 最佳化編譯後, abs1 和 abs3 的結果完全一樣, 只差在用的 register 不同。

另外, 若用 ARM assembly 的話, abs() 可以寫成:

@@ 參數由 r0 傳入, 結果寫到 r0
cmp     r0, #0
sublt     r0, 0, r0

而用 arithmetic right shift 的版本則是:

eor     r1, r0, r0, asr #31
sub     r0, r1, r0, asr #31

都是兩個算數指令, 速度應該一樣吧。我沒裝 ARM 開發環境, 應該沒寫錯吧。藉機練習寫看看, 體會指令附加 branch 以及 barrel shifter 的便利之處。指令詳細說明, 見這篇的 "23.3.2. Data instructions"。"r0, asr #31" 那部份是 Op2, 可在同一 cycle 內取得 shift 後的值。

用 XOR 檢查加法後是否有溢位

看 Python intobject.c 時看到的:

// intobject.c
static PyObject *
int_add(PyIntObject *v, PyIntObject *w)
{
    register long a, b, x;
    CONVERT_TO_LONG(v, a);
    CONVERT_TO_LONG(w, b);
    x = a + b;
    if ((x^a) >= 0 || (x^b) >= 0)
        return PyInt_FromLong(x);
    return PyLong_Type.tp_as_number->nb_add((PyObject *)v, (PyObject *)w);
}

上面這段的

if ((x^a) >= 0 || (x^b) >= 0)

是在檢查是否有溢位, 看一陣子才想通, 概念如下:

• 兩個正數相加變成負數表示溢位
• 兩個負數相加變成正數也表示溢位
• 除了這兩者外沒有其它情況
• 綜合上述的條件, 可以簡化成: 兩數相加後的正負號和原本兩者都不同 <---> 溢位
• 目前 CPU 都採用二補數系統, 最高 bit 可看作是 sign bit
• XOR 剛好可以滿足簡化的條件, 其它 bit 無關緊要, 運算完看 sign bit 即可

位元運算的世界真神奇啊, 可以省下不少工。

計算 64 bit long 內有幾個 bit 為 1 的速算法

原本在讀 bit array 為啥威, 結果看到 Hamming weight, 即計算一個 "string" 內有幾個 1 的速算法。

基本版 (popcount_1) 的想法還算好懂, 典型的 divide & conquer, 對半切的解法, 看了開頭就猜到後面要怎麼做了。先以 2 bit 為單位算出有幾個 1, 並存在原本的那 2 bit 裡, 再用同樣的手法以 4 bit 為單位算出 4 bit 裡有幾個 1, 並存在原本的那 4 bit 裡, ..., 最後就會在較低的 32 bit 裡得到全部答案。

進階版 (popcount_2) 的第一個寫法有點妙, 我最後才參透它。8、16、32 bit 時最簡單, 若 bit 數的最大值不會超出儲存的空間, 就不用擔心相加會溢位, 而能直接位移並做相加。2 bit 的作法沒變, 4 bit 的作法和 8 bit 開始的作法一樣, 只是先清乾淨「高位元」的地方, 確保後面處理 8、16、32 bit 時不會溢位。

回頭看 x -= (x >> 1) & m1, 將式子展開會得到 x = x - ((x >> 1) & m1)。而原本的算法是 (x & m1) + ((x >> 1) & m1), 將兩式合起來看會得到:

x - ((x >> 1) & m1) = (x & m1) + ((x >> 1) & m1)
x - (x & m1) = ((x >> 1) & m1) + ((x >> 1) & m1)

整理左式:

x - (x & m1) = x & (m1 << 1)

整理右式

((x >> 1) & m1) + ((x >> 1) & m1)
= 2 * ((x >> 1) & m1) 
= ((x >> 1) & m1) << 1
= x & (m1 << 1)

所以這個算法沒錯, 將上面的過程反過來推導, 就會發覺可以用 x - ((x >> 1) & m1) 取代原本的式子以省下一個運算。

H = (x >> 1) & m1
L = x & m1
H + L 
= H + H - H + L
= 2H - H + L 
= (x - L) - H + L
= x - H
= x - (x >> 1) & m1

位元運算的世界真是太奇妙了, 後面的更進階版有機會再來參透吧。