2024q1 Homework4 (quiz3+4)

contributed by yu-hsiennn

quiz 3

題目

測驗 1

假定

$P_m$ 為欲所求之平方根

$$\{\begin{array}{l}{N}^2=(a_n+a_{n-1}+a_{n-2}+...+a_0{)}^2,a_m=2^{m}or{a}_m=0\\ P_m=a_n+a_{n-1}+...+a_m\end{array}$$
經展開後，化簡可得

$$N^2=\sum _{i=0}^n{a}_i^2+2\sum _{i=0}^{n-1}{P}_{i+1}\cdot a_i$$
令

$X_m$ 為

$N^2-P_m^2$，且此輪與上一輪差設為

$Y_m$，即

$$\{\begin{array}{l}{X}_m=N^2-P_m^2\\ Y_m=X_{m+1}-X_m=P_m^2-P_{m+1}^2=2P_{m+1}{a}_m+a_m^2\end{array}$$
我們再將式子利用

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

再觀察下一輪的變化

$$\{\begin{array}{l}{c}_{m-1}=2P_m{a}_{m-1}=2(P_{m+1}+a_m)a_{m-1}=P_{m+1}{a}_m+a_m^2={\displaystyle \frac{c_m}{2}}+d_m\\ d_{m-1}={\displaystyle \frac{d_m}{4}}\end{array}$$
因此

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) {
        // m: d_m, z: c_m
        int b = z + m;
        z >>= 1;
        if (x >= b)
            x -= b, z += m;               
    }
    return z;
}

改用第二周測驗 2 的 __ffs

-   for (int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL); m; m >>= 2)
+   for (int m = 1UL << ((31 - __ffs(x)) & ~1UL); m; m >>= 2)

疑惑

$$\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}$$

在一開始看到這個時，完全無法理解是怎麼跑出這些等號的。
後來才發現，其實有段開，只是 HackMD 沒有顯示出來

測驗 2

因為乘法及除法的運算成本很高，故想利用位移及加減的方式去達到相同的效果
此處數字範圍為 0~19，只要這一段的計算有達到預期即可
題目是用 N = 128 來逼近 1/10

$1/128+1/32+1/16=(1+4+8)/128=13/128\approx 1/10$
也就是說可以

/* 13n / 128 */
q = (n + (n << 2) + (n << 3)) >> 7;

但這邊會忽略掉一些值，如果要精確的話可以把他加回去

d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
q = (n + (n << 2) + (n << 3) + d0 + d1 + d2) >> 7;

不過在範圍 0~19 其實並不會影響結果
如果範圍是 0~19 ，這邊做了一些測試

// test1: n / 16 + n / 32 + 1 / 8
q = ((n << 1) + n + 4) >> 5;
// test2: n / 16 + n / 32 + 1 / 16
q = ((n << 1) + n + 2) >> 5;

結果為

q: 0, r: 0
q: 0, r: 1
q: 0, r: 2
q: 0, r: 3
q: 0, r: 4
q: 0, r: 5
q: 0, r: 6
q: 0, r: 7
q: 0, r: 8
q: 0, r: 9
q: 1, r: 0
q: 1, r: 1
q: 1, r: 2
q: 1, r: 3
q: 1, r: 4
q: 1, r: 5
q: 1, r: 6
q: 1, r: 7
q: 1, r: 8
q: 1, r: 9

而題目

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 >> 3);
    *mod = in - ((q & ~0x7) + (*div << 1));   
}

可以看到這邊的作法是

$$x={\displaystyle \frac{3}{4}}n+{\displaystyle \frac{3}{64}}n=0.796875n\approx {\displaystyle \frac{8}{10}}n$$
所以將 q 右移 3 可得

${\displaystyle \frac{8}{10}}n\ast {\displaystyle \frac{1}{8}}n={\displaystyle \frac{1}{10}}n$
而最後是在計算 *mod = in - *div * 10

測驗 3

利用 4 層迴圈來找尋以 2 為底的對數

• 第一個迴圈作用於處理大於等於 65536 的數值，直到小於才離開此迴圈
• 第二個迴圈作用於處理大於等於 256 的數值，直到小於才離開此迴圈
• 第三個迴圈作用於處理大於等於 16 的數值，直到小於才離開此迴圈
• 第四個迴圈作用於處理大於等於 2 的數值，直到小於才離開此迴圈，並回傳結果

舉個例子:
n = 2049
2049 < 65536, (不執行第一個迴圈)
2049 >= 256, n: 8, result: 8
8 < 16, (不執行第三個迴圈)
8 >= 2, n: 4, result: 9
4 >= 2, n: 2, result: 10
2 >= 2, n: 1, result: 11
1 < 2, 回傳 result: 11

我們也可以換個想法，利用找尋目前 n 第一個 1 的前方有幾個 0 (__builtin_clz(n))，再用 31 去減掉前方 0 的個數即可

參考 log2.h，了解到 linux 核心內部 ilog2 的實作方式，其中有以下幾點:

• 會先判斷 n 是否使用了常數時間來編譯
• 不為常數時間則判斷其大小是否小於 4 byte，來執行後續相對應的 fls
• fls 之所以需要減 1 是因為它是以 1 為底

測驗 4

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=\alpha (Y_t+(1-\alpha )Y_{t-1}+(1-\alpha {)}^2{Y}_{t-2}+\dots +(1-\alpha {)}^k{Y}_{t-k}+\dots +Y_0)$

$$S_{t+1}=((S_t\ast 2^{avg->weight}-S_t)+Y_t\ast 2^{avg->factor})\ast 2^{-(avg->weight)}$$
化簡後可得

$$\{\begin{array}{l}{S}_{t+1}=(1-2^{-(avg->weight)})S_t+Y_t\ast 2^{avg->factor-(avg->weight)}\\ \alpha =2^{-(avg->weight)}\end{array}$$
而 avg->factor 的用途在於，不是每個系統都支援浮點運算，因此藉由固定其大小來確保精確度

測驗 5

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-- 用意在於，若是 x 為 2 的冪，最後的 +1 會造成結果錯誤，所以先將 x-- 來避免這個情況發生。
而利用後面一連串的判斷，來快速找出 x 的對數。
最後的回傳會判斷 x 是否還有值，且合併 r 與 shift 的值，並加上最後的進位 1。

處理 x = 0 的狀況

可以發現到，當我們直接帶入 0 及 1 時，程式會回傳結果 1，而這並不是我們所預期的，於是可以如下修改

-   x--;
+   x -= likely(x); /* x -= !!x */

-   return (r | shift | x > 1) + 1;
+   return (r | shift | x > 1) + !(x == 0);

想法是:

• 先判斷是否為 0 ，為 0 則不做 -1
• 最後回傳時再判斷 x 是否為 0，是則不做 +1

如此，即可以避免當 x 為 0 和1，回傳的錯誤結果 1，並維持 branchless

---

quiz 4

題目

測驗 1

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

計算方式採取每 4 個位元為單位，

$v=(b_{31}{b}_{30}{b}_{29}{b}_{28})...b_0$
右移 1 且只把最高位元(以 4 位元為單位)給變成0，以防其影響計算結果，

$n=(0b_{31}{b}_{30}{b}_{29})(0b_{27}...b_1)$

$v-n=v-\lfloor {\displaystyle \frac{v}{2}}\rfloor$
跑完這段即可得，

$v-\lfloor {\displaystyle \frac{v}{2}}\rfloor -\lfloor {\displaystyle \frac{v}{4}}\rfloor -\lfloor {\displaystyle \frac{v}{8}}\rfloor$
接下來會執行 v = (v + (v >> 4)) & 0x0F0F0F0F
令

$v=B_7{B}_6{B}_5{B}_4{B}_3{B}_2{B}_1{B}_0$, (

$B_n=b_{4n+3}{b}_{4n+2}{b}_{4n+1}{b}_{4n+0}$)
所以

$v=0(B_7+B_6)0(B_5+B_4)0(B_3+B_2)0(B_1+B_0)$
最後將 v *= 0x01010101, v >> 24，即可得

$popcount=B_7+B_6+B_5+B_4+B_3+B_2+B_2+B_0$

撰寫出更高效的程式碼 477. Total Hamming Distance

int totalHammingDistance(vector<int>& nums) {
    int res = 0, n = nums.size();
    /* TLE: O(n^2)
     *  for(int i=0;i<nums.size();i++){
     *       for(int j=i+1;j<nums.size();j++){
     *           int xOR = nums[i]^nums[j];
     *           ans += __builtin_popcount(xOR);
     *       }
     *   } 
     */
    for (int i = 0; i < 32; i++) {
        int count = 0;
        for (int j = 0; j < n; j++)
            count += (nums[j] >> i) & 1;
        res += count * (n - count);
    } 
    return res;
}

時間複雜度為

$O(32n)=O(n)$

測驗 2

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

$N=b_{n-1}{b}_{n-2}\cdots b_1{b}_0$，
可以改寫成

$N\equiv \sum _{i=0}^{n-1}{b}_i\cdot (-1{)}^i(mod3)$
我們可以利用 n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)，來達到相同效果
其中還可以利用以下定理來進一步的化簡

$$popcount(x\mathrm{\&}\bar{m})-popcount(x\mathrm{\&}m)=popcount(x\oplus m)-popcount(m)$$
即 popcount(m) = 16, (m = 0xAAAAAAAA)，而為了避免 n = popcount(n ^ 0xAAAAAAAA) - 16 結果為負數，我們可以加入一個足夠大的 3 的倍數的數字給 -16，範例為 39 - 16 -> 23
再來執行 n = popcount(n ^ 0x2A) - 3 來限縮範圍至 -3 ~ 3，最後的 return n + ((n >> 31) & 3) 再將負的轉為正的
其中，參考 mesohandsome 的筆記，發現到要將 n 做轉型才會得到正確結果，即

-   return n + ((n >> 31) & 3);
+   return n + ((int)n >> 31) & 3;

static inline uint32_t is_win(uint32_t player_board)
{
    return (player_board + 0x11111111) & 0x88888888;
}

當連線時，會是 0x7，為了 0x7 + ? == 0x8 要為 true，所以 ? 應為 0x1

static inline int mod7(uint32_t x)
{
    x = (x >> 15) + (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);
}

先將 x 右移 15 來保留其商數，再利用 x & UINT32_C(0x7FFF) 來保留最低的 15 個位元

測驗 3

函式功用

• xt_create: 初始化一棵新的 XTree
• xt_destroy: 利用遞迴的方式來將樹的每個節點的資源給釋放掉
• xt_frist: 找到樹中最小值的節點(最左邊)
• xt_last: 找到樹中最大值的節點(最右邊)
• xt_rotate_left: 執行下圖操作

• xt_rotate_right: 執行下圖操作

• xt_balance: 回傳左子樹樹高減去右子樹樹高
• xt_update: 確認樹是否平衡，不平衡則執行對應的旋轉操作
• xt_insert: 將新節點插入進樹中
• xt_remove: 將對應的節點從樹中移除

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

想法是若欲刪除之節點在右子樹，則用右子樹的最小來填補，即 least ，而填補完需要確認樹是否平衡，所以執行 xt_update(root, xt_right(least))，
反之，則用左子樹的最大，即 most，且更新 xt_update(root, xt_left(most))
最後在判斷整棵樹有沒有平衡，即 xt_update(root, parent)