2024q1 Homework4 (quiz3+4)

contributed by < wu81177 >

第 3 週測驗題

測驗一

版本一當中，可以理解為使用了 binary search 來找平方根，程式先找到 msb 作為初始值，在 msb 和兩倍 msb 之間不斷二分來勘根，而因為整數型會截斷小數部分，因此 a 就會是小於等於 n 二的冪次最大值，之後若 result 加 a 後平方超過 n 就把 a 除以二查找，而由於 a 最小單位是 1 ，對於非平方數只能找到平方根的整數部分。
版本二和一的差別在於 msb 起始值的計算方式，版本一是接使用 math.h 的 log2() ，而版本二是將 n 值不斷除以二直到小於一 ( bitwise 的操作會使 n 最終變 0)，過程中用 msb 記數，一樣也能算出

$lgN$ 整數部分。
而版本三的部分，根據 Digit-by-digit calculation ，所求

$N$ 可以寫成

$N^2=(a_n+a_{n-1}+a_{n-2}+...+a_0{)}^2$ ，其中

$a_n$ 為二的

$n$ 次冪或0，可整理成:

$$\begin{array}{r}\\ 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({\displaystyle \sum _{i=1}^n{a}_i)+a_0]a_0}\\ =& 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=& a_n+a_{n-1}+...+a_m\end{array}$$
可以發現

$P_0$ 肯定是所求，同時也存在某些

$m$ 滿足

$P_m=P_0$ ，因為某些

$a_m$ 是 0 ，令：

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

上一輪底數為

$m+1$ 都是已知，同時將

$Y_m$ 拆成

$c_m,d_m$

$$\begin{gathered} c_m=P_{m+1}{2}^{m+1}, \\ {d}_m=(2^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} \end{gathered}$$
就可以藉由位元運算來計算

$c_m,d_m$，迭代求出

$c_{-1}$ 即為所求

$N$ ，了解原理後，版本三 AAAA 應為 2 ，因為

$d_m=4^m$， BBBB 應為 1 ，因為

$c_m=P_{m+1}{2}^{m+1}$

版本三

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;
}

測驗二

為了使用 bitwise 取代 / 和 % ，由於 0.1 寫成二進位會是無窮小數，需要找到一個二進位值近似 0.1 ，目標是在 0~19 範圍內的 temp 乘上它到小數點後第一位以前都正確，而測驗說明中證明了只需要檢查範圍內 temp 最大值是否達標就好，因為近似值乘上 temp 誤差會隨著 temp 等比變化，若是最大值在誤差範圍內，絕對值較小值想必也會在範圍內

證明

假設目標的最大數是 n ， l 則是比 n 還要小的非負整數。假設

$n=ab$ (

$a$ 是十位數 b 是個位數)，且

$l=cd$ (

$c$ 是十位數，

$d$ 是個位數)，則只要以

$n$ 算出來的除數在

$n$ 除以該除數後在精確度內，則

$l$ 除與該除數仍然會在精確度內，證明如下：
假設

$x$ 是近似值，

$$\begin{gathered} a.b\le \frac{n}{x}\le a.b9 \\ \Rightarrow \frac{n}{a.b9}\le x\le \frac{n}{a.b} \end{gathered}$$

1. 可見
   $x$ 的右邊是
   $10$，因此一定在精確度內。
2. $x$ 的左邊：
   
   $$\begin{gathered} c.d\le l\times \frac{a.b9}{n}\le c.d9 \\ \Rightarrow c.d\le cd\times \frac{a.b9}{ab}\le c.d9 \end{gathered}$$
   • $cd\times \frac{a.b9}{ab}$ 的左邊顯然成立
   • $cd\times \frac{a.b9}{ab}$ 的右邊：
     
     $$c.d+\frac{0.09cd}{ab}\le c.d+0.09$$
     因為
     $ab>cd$ 因此上述式子依然成立。

總之在這個例子中只需要檢查 temp 為 19 的情況就好。
0.1 近似值必須滿足一些條件，以有理數

$a/b$ 表示，其中

$b$ 為二的冪，就能用 bitwise 來運算，運算方式如下:
將

$a$ 轉為二進位

$a_n...a_2{a}_1{a}_0$ ，

$b$ 寫為

$2^N$，可以將

$tmp\ast (a/b)$ 寫為以下通式:

$$\frac{a_{n}tmp}{2^{N-n}}+...+\frac{a_{2}tmp}{2^{N-2}}+\frac{a_{1}tmp}{2^{N-1}}+\frac{a_{0}tmp}{2^N}$$
分母 2 的指數就是要右移的位數，以

$13/8$ 為例，

$N=3$，

$a=(1101{)}_2$ ，可以寫為， d 是為了保留右移被捨棄的位元 :

d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;

((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2)

從上述通式能看出若要找 0.1 近似值，

$a$ 取

$b/10$ 四捨五入至整數位誤差最小，但為了滿足小數點後第一位以前正確，

$a$ 應該取

$b/10$ 無條件進位值，這樣誤差才會是正的，而

$N$ 越大誤差也就越小，如圖，橫軸為 N 縱軸為 0.1近似值*19的誤差:

可以看到當

$N>=7$ 就能滿足要求，然而 a 實際上可為小數。
接下來看題目實作的函式:

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) 是將 in 減去 1/4 存為 x ，但不懂為何要 in | 1，而第二行 q 約為

$in(\frac{3}{4}\ast \frac{1}{16}+\frac{3}{4})$ 等於

$0.796875in$ ，之後的5行讓 q 超級接近 0.8 ，這時我才發現它想讓近似值的分子為 0.8，我在前面的討論沒有意識到分子可為小數。
因為分子是 0.8 ，分母應該要為 8 ，所以依照上面的規則 CCCC 應為 3 。
q & ~0x7 的效果是將後面三位設成0，相當於

$0.8in/8\ast 8=8\ast div$ ，因此後面那項應該要是

$2\ast div$ 加起來才會是

$10\ast div$， DDDD 應為 1 ，我的疑惑是不知為何 (q & ~0x7) + (*div << 1) 兩項不統一為其中一種寫法就好，不知道這樣有什麼好處。

測驗三

版本一的做法是逐位右移，相當於是一直除以二並且計數，直到等於0就能求出二為底的對數
版本二的作法也一樣，差別是對於很大的數檢查如果大於某數就一次右移很多位，這樣計算的次數比起版本一會少很多，所以 AAAA 應為

$2^{16}=65536$，BBBB 應為

$2^8=256$ ， CCCC 為

$2^4=16$

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;
}

版本三使用了 __builtin_clz 來改寫 ，他是 GCC 中的 extention，功能是會計算括號內數左邊0的個數
根據 GCC 手冊

Built-in Function: int __builtin_clz (unsigned int x)
Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.

所以括號中應該要填 v ，但為了避免 v 為 0 時 __builtin_clz 沒定義的問題， DDDD 應該要填 v|1 這樣最小值就是 1 不是 0

測驗四

Exponentially Weighted Moving Average EWMA 是一種 moving average ，越新的樣本權重越大，隨著時間增加權重等比下降，但定義上永遠不會變0，所以和 simple moving average 不同，EWMA 是一種 IIR filter，然而得到新的點不用經過太多計算，因為可以用遞迴關係表示:

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

$0<\alpha \le 1$ ，如同其他 ME ，它是個低通濾波器，可以觀察到若

$\alpha$ 值越大，過去的樣本影響越小，因此留下的高頻成分越多，若

$\alpha$ 為最大值 1 ，EWME 將不發揮作用
題目實作程式碼如下：

#include <assert.h>

/* Exponentially weighted moving average (EWMA) */
struct ewma {
    unsigned long internal;
    unsigned long factor;
    unsigned long weight;
};

/* This section contains universal functions designed for computing the EWMA.
 * It maintains a structure housing the EWMA parameters alongside a scaled
 * internal representation of the average value, aimed at minimizing rounding
 * errors. Initialization of the scaling factor and the exponential weight (also
 * known as the decay rate) is achieved through the ewma_init(). Direct access                                                                                
 * to the structure is discouraged; instead, interaction should occur
 * exclusively via the designated helper functions.
 */

/**
 * ewma_init() - Initialize EWMA parameters
 * @avg: Average structure
 * @factor: Scaling factor for the internal representation of the value. The
 *     highest achievable average is determined by the formula
 *     ULONG_MAX / (factor * weight). To optimize performance, the scaling
 *     factor must be a power of 2.
 * @weight: Exponential weight, also known as the decay rate, dictates the rate
 *     at which older values' influence diminishes. For efficiency, the weight
 *     must be set to a power of 2.
 *
 * Initialize the EWMA parameters for a given struct ewma @avg.
 */
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_add() - Exponentially weighted moving average (EWMA)
 * @avg: Average structure
 * @val: Current value
 *
 * Add a sample to the average.
 */
struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
    avg->internal = avg->internal
                        ? (((avg->internal << EEEE) - avg->internal) +
                           (val << FFFF)) >> avg->weight
                        : (val << avg->factor);
    return avg;
}

/**
 * ewma_read() - Get average value
 * @avg: Average structure
 *
 * Returns the average value held in @avg.
 */
unsigned long ewma_read(const struct ewma *avg)
{
    return avg->internal >> avg->factor;
}

一個 ewma 結構代表一個 EWMA 濾波器， internal 表示要輸出的樣本值， factor 是一個放大因子用來放大輸入樣本降低誤差， weight 則和

$\alpha$ 有關。
函式部分有 ewma_init ，用來初始化一個 EWMA 結構，它先檢查 weight 和 factor 是不是二的冪，然後將它們取

$lg$ ，方便後面用平移運算取代乘法，還有將 internal 設為 0 。
再來是 ewma_add 函式:

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

函式首先判斷 internal 是不是非 0，也就是在判斷是不是初始狀態，如果是就將 val 放大取

$lg$ 前的 factor 倍; 如果不是初始狀態則要進行前面提到的遞迴運算

$\alpha Y_t+(1-\alpha )S_{t-1}$ ，(avg->internal << EEEE) - avg->internal 可以對應到

$(1-\alpha )S_{t-1}$ 的部分，整段可以寫成(忽略 factor)

$((Y_t\ast weight-Y_t)+Y_{t-1})/weight$ ，因此 EEEE 應為 avg->weight ，FFFF 應為 avg->factor。
從這邊也就能看的出來 weight 和

$\alpha$ 應為倒數關係。

測驗五

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 > GGG) + 1;       
}

函式一開始的 x-- 是為了破壞 x 剛好是二的冪的情況，函式中會達到計算

$lg$ 並且捨棄小數的效果，最後 return 的地方會再加上一達成向上取整的功能，如果沒有 x-- 二的冪的輸出就會多一。
再來的原理和測驗版本二類似，只是寫法更緊湊，需要好好思考; 由於 x 是 32 bits 的無號整數型，因此可以不使用迴圈逐步將檢查範圍覆蓋到八個 byte 就好， 依序檢查 x 是否大於

$2^{2^4}$,

$2^{2^3}$,

$2^{2^2}$,如果是就除掉，用 or operator 將結果依序更新到 r ，最後在 return 之前已經檢查了

$2^4+2^3+2^2+2=30$ 個 bits ，還剩兩個要檢查，於是我代入函式傳入時 x=1~4 推定 GGG 為 1 。

第 4 週測驗題

測驗一

Popcount 用來計算二進位數值中有幾個 1 ，實作版本一如下:

unsigned popcount_naive(unsigned v)
{
    unsigned n = 0;
    while (v)
        v &= (v - 1), n = -(~n);
    return n;
}

v &= (v - 1) 的效果是將最低位的 1 變成 0 ， n = -(~n) 的作用相當於 n++ ，因為在二補數系統中負號是先取反再加一，所以 －(~n) = ~(~n)+1 = n+1， n 會在迴圈中計數直到 v 為 0 終止，但這樣程式不是常數時間，可能有資安疑慮，因此可以改進成版本二:

unsigned popcount_branchless(unsigned v)
{
    unsigned n;
    n = (v >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;

    v = (v + (v >> 4)) & 0x0F0F0F0F;
    v *= 0x01010101;                                     

    return v >> 24;
}

這個方式是將每 4 個 bits 平行計算， (v >> 1) & 0x77777777 作用是先除以二，且為了保持每 4bits 的獨立性，將每四位第一位改成 0 。
考慮到計算 n bits popcount 公式如下:

$$nbitspopcount(x)=x-\lfloor \frac{x}{2^1}\rfloor -\lfloor \frac{x}{2^2}\rfloor -...-\lfloor \frac{x}{2^n}\rfloor$$
所以 4 bits 就是

$x-\lfloor \frac{x}{2^1}\rfloor -\lfloor \frac{x}{2^2}\rfloor -\lfloor \frac{x}{2^3}\rfloor$，第 4~9 行用遞迴的方式達成此式。
但這樣還不是最終結果，還需要把 v 每 4bits 分開再加起來，第 11 行會讓第奇數個 nibble 為原本的每兩個相加，第偶數個變 0 ，因為原本每個 nibble 最大是 4 所以不會溢位，然後再乘以 0x01010101 再右移 24 bits 寫成直式就能看出會是所有 nibble 相加的效果。

Hamming distance

兩數的 Hamming distance 可以理解為將一個數變成另一個數需要翻轉的位元個數，實作上可以先做 bitwise XOR 再用 popcount 來數:

int hammingDistance(int x, int y)
{
    return __builtin_popcount(x ^ y);
}

其中__builtin_popcount 是 GCC 內建函式。
考慮 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;
}

題目要求加總陣列每兩個數的 Hamming distance ，for loop 算的每一組會重複配對一次，所以結果要除以二， AAAA 為 2 。

測驗二

若

$a=b(modm)$ ,

$c=d(modm)$ 則

$ac=bd(modm)$ ,

$a+c=b+d(modm)$ 。
若除數為 3 ，則套用上面第二個性質可以推導出

$$2^k=\{\begin{array}{l}1(mod3)\text{,k is even }\\ -1\text{, k is odd }\end{array}$$
考慮正整數

$n$ 以二進位表示 :

$b_{n-1}{b}_{n-2}...b_0$，

$n=b_{n-1}{2}^{n-1}+b_{n-2}{2}^{n-2}+...+b_0{2}^0$ ，套用上面性質得到

$n={\mathrm{\Sigma }}_{i=0}^{n-1}{b}_i(-1{)}^i(mod3)$ ，求和的部分相當於奇數位和減偶數位和: popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA) ，

$mod3$ 的部分則是重複操作前面直到 n = 0~2 ，類似的方法也可以推廣到任何除數為

$2^k\pm 1$ 的時候。
考慮 tictactoe.c 中 mod 7 的部分:

static inline int mod7(uint32_t x)
{
    x = (x >> CCCC) + (x & UINT32_C(0x7FFF));
    /* Take reminder as (mod 8) by mul/shift. Since the multiplier
     * was calculated using ceil() instead of floor(), it skips the
     * value '7' properly.
     *    M <- ceil(ldexp(8/7, 29))
     */
    return (int) ((x * UINT32_C(0x24924925)) >> 29);
}

不同於之前的討論，這裡是用 Hacker's Delight 中 14 頁的方法，可以避免使用 popcount ，對於每個

$2^k\pm 1$ 的除數，都可以從某些地方拆分成前後兩項再加起來而不影響求模結果，在這題由於後面是 & 0x7FFF 是 15 個 1，可以推測 CCCC 是 15。

/* Determine if the tic-tac-toe board is in a winning state. */
static inline uint32_t is_win(uint32_t player_board)
{
    return (player_board + BBBB) & 0x88888888;
}

如果想知道 BBBB 的答案就要先知道如何描述棋盤，它使用一個 32 bits 無號整數來描述棋盤狀態，經過比對 move_masks 可以發現，每四個 bits 描述一種連線，總共有八條剛好 32 bits，第一個 bit 為 0，第 2~4 分別代表那條線上三個位置有沒有值，所以 BBBB 是 0x11111111，因為這樣只要存在一條連線那條線的四個 bit 就會是 0111 + 0001 然後進位，接著 & 1000 看有沒有進位就完成檢查。

測驗三

AVL tree 為了極致的降低搜尋速度，限制任何左右子樹高度差不能超過一，因此常常需要透過旋轉 (LL, LR, RR, RL)來平衡子樹，過程會消耗許多資源，因此 XTree 被設計出來，它比較不會消耗資源去平衡，但樹的分布較不平衡，平衡能力和搜尋時間介於 red-black tree 和 AVL tree 之間。
接著來看移除節點的函式:

static void __xt_remove(struct xt_node **root, struct xt_node *del)
{
    if (xt_right(del)) {
        struct xt_node *least = xt_first(xt_right(del));
        if (del == *root)
            *root = least;

        xt_replace_right(del, AAAA);
        xt_update(root, BBBB);
        return;
    }

    if (xt_left(del)) {
        struct xt_node *most = xt_last(xt_left(del));
        if (del == *root)
            *root = most;

        xt_replace_left(del, CCCC);
        xt_update(root, DDDD);
        return;
    }

    if (del == *root) {
        *root = 0;
        return;
    }

    /* empty node */
    struct xt_node *parent = xt_parent(del);

    if (xt_left(parent) == del)
        xt_left(parent) = 0;
    else
        xt_right(parent) = 0;

    xt_update(EEEE, FFFF);
}

根據註解，如果要刪除的節點有右子節點，則待刪除的節點會先替換為右子樹中遇到的第一個節點，因此 AAAA 為 least ，它在前面被設為 xt_first ，接著會對新接上節點的右子節點做 xt_update 來更新 hint ，並且決定要不要平衡，因此 BBBB 是 xt_right(least) ，同理 CCCC 為 most ， DDDD 為 xt_left(most)。
如果要刪除的節點是 leaf 的情況可以直接刪除，但就要檢查它的親代，確保平衡狀況，所以 EEEE 是 root ， FFFF 是 parent 。