2022q1 Homework3 (quiz3)

contributed by < tinyynoob >

作業要求

測驗 3

解釋程式碼








#include <stdint.h>
uint8_t rev8(uint8_t x)
{
    x = (x >> 4) | (x << 4);               
    x = ((x & 0xCC) >> 2) | (EXP2);
    x = ((x & 0xAA) >> 1) | (EXP3);
    return x;
}

為了對 8 bit 數字進行翻轉，這裡的作法以類似 top-down 的方式往下處理，每次交換左右半部，直到每個半部的 size 都為 1 即完成所有翻轉，圖解如下：

因為 0xcc = 11001100 ，所以第 5 行需要補上的就是另一半的操作，因此 EXP2 為 (x & 0x33) << 2 。

同理 EXP1 為 (x & 0x55) << 1 。

測驗 5

解釋程式碼
































#include <limits.h>
int divide(int dividend, int divisor)
{
    int signal = 1;
    unsigned int dvd = dividend;
    if (dividend < 0) {
        signal *= -1;
        dvd = ~dvd + 1;
    }

    unsigned int dvs = divisor;
    if (divisor < 0) {
        signal *= -1;
        dvs = ~dvs + 1;
    }

    int shift = 0;
    while (dvd > (EXP6))
        shift++;

    unsigned int res = 0;
    while (dvd >= dvs) {
        while (dvd < (dvs << shift))
            shift--;                         
        res |= (unsigned int) 1 << shift;
        EXP7;
    }

    if (signal == 1 && res >= INT_MAX)
        return INT_MAX;
    return res * signal;
}

可以看出第 4 至 15 行在把 dvd 跟 dvs 都轉成非負的值，並把 sign 值存放到 signal 。另外也可以看出缺乏對 divisor 為 0 的檢查。

這裡由於是 int 轉成 unsigned int ，所以取二補數並不會因 INT_MIN 導致 undefined behaviour ，以 4 bit 舉例說明：

1000    # INT_MIN of 4-bit signed integer, which is -8
1000    # two's complement of 1000, it is unsigned 8

因此第 8 及 14 行會得出正確的結果。

接著後面就是要想辦法用 shift 和 subtract 來做出除法。因為除法是從最高位開始，所以這裡的第 17~19 行即是要推出最高位並以 shift 紀錄，大概會想到兩個可能的答案

dvs << shift
dvs << shift + 1

有例子會比較容易推測，所以我們先假定 dvd = 16 和 dvs = 3 ，則相除結果應為 0101 。除了第 26 行未知以外，唯一在更新答案的就是第 25 行，其功能為將特定位元設為 1 ，因此可以推斷第一次執行到第 25 行時 shift 應該要為 2 。而因為第 23 ~ 24 行結束後會保證 dvd >= (dvs << shift) ，所以第 17~19 行把 shift 起始值稍微設超過也沒關係，因此 EXP6 可以選成較為簡潔的 dvs << shift 。

       1 
   --------
11 / 10000
      1100    # (11) << 2
    ------

在繼續除之前，我們需要更新 dvd ，所以 EXP7 就是 dvd -= (dvs << shift) ，得：

       1 
   --------
11 / 10000
      1100
    ------
       100

接著繼續想，會發現 23~24 行幫助我們找到哪一位可以繼續除，全部做完就會像這樣：

不過因為我們是用 unsigned int 在除，而商的回傳型態是 int ，所以前面所提及的 INT_MIN 還是可能造成問題。一樣以 4-bit 為例， 1000 除以 1 得 1000 ，但因為 1000 原本代表的是負數所以 signal 會被標記為 -1 ，那麼到第 31 行就會發生 1000 * (-1) 的狀況

abs(INT_MIN) 是未定義行為，那 (INT_MIN) * (-1) 呢？
測試似乎會得到原值

補齊後程式碼































int divide(int dividend, int divisor)
{
    int signal = 1;
    unsigned int dvd = dividend;
    if (dividend < 0) {
        signal *= -1;
        dvd = ~dvd + 1;
    }

    unsigned int dvs = divisor;
    if (divisor < 0) {
        signal *= -1;
        dvs = ~dvs + 1;
    }

    int shift = 0;
    while (dvd > (dvs << shift))
        shift++;

    unsigned int res = 0;
    while (dvd >= dvs) {
        while (dvd < (dvs << shift))
            shift--;
        res |= (unsigned int) 1 << shift;
        dvd -= (dvs << shift);
    }

    if (signal == 1 && res >= INT_MAX)
        return INT_MAX;
    return res * signal;
}

不理解 28, 29 行的功能

改進並實作

錯誤處理

我發現它沒有對 divisor 為 0 的檢查，結果會使程式陷入無窮迴圈。

使用一些 corner cases 測試可以發現：

2147483647 / 1 = 2147483647
2147483647 / -1 = -2147483647
-2147483648 / 1 = -2147483648
-2147483648 / -1 = 2147483647

在 INT_MIN / -1 的情境中，也因為結果出現 32 位有號數無法表示的

2^{32}

導致答案不如預期。

為了了解在慣例上這兩種情況該如何處置，我去測試了 standard library 中的 div() 。
呼叫 div(7, 0) 以及 div(INT_MIN, -1) 皆得到相同結果：

浮點數例外 (核心已傾印)

看來它們都是直接用 exception 的方式處理，查詢 man signal 找到這段敘述：

Integer division by zero has undefined result. On some architectures it will generate a SIGFPE signal. (Also dividing the most negative integer by -1 may generate SIGFPE.) Ignoring this signal might lead to an endless loop.

我想仿造使用相同的方式，因此在 divide() 函式的最前方加入：

    if (divisor == 0)
        raise(SIGFPE);
    else if (dividend == INT_MIN && divisor == -1)
        raise(SIGFPE);

測試後 signal 可被正確地觸發。

效能改善

我認為這裡 while 迴圈判斷 shift 的方式，可以使用 clz() 來取而代之，程式碼更動如下：

+   int shift = __builtin_clz(dvs) - __builtin_clz(dvd);
-   int shift = 0;
-   while (dvd > (dvs << shift))
-       shift++;

    unsigned int res = 0;

    while (dvd >= dvs) {
+       if (dvd < (dvs << shift))
+           shift--;
-       while (dvd < (dvs << shift))
-           shift--;
        res |= (unsigned int) 1 << shift;
        dvd -= (dvs << shift);
    }

嘗試測量 10000 組隨機數進行除法所耗費的時間，得：

number of cases: 10000
origial consumption: 966675
v1 consumption: 148825867

結果顯示改進版程式碼的耗費時間竟為原版的

10^{2}

倍，與我的預期截然不同。

working on…

測驗 7

解釋程式碼

給定 32 位的 unsigned v ，這題想計算出需要幾個 bit 才能表示 v ，顯然 ret

\in [0, 32]

，且 ret 是由 v 的 most significant set bit 所決定。

根據二進制的原理，我們知道任意 v 必定能表示為

\sum_{i = 0}^{31} b_{i} \cdot 2^{i} where b_{i} \in {0, 1}

i

值最大且為 1 的

b_{i}

就是我們要的 ret ，這裡想透過類似二分搜尋的方式來進行查找。

首先把 32 個 bit 分為兩半，數字代表 most significant set bit 在由下往上的第幾位。

v 至少要是

binary 1 \underset{16 個 0}{\underset{⏟}{00 \dots 00}}

才能落在左邊，因此當 v

\geq 2^{16}

時 v 屬於左半，否則 v 屬於右半。

在屬於左半的情況若是我們先把數字 16 存到某個變數，那麼左右兩個 branch 又變成一樣的問題，那就是求「最大的 1 在 16 位中的何處？」，於是繼續：

完成後一樣把結果加到 m 並繼續 4 位的搜尋。


















int ilog32(uint32_t v)
{
    int ret = v > 0;
    int m = (v > 0xFFFFU) << 4;
    v >>= m;
    ret |= m;
    m = (v > 0xFFU) << 3;
    v >>= m;
    ret |= m;
    m = (v > 0xFU) << 2;
    v >>= m;
    ret |= m;
    m = EXP10;
    v >>= m;
    ret |= m;
    EXP11;
    return ret;
}

在表示上，我們將前面的 v

\geq 2^{16}

更改為等價的 v

> 2^{16} - 1 = 0 xFFFFU

以利程式看起來的對稱性和二分關係。

在每一輪中我們藉由比較和左移把結果暫存到 m ，並用 | 來立起 ret 中對應的 bit ，接著右移 v 來讓兩種結果回歸同一個新問題。

繼續二分的結果， EXP10 顯然為 (v > 0x3U) << 1 。

最後第 16 行的 EXP11 就比較特別了，在程式走到這裡的時候僅剩下 2 個 bit ，所以要判斷是否

> 1

，我們可以列出 v 為 00, 01, 10, 11 的情況來思考，結果應該要分別對應到 0, 1, 2, 2。
然而我們目前不能產生兩位數 (已於 13 至 15 行佔用第 2 位)，此外仔細觀察可以發現第 1 位也早在第 3 行就已經被佔用，因此 0 和 1 的 case 其實已經判斷完畢。
那麼 EXP11 就應該以加法更新，並只在最後是 10 或 11 時更新，即為 ret += (v > 1) 。

探索 kernel 中的應用

研讀《Using de Bruijn Sequences to Index a 1 in a Computer Word》

雖然目前 architecture 上大都有提供 fls 或相似的 instruction，但為了 cross-platform 和 cross-compiler 的情境，我們有時仍需要 software 的 fls 處理。上方的 ilog32() 程式碼仍可能有 branch 產生，因此我們想找到更好的 branch-free 版本。

考慮

n

-bit 無號整數，此篇文章以

n = 8

、求 ffs 為例。

首先，我們需要一組長度為

n

的 de Bruijn sequence，這類特殊 binary sequence 所有

\log n

長度的 cyclic subsequence 完全不重複，也就是每個

\log n

長度的整數都會剛好出現一次。

例如這裡給定為 debru = 00011101，它所有長

\log n

的 cyclic subsequence 為：000, 001, 011, 111, 110, 101, 010, 100 (由左而右)，我們可以建立 look-up table 以快速查詢位置：

subsequence	index
`000`	0
`001`	1
`010`	6
`011`	2
`100`	7
`101`	5
`110`	4
`111`	3

接著要如何實現 ffs 呢？
我們先針對輸入值 weight 為 1 的 case，做法如下：

輸入值 x 乘以 debru
>> 5
& 7
查表

原理

以 x = 00010000 為例，首先因為 x 的 weight 為 1，將其乘以 debru 等同於對 debru 左移 ffs(x)

        00010000
      * 00011101
-----------------
0000000111010000
        @@@

步驟 2 和 3 固定取出 @@@ 位置的 bits

可以發現這些操作就等同取出 debru 的第 ffs(x) 個 cyclic subsequence，於是反查就可以輕易找到 ffs(x)。

以上僅適用於 weight 為 1 的 x，那要如何擴充到任意 x 呢？只要先對 x 使用 isolate first set 操作就行了，一個經典的方法是 x = x & (-x)。

需要注意到我們這裡的 ffs 並未考量 x = 0

解釋至此， fls() 也可以用相似的技巧配合 isolate last set 進行快速求解。

我們取 stackoverflow 討論串裡的 debru，並嘗試自己撰寫出 fls() 。

以 32-bit 為目標，首先撰寫輔助程式建出 look-up table

static const int debruijn_fls32_table[32] = {
    32, 31, 22, 30, 21, 18, 10, 29, 2,  20, 17, 15, 13, 9, 6,  28,
    1,  23, 19, 11, 3,  16, 14, 7,  24, 12, 4,  8,  25, 5, 26, 27};

於是可以寫出：

/* most significant = 1
 * least significant = 32
 * 0 = 0
 */
int fls32(uint32_t x)
{
    return (0 - !!x) & debruijn_fls32_table[(ils(x) * 0x07C4ACDDu) >> 27 & 31];
}

ils() function

// isolate last set
uint32_t ils(uint32_t n)
{
    n |= n >> 1;
    n |= n >> 2;
    n |= n >> 4;
    n |= n >> 8;
    n |= n >> 16;
    return n ^ (n >> 1);
}

依據之前 branch 版本的慣例，我將 MSB 視為 1、LSB 為 32，並補上 fls(0) = 0。

簡單用亂數進行測試，將 fls32() 的結果與 __builtin_clz() + 1 比對，完全沒出現錯誤回報，唯兩者對 0 會產生不同結果。

test code

#define TESTNUM 10000000

int main()
{
    srand(time(NULL));
    puts("start test:");
    for (int i = 0; i < TESTNUM; i++) {
        uint32_t x = rand() << 16 | rand();
        if (__builtin_clz(x) + 1 != fls32(x))
            printf("different answer when x = %d\n", x);
    }
    for (int i = 0; i < TESTNUM; i++) {
        uint32_t x = rand();
        if (__builtin_clz(x) + 1 != fls32(x))
            printf("different answer when x = %d\n", x);
    }
    puts("end test.");
    return 0;
}

編譯時期的 `ilog2()` 實作

測驗 10

解釋程式碼

DIVIDE_ROUND_CLOSEST(x, divisor)

此題想在整數除法中選擇較接近的整數作為答案，而非原本的無條件捨去。

先考慮正數的例子，從數學上來思考，假設除法結果落在 @ 會得到 1，而落在 # 會得 2，從數線看起來如：

0      1      2      3
|------|------|------|
       @@@@@@@#######    traditional computer result
    @@@@@@@#######       our goal

假設除出來位在

0      1      2      3
|------|------|------|
     x

如果我們能夠把它右移 0.5 ，那麼 x 就會落入傳統 @ 的範圍並被電腦映射到 1。

然而在整數中根本沒有 0.5 的概念，而追究 0.5 的來源即是整數除法中無法整除的部份，因此我們應想辦法去調整被除數。

假設正整數

a, b

在電腦上做除法得

a / b = c

，我們希望在算出

c

之後可以將它右移

0.5

，那麼可以將原算式改造成

(a + b * 0.5) / b

達成目的，

* 0.5

在電腦上可以使用 >> 1 。

接著我們需要將負整數除法也納入討論，在本機使用 gcc 編譯測試顯示 / 運算的結果會對稱於原點 0 ：

-6 / 3 = -2
-5 / 3 = -1
-4 / 3 = -1
-3 / 3 = -1
-2 / 3 = 0
-1 / 3 = 0
0 / 3 = 0
1 / 3 = 0
2 / 3 = 0
3 / 3 = 1
4 / 3 = 1
5 / 3 = 1
6 / 3 = 2

注意負數部份，例如 -4/3 得 -1 而非 -2 。

嘗試搜尋 C99 規格書，看到除法運算寫在 6.5.5 章節，但沒有查到對於商或餘數選擇的敘述

由於這種對稱性，因此前述對於正數右移 0.5 的工作，到了負數情況便改成要左移 0.5 。

回頭來看題目：









#define DIVIDE_ROUND_CLOSEST(x, divisor)                       \
    ({                                                         \
        typeof(x) __x = x;                                     \
        typeof(divisor) __d = divisor;                         \
        (((typeof(x)) -1) > 0 || ((typeof(divisor)) -1) > 0 || \
         (((__x) > 0) == ((__d) > 0)))                         \
            ? ((RRR) / (__d))                  \
            : ((SSS) / (__d));                 \
    })

因為 -1 = UINT_MAX > 0 ，所以表達式 ((typeof(x)) -1) > 0 可用於判斷該變數是否為 unsigned ；而 (((__x) > 0) == ((__d) > 0))) 檢驗兩者是否同號。

因為 C 語言會將 unsigned 和 signed 混合運算的結果轉為 unsigned ，因此只要被除數或除數中有 unsigned number ，那結果必定非負。若兩者同號那相除結果也會非負。
於是乎我們可以推測 RRR 處理

\geq 0

的結果， SSS 處理

< 0

的結果。根據之前的討論，填入：

RRR : __x + (__d >> 1)
SSS : __x - (__d >> 1)

探索 kernel 中的應用

待補

測驗 11

解釋程式碼

這裡我跳過 fls() 的部份，直接針對開平方根的部份進行解釋，其功能等同

⌊ \sqrt{x} ⌋

。





















unsigned long i_sqrt(unsigned long x)
{
    unsigned long b, m, y = 0;

    if (x <= 1)
        return x;

    m = 1UL << (fls(x) & ~1UL);
    while (m) {
        b = y + m;
        XXX;

        if (x >= b) {
            YYY;
            y += m;
        }
        ZZZ;
    }

    return y;
}

對於

x

滿足

fls (x) = n

，我可以猜到

\sqrt{x} \geq ⌊ n / 2 ⌋

，但是我不知道要怎麼算出更細節的結果。

我在 wiki 找到相關的說明，它的概念是利用

(a + b)^{2} = a^{2} + 2 a b + b^{2}

公式，想辦法在已知

a

的情況下去湊出

b

。

這題的變數和要填空的格子有點多，很難依可見的資訊判斷作答，我研究了好一陣子都很茫然，所以我直接去看答案。



















unsigned long i_sqrt(unsigned long x)
{
    unsigned long b, m, y = 0;

    if (x <= 1)
        return x;

    m = 1UL << (fls(x) & ~1UL);
    while (m) {
        b = y + m;
        y >>= 1;
        if (x >= b) {
            x -= b;
            y += m;
        }
        m >>= 2;
    }
    return y;
}

以

x = (120)_{10} = (1111000)_{2}

為例來說明：

首先fls(x) 會算出 6 ，接著 & ~1 的作用是把最低的位元改成 0 ，也就是取出小於等於 fls(x) 的最大偶數，以我們的例子會輸出 6 。最後把 m 初始化為

2^{6}

，稍候就會發現 m 變數是用來指示目前運算到第幾位。

因為 return y 所以 y 應是存開根號結果的地方。

 y =  1
    -----------------
   \/ 1 11 10 00   = x
 b =  1
     -----------
      0 11

每次它會算出 b 之後檢查 x >= b ，如果成立那 x 便可以減 b 並得該位的 y 是 1 ，否則該位 y 是 0 。

繼續 loop 之前執行第 16 行，目的是把 m 對齊下一個要運算的位。

 y =    1
    -----------------
   \/ 1 11 10 00
      1 
     -----------
      0 11       = x
 b =  1 01

2022q1 Homework3 (quiz3)

測驗 3

解釋程式碼

測驗 5

解釋程式碼

改進並實作

錯誤處理

效能改善

測驗 7

解釋程式碼

探索 kernel 中的應用

研讀《Using de Bruijn Sequences to Index a 1 in a Computer Word》

原理

編譯時期的 ilog2() 實作

測驗 10

解釋程式碼

探索 kernel 中的應用

測驗 11

解釋程式碼

Read more

讀 CS:APP 簡中版

學習 Linux 流水筆記

2022q1 Homework3 (fibdrv)

fibdrv 2

編譯時期的 `ilog2()` 實作