2024q1 Homework4 (quiz3+4)

contributed by < yuyuan0625 >

第三週測驗

2024q1 第 3 週測驗題

測驗一

i_sqrt 版本三利用 Digit-by-digit calculation，將要開平方的數

$N$ 拆成 2 的冪相加。

例如

$N=(19{)}^2$ 就會轉換為

$N=(10011{)}_2$

$$N^2=(a_n+a_{n-1}+a_{n-2}+...+a_0{)}^2，a_m=2^{m}or{a}_m=0$$

根據

$(x+y{)}^2=x^2+2xy+y^2$ 的規律展開為：

$N^2={a_n}^2+[2a_n+a_{n-1}]a_{n-1}+[2(a_n+a_{n-1}+a_{n-2})]a_{n-2}]+...+[2(\sum _{i=1}^n{a}_i)+a_0]a_0$

假設

$P_m=a_n+a_{n-1}+...+a_m$，

$P_0=a_n+a_{n-1}+...+a_0$ 即為所求平方根

$N$
則原式可代換為

$N^2={a_n}^2+[2P_n+a_{n-1}]a_{n-1}+[2P_{n-1}+a_{n-2}]a_{n-2}+...+[2P_1+a_0]a_0$
而

$P_m=P_{m+1}+a_m$

我們需要從

$m=n$ 一直往右測到

$m=0$ ，每一輪則透過檢查

${P_m}^2\le N^2$ 是否成立來求得

$a_m=2^{m}or{a}_m=0$。
但如果每一輪都要計算

$P^2-N^2$ 來求得

$a_m$ 的運算成本太高，因此使用上一輪的差值

$X_{m+1}$ 減去

$Y_m$ 得到

$X_m$

• $X_m=N^2-P_m^2=X_{m+1}-Y_m$
• $Y_m=P_m^2-P_{m+1}^2=(P_{m+1}+a_m{)}^2-P_{m+1}^2=2P_{m+1}{a}_m+a_m^2$

就可以透過紀錄

$P_{m+1}$ 來取代

$P_m^2$ 的計算。
接著將

$Y_m$ 拆成

$c_m，d_m$ 兩個部分，

• $c_m=2P_{m+1}{a}_m$
• $d_m={a_m}^2$
• $Y_m=\{\begin{array}{ll}{c}_m+d_m& \text{if }{a}_m=2^m\\ 0& \text{if }{a}_m=0\end{array}$

藉由位元運算從

$c_m，d_m$ 推出下一輪

$c_{m-1}，d_{m-1}$

• $c_{m-1}=P_m{2}^m=(P_{m+1}+a_m)2^m=P_{m+1}{2}^m+a_m{2}^m=\{\begin{array}{ll}{c}_m/2+d_m& \text{if }{a}_m=2^m\\ c_m/2& \text{if }{a}_m=0\end{array}$
• $d_m=a_{m-1}^2=(2^{m-1}{)}^2={\displaystyle \frac{(2^m{)}^2}{2^2}}={\displaystyle \frac{d_m}{4}}$

結合上述方法，求

$a_n+a_{n-1}+...+a_0$，假設從

$a_n$ 開始往下測試，所以

• $X_{n-1}=N$
• $c_n=0$
• $c_{-1}=P_0{a}^0=P_0=a_n+a_{n-1}+...+a_0$ 即為所求
  $N$

int i_sqrt(int x)
{
    if (x <= 1) /* Assume x is always positive */
        return x;

    int z = 0;
    for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
        int b = z + m;
        z >>= 1;
        if (x >= b)
            x -= b, z += m;               
    }
    return z;
}

首先要確認 x 不為 0 或負數(涉及虛數的處理，不討論)

if (x <= 1) /* Assume x is always positive */
    return x;

int m 即為

$d_m$ ，由於首項

$d_n=(2^n{)}^2$，因此要先利用 __builtin_clz 找到 x 最高有效位元前面有幾個 0， x 為 0 時，未定義。

b 為

$Y_m$ 而 z 即為

$c_m$ ，由上述推導

$c_m$ 每一輪都要除以2， 因此 z 要向右位移 1 位
那因為

$d_{m-1}={\displaystyle \frac{d_m}{4}}$ 因此 m 每次迴圈往右位移 2 位

嘗試用 ffs / fls 取代 __builtin_clz

__builtin_clz(x) 函式回傳 x 的最高有效位元前面的連續 0 位元數量，因此 31 - __builtin_clz 即為最高有效位元位置。
而 fls()也是相同概念，但是其索引值是由 1 開始計算，因此需要將 fls() - 1 才能達到和 __builtin_clz(x) 一樣的效果。

因此可以將程式碼修改：

int i_sqrt(int x)
{
    if (x <= 1) /* Assume x is always positive */
        return x;

    int z = 0;
+   for (int m = m = 1U << ((fls() - 1) & ~1U); m; m >>= 2) {
-   for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
        int b = z + m;
        z >>= 1;
        if (x >= b)
            x -= b, z += m;               
    }
    return z;
}

測驗二

static void str_add(char *b, char *a, char *res, size_t size)
{
    int carry = 0;

    for (int i = 0; i < size; i++) {
        int tmp = (b[i] - '0') + (a[i] - '0') + carry;
        carry = tmp / 10;
        tmp = tmp % 10;
        res[i] = tmp + '0';
    }
}

餘式定理：
被除數 ＝ （商*除數）+ 餘數

對應程式碼：

carry = tmp / 10;
tmp = tmp - carry * 10;

若採用 bitwise operation 來實作上述除法，會因為

$10$ 包含

$5$ 這個因數無法全用

$2$ 的冪項來表示，進而產生誤差。

由上述程式碼可發現 tmp 的值不會大於 19

• (b[i] - '0') 和 (a[i] - '0') 皆為 0~9的整數
• carry 為進位值，範圍由 0~1
• tmp 最大為 9 + 9 + 1 = 19

$1.9\le {\displaystyle \frac{19}{x}}\le 1.99\Rightarrow 9.55\le x\le 10$

找除數的方法是使用 bitwise operation

$\frac{2^N}{a}$ 找到介於

$9.55\le x\le 10$ 的除數，若被除數為

$n$ ，商式可以寫成

$\frac{an}{2^N}$ ，因此只需查看

$2^N$ 再配對適合的

$a$ 即可。
其中，

$2^N=128，a=13，\frac{128}{13}\approx 9.84$ 為一個可用的除數，由於 13 可以拆成

$13=8+4+1=2^3+2^2+2^0$ ，因此範例程式中透過 (tmp >> 3) + (tmp >> 1) + tmp 得到

$\frac{tmp}{8}+\frac{4tmp}{8}+tmp=\frac{13tmp}{8}$ ，再將此式乘上 8 (向左位移 3 bits) 即可得到

$13tmp$ ，只要再將其除以

$128$ (

$2^7$) 即可得到目標商式

$\frac{13tmp}{2^7}$ 。

(((q << 2) + q) << 1) 這部分是將 q * 10 透過 (q*4 + q) * 2 實作

包裝後函式：

#include <stdint.h>
void divmod_10(uint32_t in, uint32_t *div, uint32_t *mod)
{
    uint32_t x = (in | 1) - (in >> 2); /* div = in/10 ==> div = 0.75*in/8 */
    uint32_t q = (x >> 4) + x;
    x = q;
    q = (q >> 8) + x;
    q = (q >> 8) + x;
    q = (q >> 8) + x;
    q = (q >> 8) + x;

    *div = (q >> CCCC);
    *mod = in - ((q & ~0x7) + (*div << DDDD));   
}

uint32_t x = (in | 1) - (in >> 2) 將

$x=\frac{3}{4}in$
再透過 uint32_t q = (x >> 4) + x ，

$q=\frac{3}{4\ast 2^4}in+\frac{3}{4}in=\frac{102}{128}in$ ，而其中

$\frac{102}{128}\approx 0.797\approx \frac{8}{10}$

後續再用 q = (q >> 8) + x 對 q 做逼近增加精度更靠近

$\frac{8}{10}$

*div = (q >> CCCC) 是計算商，因為我們前面計算的 q 為

$\frac{8}{10}in$ ，所以要再除以 8 才會等於

$\frac{1}{10}in$ ， 故 q >> 3

*mod = in - ((q & ~0x7) + (*div << 1)) 是計算餘數，其中的 q & ~0x7 等於是 *div << 3 ，因此程式碼可以轉換成 *mod = in - (*(div << 3) + (*div << DDDD))
根據餘數定理

$\mathrm{餘}\mathrm{數}=\mathrm{被}\mathrm{除}\mathrm{數}(in)-\mathrm{商}(div)\times 10\mathrm{除}\mathrm{數}(10)$

$\mathrm{商}\times 8+\mathrm{商}\times 2=\mathrm{商}\times (8+2)=\mathrm{商}\times 10$ ，故後半部應為 *(div << 3) + (*div << 1)

測驗三

ilog2 計算以 2 為底的對數，且其輸入和輸出皆為整數。

版本一

int ilog2(int i)
{
    int log = -1;
    while (i) {
        i >>= 1;
        log++;
    }
    return log;
}

從最低位元往高位元尋找最高的有效位元位置，最初將 log 設為 -1 ，讓函式傳入 0 時輸出 -1 。每一次迴圈內 i >>= 1 ，相當於 i 除以 2 ，並且將 log 的值加 1 。當 i == 0 時迴圈停止， log 即為所求。

版本二

static size_t ilog2(size_t i)
{
    size_t result = 0;
    while (i >= AAAA) {
        result += 16;
        i >>= 16;
    }
    while (i >= BBBB) {
        result += 8;
        i >>= 8;
    }
    while (i >= CCCC) {
        result += 4;
        i >>= 4;
    }
    while (i >= 2) {
        result += 1;
        i >>= 1;
    }
    return result;
}

依照 i 的大小，提供四種右移的方式，只要 i 大於

$2^k$，就一次右移 k 個位元，讓迴圈的執行次數可以小於 n。
AAAA 即為

$2^{16}=65536$
BBBB 等於

$2^8=256$
CCCC 等於

$2^4=16$

版本三

int ilog32(uint32_t v)
{
    return (31 - __builtin_clz(DDDD));
}

利用 GNU extension __builtin_clz ，找出最高有效位元前面 0 的數量，因此 31 - __builtin_clz(v | 1) 即為最高有效位元的位置。另外，__builtin_clz 輸入若是 0 則無定義，所以需使用 v | 1 確保輸入不為 0 。

測驗四

Exponentially Weighted Moving Average (EWMA; 指數加權移動平均) 是種取平均的統計手法，並且使經過時間越久的歷史資料的權重也會越低，以下為 EWMA 的數學定義:

$$S_t=\{\begin{array}{l}{Y}_0,t=0\\ \alpha Y_t+(1-\alpha )\cdot S_{t-1},t>0\end{array}$$

• $\alpha$ 表示歷史資料加權降低的程度，介在 0 ～ 1 之間，越高的
  $\alpha$ 會使歷史資料減少的越快
• $Y_t$ 表示在時間
  $t$ 時的資料點
• $S_t$ 表示在時間
  $t$ 時計算出的 EWMA

struct ewma {
    unsigned long internal;
    unsigned long factor;
    unsigned long weight;
};

首先看 ewma 結構，

• internal 儲存平均值，也就是
  $S_t$
• factor 為 scalinf factor
• weight 為 decay rate
  $\alpha$

在 ewma_init() 註解有提到 internal 可記錄最大值的計算公式為 ULONG_MAX / (factor * weight) ，其中 factor 是因為做 scaling 提高數值精度變相會犧牲可容納的對大值。

static inline int is_power_of_2(unsigned long n)
{
    return (n != 0 && ((n & (n - 1)) == 0));
}

void ewma_init(struct ewma *avg, unsigned long factor, unsigned long weight)
{
    if (!is_power_of_2(weight) || !is_power_of_2(factor))
        assert(0 && "weight and factor have to be a power of two!");

    avg->weight = ilog2(weight);
    avg->factor = ilog2(factor);
    avg->internal = 0;
}

接著看到 ewma_init 函式，用來初始化 ewma 結構得初始值。其中使用 is_power_of_2 來確保 weight 跟 factor 是 2 的冪，因為 internal 是定點數若是與 2 的冪做乘除法，可以用位元運算。
is_power_of_2 是運用當數值 n 是 2 的冪時，n 跟 n-1 在二進位時，不會有相同的位數的特性，因此做 ＆ 運算後值會為 0 。
例如：

$8={1000}_2,7={0111}_2$ ，兩者做 & 為

${0000}_2$

struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
    avg->internal = avg->internal
                        ? (((avg->internal << avg->weight) - avg->internal) +
                           (val << avg->factor)) >> avg->weight
                        : (val << avg->factor);
    return avg;
}

$$S_t=\{\begin{array}{l}{Y}_0,t=0\\ \alpha Y_t+(1-\alpha )\cdot S_{t-1},t>0\end{array}$$

最後看 ewma_add，當 avg->internal 為 0 時，不需考慮歷史平均，直接將 val作為輸入，但因為有 scaling 所以還需要 val << avg->factor

原先我看 (avg->internal << avg->weight) - avg->internal) 這邊操作好像和想像中的

$1-\alpha$ 剛好相反有點困惑，後來參考SHChang-Anderson 同學的筆記後才看懂，此程式將

$\alpha$ 設為

$\frac{1}{2^{avg->weight}}$，所以程式是先將

$$\alpha Y_t+(1-\alpha )\cdot S_{t-1}=[\alpha Y_t+(1-\alpha )\cdot S_{t-1}]\times 2^{avg->weight}\times \frac{1}{2^{avg->weight}}$$
將

$2^{avg->weight}$ 乘進去原始公式，再簡單地移項方便對證程式碼，可以得到

$$[2^{avg->weight}\cdot S_{t-1}-S_{t-1}+Y_t]\times \frac{1}{2^{avg->weight}}$$
其中先乘以

$2^{avg->weight}$ 再乘

$\frac{1}{2^{avg->weight}}$ 是為了提高精度，程式即是使用上述的公式做計算，並且在每次輸入 val 都會 scaling

測驗五

以下程式碼可計算

$\lceil lo{g}_2(x)\rceil$ ，對於傳入的參數

$x$ ，回傳最小的整數

$n$，滿足

$x\le 2^n$ 。

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;

    x--;
    r = (x > 0xFFFF) << 4;
    x >>= r;

    shift = (x > 0xFF) << 3;
    x >>= shift;
    r |= shift;

    shift = (x > 0xF) << 2;
    x >>= shift;
    r |= shift;

    shift = (x > 0x3) << 1;
    x >>= shift;
    
    return (r | shift | x > 1) + 1;  

函式剛開始時會將 x--，因為是取 ceil ，若 x 剛好為 2 的冪次時不需要進位，所以就直接將 x-- 讓 x 變小後再進位即可。
r = (x > 0xFFFF) << 4 這邊和測驗 3 是相同的概念 0xFFFF 等於

$2^{16}$，若x > 0xFFFF，r = 1 << 4 等於 16 ，接著將 x 向右位移16 位。
而此處 r | shift 等效於 r + shift （result +=）。
因此位移後程式碼將持續累加位移量，以找到最高位元位置。

最後 return (r | shift | x > 1) + 1 的部分我們分開看，其中 r | shift 是將前一部分的 r |= shift 合併進來。
對照測驗三程式碼會發現，此函式少了一個判斷：

while (i >= 2) {
    result += 1;
    i >>= 1;
}

因此後半部份應為 x > 1 是來處理 x= 0x2 的情況，最後再加上 1 達到取上界 (ceil) 的作用。

改進程式碼

當 x=0 時會因為減一變成 0xFFFFFFFF ，和預期結果相同。
使用!!(x) 將整數輸入結果控制在 0 跟 1 ， 那麼就可以將 x-- 變為 x = x - !!x ，當 x > 0 時減一，而當 x = 0 時則不變。

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;

-   x--;
+   x = x - !!x
    r = (x > 0xFFFF) << 4;
    x >>= r;

第四週作業

2024q1 第 4 週測驗題

測驗一

針對 LeetCode 477. Total Hamming Distance，考慮以下程式碼:

int totalHammingDistance(int* nums, int numsSize)
{
    int total = 0;;
    for (int i = 0;i < numsSize;i++)
        for (int j = 0; j < numsSize;j++)
            total += __builtin_popcount(nums[i] ^ nums[j]); 
    return total >> AAAA;
}

上述程式會計算兩兩數字的漢名距離，如果有兩個數 a 和 b 就會計算 ab 和 ba 的漢明距離的總和。 由此可知最後的輸出應該要再除以 2 ，因此 total >> 1

測驗二

Remainder by Summing digits，若除數符合

$2^k\pm 1$ ，則可以運用以下手法來達成不使用任何除法就算出某數除以另一個數的餘數。

若

$a\equiv b(modm)$ 且

$c\equiv d(modm)$ ， 則

$a+c\equiv b+d(modm)$ 且

$ac\equiv bd(modm)$

以除數 3 為例，

$1\equiv 1(mod3)$ 且

$2\equiv -1(mod3)$

$2^k\equiv \{\begin{array}{l}1(mod3),keven\\ -1(mod3),kodd\end{array}$

若 n 的二進位表示為

$b_{n-1}{b}_{n-2}{b}_{n-3}...b_1{b}_0$

$n=b_{n-1}\cdot 2^{n-1}+b_{n-2}\cdot 2^{n-2}+b_{n-3}\cdot 2^{n-3}+...+b_1\cdot 2^1+b_0\equiv b_{n-1}+\dots -b_3+b_2-b_1+b_0(mod3)$

位元和可以利用 population count 這類的函式來得到

• $5={0101}_2$ ，因此 0x55555555 即為所有奇數位
• $A={1010}_2$ ，因此 0xFFFFFFFF 即為所有偶數位

因此寫成程式的話可以將上式表示為 n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)

接著，使用以下定理進行化簡：

$popcount(x\wedge \bar{m})-popcount(x\wedge m)=popcount(x\oplus m)-popcount(m)$

因此，n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA) 可以寫為 n = popcount(n ^ 0xAAAAAAAA) - 16。

但此作法的計算結果會介於 -16 至 16 之間，若希望餘數為正就必須再加上一個 3 的倍數來確保餘數為正。

文中的例子是加上 39 。範例程式如下

int mod3(unsigned n)
{
    n = popcount(n ^ 0xAAAAAAAA) + 23;
    n = popcount(n ^ 0x2A) - 3;
    return n + ((n >> 31) & 3);
}

《Hacker's Delight》中說明為何要選 39 ：

We want to apply this transformation again, until n is in the range 0 to 2, if possible. But it is best to avoid producing a negative value of n, because the sign bit would not be treated properly on the next round. A negative value can be avoided by adding a sufficiently large multiple of 3 to n. Bonzini’s code, shown in Figure 10–21, increases the constant by 39. This is larger than necessary to make n nonnegative, but it causes n to range from –3 to 2 (rather than –3 to 3) after the second round of reduction. This simplifies the code on the return statement, which is adding 3 if n is negative. The function executes in 11 instructions, counting two to load the large constant.

這比必要的大，但它使得 n 在第二輪縮減後的範圍是 -3 到 2（而不是 -3 到 3）。這簡化了返回語句中的程式碼，如果 n 為負，則添加 3。

另一種變形是利用 lookup table，將 0 到 31 mod 3 的結果

int mod3(unsigned n)
{
    static char table[33] = {2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1 };
    n = popcount(n ^ 0xAAAAAAAA);
    return table[n];
}

井字遊戲程式碼

測驗三

此程式紀錄 3x3 棋盤中可能的 8 條線，而不是記錄傳統的九宮格方式。因為棋盤大小為 3x3 故總操作數只會有 9 次。每次從陣列中選出一個可以操作的位置，選擇完後將本次選擇的位置，從陣列中移除。並將其給 board | move_masks[move] ，隨後去檢查該名玩家是否勝利。

	1	2	3
1	0	1	2
2	3	4	5
3	6	7	8

假設九宮格由左至由，由上而下分別編號為 0~8 號。
8 條線分別為: (0,1,2), (3,4,5), (6,7,8), (0,3,6), (1,4,7), (2,5,8), (0,4,8), (2,4,6)，分別對應到下方 move_masks 陣列中的每個元素 16 進位由左至右的 8 個數值。

static const uint32_t move_masks[9] = {
    0x40040040, 0x20004000, 0x10000404, 0x04020000, 0x02002022,
    0x01000200, 0x00410001, 0x00201000, 0x00100110,
};

move_masks 陣列中的每個元素代表了在將棋子放置到特定位置後，對於連線狀態的影響。每個元素的二進位表示描述了在該位置放置棋子後，連線狀態的改變。

勝利的條件判斷為 player_board 以四個位元為單位出現 0111 即判斷該玩家獲勝，可以看到程式碼將 (player_board + BBBB) 與 0x88888888 做

$and$ 運算。由此可知，當出現 0111 時需要將棋結果轉為 1000 ，而將 0111 + 1 即可達成此效果，因此 BBBB 應填入 0x11111111 。