2024q1 Homework4 (quiz3+4)

contributed by < ssheep773 >

第 3 週測驗題

2024q1 第 3 週測驗題

測驗 `1`

目的: 計算 N 的整數平方根

方法一 :

通過逐步逼近來計算 N 的平方根。它首先找到 N 的最大位元位置 msb，並將最大位元位置對應的數值儲存在 a \((a = 2^{msb})\) 。然後逐步減小 a 的值，同時檢查 \((result + a)^2\) 是否小於或等於 N。如果是的話，就將 result 加上 a，以逼近 N 的平方根。這樣的操作會持續進行，直到 a 的值為 0。

int i_sqrt(int N)
{
    int msb = 0;
    int n = N;
    while (n > 1) {
        n >>= 1;
        msb++;
    }
    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}

方法二 :
利用 Digit-by-digit calculation 提及常見的直式開方法 \((x + y)^2 = x^2 + 2xy + y^2\) 的變形實作開平方根。

若要求 \(x\) 的平方根，可以假設 \(x = N^2\) ，那麼 \(N\) 即為我們要求的平方根數值，接著，將 \(N^2\) 拆成

\[ N^2=(a_{n}+a_{n-1}+a_{n-2}+...+{a_0})^2 = [(a_{n}+a_{n-1}+...+{a_1})+({a_0})]^2 \]

經過整理

\[ N^2 = (\sum_{i=0}^na_i)^2 = (\sum_{i=1}^{n} a_i)^2 + 2a_0(\sum_{i=1}^{n}a_i) + a_0^2 \]

這裡假設 \(P_m = \sum_{i=m}^n a_i\) ，那麼 \(P_0 = a_{n}+a_{n-1}+a_{n-2}+...+{a_0}\) 即為我們要求的 \(N\)， \(P_0\) 可在分解為 \(P_0 = (a_{n}+a_{n-1}+...+{a_1})+({a_0}) =P_1 + a_0\)
可以得到一般式
\[ P_m = P_{m+1} +a_m \]
我們的目的是要試出所有 \(a_m\) ，\(a_m = 2^m or 0\)。從 \(m = n\) 一路試到 \(m = 0\) ，試驗 \(P_m^2 \le N^2\) 是否成立來求得 \(a_m\)

\[ \begin{cases} P_m = P_{m+1} +2^m \ \ ,if \ \ P_m^2 \le N^2\\ P_m = P_{m+1} \ \ ,if \ \ a_m=0\\ \end{cases} \]

然而如果每次都試驗 \(P_m^2 \le N^2\) 的運算成本會太高，一種想法是透過上一輪\(N^2 - P_{m+1}^2\) 的差值推估 \(N^2 - P_m\) 相關的證明如下:

因為前面 \(P_m\) 的假設 \(N\) 可以整理為
\[ N^2 = P_0^2 = a_0^2 + 2a_0P_1 + P_1^2 \]
並可推展出一般式
\[ P_m^2 = a_m^2 + 2a_mP_{m+1} + P_{m+1}^2 \]
經過移項
\[ P_m^2 -P_{m+1}^2 = 2a_mP_{m+1} + a_m^2 \]
假設 \(Y_m = P_m^2 -P_{m+1}^2 = 2a_mP_{m+1} + a_m^2\)

以及假設 \(X_m = N^2 -P_m^2\) ，\(X_m\geq0\) 即為每次判斷 \(a^m\) 的條件，可以推廣 \(X_{m+1} = N^2 -P_{m+1}^2\) 為前一次的差值
\[ Y_m = P_m^2 -P_{m+1}^2 = (N^2 -X_m) - (N^2 -X_{m+1}) = X_{m+1}-X_m \]
我們就可以得到，從上一輪的差值 \(X_{m+1}\) 減去 \(Y_m\) 得到 \(X_m\)
\(X_m = N^2 -P_m^2 = X_{m+1} - Y_m\)

將原本 \(X_m \geq 0\) 的問題轉換為 \((X_{m+1} - Y_m) \geq 0\) 即使用到上一輪的差值來計算。
進一步調整，將 \(Y_m\) 拆成 \(c_m\) 和 \(d_m\)

\(c_m = 2P_{m+1}a_{m} = P_{m+1}2^{m+1}\)
\(d_m = a_{m}^2 = (2^{m})^2\)
\(Y_m = \begin{cases} c_m + d_m 　　, if 　a_m = 2^m \\ 0　　　　　　, if 　a_m = 0 \\ \end{cases}\)

以此推出下一輪

\(c_{m-1} = 2P_m a_{m-1} = P_m 2^m = (P_{m+1} + a_m)2^m = P_{m+1}2^m + a_m2^m = \begin{cases} c_m/2 +d_m 　, if 　a_m = 2^m\\ c_m/2　　　　, if 　a_m = 0 \\ \end{cases}\)
\(d_{m-1} = (2^{m-1})^2 = \dfrac{d_m}{4}\)

這裡說明初始與中止條件，因為是判斷 \((X_{m+1} - Y_m) \geq 0\) ，所以我們需要準備兩者的初始條件，在程式進行時，會由 \(a_n\) 往下試

\(X_{n+1} = N\)
\(n\) 是最大的位數，使的 \(d_n^2 = a_n^2 = (2^n)^2 = 4^n \leq N^2\) 所以用迴圈逼近 \(d_n\)
\(c_n = 0\) 因為 \(c_n = 2P_{n+1}a_{n}\) 而 \(a_n\) 是最高位數所以 \(P_{n+1} = 0\)

因為 \(d_m\) 恆大於 \(0\) 所以終止條件可以設為 d != 0 ，又 \(c_{-1} = P_02^0 = P_0 = a_n + a_{n-1} +...+ a_0\) 故求出 \(c_{-1}\) 則等於求出 \(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 或 1 的情況

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

變數 m 對應到 \(d_m\) ，其中因為 m = \(d_m = a_{m}^2 = (2^{m})^2\) ，所以 m 不可能為奇數，因此要 & ~1UL ，同時再前一個步驟也先處理 x = 1 的情況
因為 \(d_{m-1} = \dfrac{d_m}{4}\) 所以每次迴圈 m 右移 2 。

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

x 對應到 \(X_m\) 呼應初始條件 \(X_{n+1} = N\)，變數 b 對應到 \(Y_m\) ，z 對應到 \(c_m\) 同時也說明 \(Y_m = c_m + d_m\)

而計算 \(c_{m-1} = \begin{cases} c_m/2 +d_m 　, if 　a_m = 2^m\\ c_m/2　　　　, if 　a_m = 0 \\ \end{cases}\) ，不管 \(a_m = 2^m\) 是否成立 \(c_m\) 都需要除以二，這裡以 z >>= 1 實作

且若則減去 \(Y_m\) 得 \(X_{m-1}\)。

前面提到我們將原本 \(X_m \geq 0\) 的問題轉換為 \((X_{m+1} - Y_m) \geq 0\) 同時也等價 \(X_{m+1} \geq Y_m\) ，若成立則 \(X_{m-1} = X_m - Y_m\)

int b = z + m;
z >>= 1;
if (x >= b)
    x -= b, z += m;

ffs / fls 取代 __builtin_clz

ffs 是找最低為 1 的位元，所以透過迴圈的方式由低往高位元找，到最高的有效位元的位置為止，並且還需要將其減一才會與 (31 - __builtin_clz(x) 相等，若沒有減一，雖對於結果沒有影響，但是會多花費一輪的運算成本。

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

+       int tmp = x, 
+       int msb = 0;
+       while (tmp) {
+           int i = ffs(tmp);
+           msb += i;
+           tmp >>= i;
+       }
        
    int z = 0;
+   for (int m = 1UL << ((msb-1) & ~1UL); 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;
}

測驗 `2`

針對正整數在相鄰敘述進行 mod 10 和 div 10 操作，減少運算成本
原始的方法:

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

可利用除法原理將 mod 10 和 div 10 合併。根據除法原理：
\(f = g \times Q +r\)
其中 \(f:\) 被除數， \(g:\) 除數， \(Q:\) 商， \(r:\) 餘數
將程式碼改寫為

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

這邊參考《Hacker's Delight》，採用 bitwise operation 來實作 10 的除法，但是因為 10 有 5 這個因數，無法完全用 2 的冪項來表示，因此結果會是不精確的。

我們的目標是，得到 tmp / 10 且直到小數點後第一位都是精確的。

為此，提出一個猜想，假設我目標的最大數是 \(n\) ，而 \(l\) 是小於 \(n\) 的非負整數，則只要以 \(n\) 算出來的除數在 \(n\) 除以該除數後在精確度內，則 \(l\) 除與該除數仍然會在精確度內，證明如下：
假設 \(n=ab\) ( \(a\) 是十位數 \(b\) 是個位數)，\(l=cd\) ( \(c\) 是十位數 \(d\) 是個位數)。

首先看到 \(\frac{n}{x}\)，因為我們是要確保小數點後第一位是精確的，所以 \(\frac{n}{x}\) 可以表示如下 :

\(a.b\leq\frac{n}{x}\leq a.b9\\ \Rightarrow\frac{n}{a.b9}\leq x\leq\frac{n}{a.b}\)

透過上面的轉換，我們將問題轉換為確保 \(x\) 的範圍，來達到 \(\frac{n}{x}\) 在精度以內，下面是對 \(x\) 範圍的探討

\(x\) 的右不等式 \(x \leq \frac{n}{a.b} = \frac{ab}{a.b} = 10\) ，因此一定在精確度內。

\(x\) 的左不等式 \(\frac{n}{a.b9}\leq x\) 無法直接證明，透過代入 \(l\) 來證明，\(l\)的不等式如下:

\(c.d\leq\frac{l}{x}\leq c.d9\) ，透過將 \(\frac{n}{a.b9}\leq x\) 代入可表示為

\[ c.d\leq l \times \frac{a.b9}{n} \leq c.d9 \Rightarrow c.d\leq cd \times \frac{a.b9}{ab} \leq c.d9 \]

左不等式：
\[c.d\leq cd\times\frac{a.b9}{ab} = cd\times (\frac{a.b}{ab}+ \frac{0.09}{ab}) = c.d + \frac{0.09cd}{ab} \]
最後推得 \(c.d\leq c.d + \frac{0.09cd}{ab}\) 左式顯然成立。

右不等式：
\[cd \times \frac{a.b9}{ab} \leq c.d9 \Rightarrow c.d + \frac{0.09cd}{ab}\leq c.d9 \Rightarrow c.d + \frac{0.09cd}{ab}\leq c.d + 0.09 \]
其中 \(\frac{0.09cd}{ab}\leq 0.09\) ，因為 \(ab > cd\) 所以右不等式也成立

於是前述猜想也成立，接下來只需要針對 \(n\) 來做討論即可，而因為 int tmp = (b[i] - '0') + (a[i] - '0') + carry 所以 tmp 最大是 9+9+1=19，\(n\) 也即為 19 ，下面針對 19 做討論即可

\[ 1.9\leq \frac{19}{x}\leq 1.99\\ \Rightarrow 9.55\leq x\leq10 \]

上述不等式可得知除數介於 \(9.55\) 至 \(10\) 之間皆可程式中達到相同效果，即除數直到小數點後第一位都是精確的，而我們的商會比實際的還要大有正偏差。

我們透過 bitwise operation 找到介於 \(9.55\) 至 \(10\) 之間的除數，除數的算式必定是 \(\frac{2^N}{a}\) ，其中 \(N\) 是非負整數， \(a\) 是正整數，而計算 \(n\) 的商可以寫成 \(\frac{an}{2^N}\)。

因此只需要透過查看 \(2^N\) 在配對適合的 \(a\) 即可
這裡 \(2^N = 128, a= 13 ,\frac{128}{13} \approx 9.84\) 即為可用的除數解

接著看到實際的 bitwise operation

unsigned d2, d1, d0, q, r;

d0 = q & 0b1;
d1 = q & 0b11;
d2 = q & 0b111;
q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7;
r = tmp - (((q << 2) + q) << 1);

q 是要求的商 \(\frac{an}{2^N}\) ，在這裡是 \(\frac{13tmp}{2^7}\)，我們首先要得到 13 ，透過 (tmp >> 3) + (tmp >> 1) + tmp) << 3 得到 \(\frac{13tmp}{8} = \frac{tmp}{8}+\frac{4tmp}{8}+tmp\) 再 << 3 得到 \(13tmp\)，最後再除以 128 即得到我們要的商 \(\frac{13tmp}{128}\)

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

而因為原有的 tmp 在右移後，會有部分位元被捨棄，因此要另行處理。

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

然而我覺得這步驟可以省略，因為後續的操作會在將 tmp >> 7 等於是補回來的位元馬上又被捨棄，而我思考若是在產生 13 而右移的位數大於 \(2^N\) 的 \(N\) 時，則需要處理，而目前分別是 3 跟 7 ，所以可以省略不處理。

在來看到包裝的函式

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

uint32_t x = (in | 1) - (in >> 2) 將 \(x =\frac{3}{4}in\)
再透過 uint32_t q = (x >> 4) + x 等於是 \(q = \frac{3}{4*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 >> 3) 是計算商，而我們前面求的 q 是 \(\frac{8}{10}in\) ，所以需要再除以 8 也就是 q >> 3 。

*mod = in - ((q & ~0x7) + (*div << 1)) 是計算餘數，其中的 q & ~0x7 等於是 *div << 3 ，那麼就程式碼可以轉換成 *mod = in - (*(div << 3) + (*div << 1)) 後半部的 *(div << 3) + (*div << 1) 等於是 \(商\times8 + 商\times2 = 商 \times (8+2) = 商 \times10\) ，根據餘數定理 \(餘數 = 被除數(in) - 商(div) \times10 除數(10)\)

測驗 `3`

ilog2 可計算以 2 為底的對數，且其輸入和輸出皆為整數，這裡注意到因為輸出為整數，所以輸入的最高二進位有效位元位置即為結果。

方法一 :
由最低位元往高位元處找最高二進位有效位元位置

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

方法二 :
透過二元搜尋的方法尋找

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

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

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

在 linux/log2.h 中使用 fls 找出最高有效位元並且說明輸入 0 是未定義的。

/*
 * non-constant log of base 2 calculators
 * - the arch may override these in asm/bitops.h if they can be implemented
 *   more efficiently than using fls() and fls64()
 * - the arch is not required to handle n==0 if implementing the fallback
 */
#ifndef CONFIG_ARCH_HAS_ILOG2_U32
static __always_inline __attribute__((const))
int __ilog2_u32(u32 n)
{
	return fls(n) - 1;
}
#endif

#ifndef CONFIG_ARCH_HAS_ILOG2_U64
static __always_inline __attribute__((const))
int __ilog2_u64(u64 n)
{
	return fls64(n) - 1;
}
#endif

測驗 `4`

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

\[ S_t = \begin{cases} \ Y_0 \qquad\qquad\qquad\qquad ,t=0\\ \ \alpha Y_t +(1-\alpha) \cdot S_{t-1} \ \ ,t>0\\ \end{cases} \]

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

首先看到 ewma 結構，internal 儲存我們的 EWMA 值也就是 \(S_t\) ，值得注意的是它的型態是 unsigned long 是無號整數，所以我們可以知道這裡是使用定點數，並在 ewma_init() 註解處中有說明 factor 是 Scaling factor ，其中註解處有提到internal 可記錄最大值的計算公式為 ULONG_MAX / (factor * weight) ，其中 factor 是因為做 scaling 提高數值精度變相會犧牲可容納的對大值，而 weight 的部分不太理解，只能夠推測因為程式中有做 avg->internal << avg->weight 與 avg->internal >> avg->weight 與 factor 一樣。

並在 ewma_read 可以看到輸出時還原 scaling ，可見 avg->internal >> avg->factor

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

接著看到 ewma_init 函式，目的是初始化 ewma，其中使用 is_power_of_2 來確保 weight 跟 factor 是 2 的冪，因為 internal 是定點數若是與 2 的冪做乘除法，可以用位元運算。
is_power_of_2 是運用當數值 n 是 2 的冪時，n 跟 n-1 在二進位時，不會有相同的位數會是錯開的，利用這個特點檢測

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

\[ S_t = \begin{cases} \ Y_0 \qquad\qquad\qquad\qquad ,t=0\\ \ \alpha Y_t +(1-\alpha) \cdot S_{t-1} \ \ ,t>0\\ \end{cases} \]

最後看到 ewma_add ，這裡可想而知會用到前面提到的 \(S_t\) ，先說明 avg->internal 等於 0 的情況，將 avg->internal 設為輸入的 val 對應到 \(Y_0\) ，但別忘了 scaling 所以還要 val << avg->factor

在來看到大於 0 的情況，第一次看到 (((avg->internal << avg->weight) - avg->internal) +(val << avg->factor)) >> avg->weight 時感覺跟 \(\alpha Y_t +(1-\alpha) \cdot S_{t-1}\) 公式有關，但是仔細看時發現正負號會有錯，這邊再看過 ShchangAnderson 同學的筆記後才了解，其實 \(\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

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

測驗 `5`

這個程式碼實做一個函式，用來計算輸入的 32 位元無號整數 \(x\) 的 2 次方指數值向上進位的結果。也就是說，對於傳入的參數 \(x\) ，回傳最小的整數 \(n\)，滿足 \(x\leq2^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 減 1，這樣可以確保當 x 是 2 的冪時的正確性，因為 \(x\leq2^n\)，當等於時不用進位。

接著看到

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

這裡其實跟測驗 3 是相同的概念，將測驗三的變數調整方便比較，其中 0xFFFF = 2^16 ，並且將數字都以 2 的冪表示。
此方法與測驗三的差異在於紀錄 r 的方式是透過 bitwise 操作而非加法運算， r | shift 等效於 r + shift 。

/*   ilog2   */
while (x >= 2^8) {
    r += 2^3;
    x >>= 2^3;
}

最後 return (r | shift | x > 1) + 1 的部分我們分開看，其中 r | shift 是將前一部分的 r |= shift 合併進來， x > 1 是為了要處理前面 shift = (x > 0x3) << 1 會忽略 x= 0x2 的情況，最後再加上一作為取上界 (ceil) 。

改進程式碼，使其得以處理 x = 0 的狀況，並仍是 branchless

當 x=0 時會因為減一變成 0xFFFFFFFF ，而這不是我們想要的結果。
我想到在 likely / unlikely 中有用到 !!(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;

第 4 週測驗題

2024q1 第 4 週測驗題

測驗 `1`

population count 簡稱 popcount 或叫 sideways sum，計算數值的二進位表示中 1 的個數。

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

v &= (v - 1) : reset LSB 是將最低有效位 (LSB) 設為 0 。這讓我想到這跟前面提到的 is_power_of_2 觀念類似，都是判斷 (n & (n - 1)) == 0 只是 is_power_of_2 只判斷一次，而 popcount_naive 是重複判斷來計算 1 的個數。
n = -(~n) 等同 n++，因為在二補數系統中 \(-n = \sim n+1\) ，經過移項 \(-(\sim n) = n +1\)

針對 LeetCode 477. Total Hamming Distance 撰寫出更高效的程式碼。

透過觀察 Total Hamming Distance 的規則

可以得到一般式:
觀察每個 Number 的第 i 個 bit ，若有 cnt 個 1 ，其餘皆為 0 ，則 Hamming Distance : (cnt) * (numsize - cnt)

我們要走訪所有的 bit ，在每次走訪時計算第 i 個 bit 為 1 的 number 個數(cnt)，其中使用 (nums[j] >> i) & 1 得知目前的number 第 i 個 bit 是否為 1 ，最後透過公式 (cnt) * (numsize - cnt) 計算得到第 i 個 bit 的 Hamming Distance ，最後加總所有位元的 Hamming Distance 即為 nums 中所有整數的 Hamming Distance 。

int totalHammingDistance(int* nums, int numsSize) {
    int total = 0;;
    for (int i = 0;i < 32;i++){
        int cnt = 0;
        for (int j = 0; j < numsSize;j++)
            cnt += ((nums[j] >> i) & 1);
        total += cnt *(numsSize - count);
    }
    return total;
}

測驗 `2`

Remainder by Summing digits
如何不使用任何除法就算出某數除以另一個數的餘數呢？若除數符合 \(2^k \pm 1\) ，則可運用以下手法來達成這個目的。

以除數為 3 為例，我們可以觀察到 \(1 \equiv 1 (mod \ \ 3)\) 且 \(2 \equiv -1 (mod \ \ 3)\)。將後者不斷的乘上前者可得

\[2^k \equiv \begin{cases} 1 (mod \ \ 3), \ \ k \ even\\ -1 (mod \ \ 3), \ \ k \ odd \\ \end{cases}\]

因此若 n 的二進位表示為 \(b_{n-1}b_{n-2}b_{n-3}...b_1b_0\) ，則 \(n = \sum_{i=0}^{n-1} b_i \cdot 2^i\) ，由以上的推論可得 \(n \equiv \sum_{i=0}^{n-1} b_i \cdot (-1)^i \ \ (mod \ \ 3)\) 將奇偶數分開 \(n \equiv \sum_{i=even} b_i - \sum_{i=odd} b_i \ \ (mod \ \ 3)\)

可以利用 population count 來計算位元和，因此可以將上式表示為 n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA) ，其中 & 0x55555555 是留下 n 的奇數位元， & 0xAAAAAAAA 是留下 n 的偶數位元，並利用以下定理可以進一步化簡。

\[ popcount(x \& \overline{m}) - popcount(x \& m) = popcount(x \oplus m) - popcount(m) \]

於是可以寫為 n = popcount(n ^ 0xAAAAAAAA) - 16 ，計算結果的範圍會落在 -16 至 16 之間，但為了避免 n 變成負數，我們要將它加上一個足夠大的 3 的倍數，文中的例子是加上 39 讓數值介於 23 到 55 ，我們再做一次與上面相同的操作，只是這次的數值較小，相差 32 所以跟 0x2A = 0b101010 做運算而 popcount(0x2A) = 3 所以是減三，最後區間會介於 -3 <= n <= 2 ，還需要對負數加 3 ，所以最後 (n >> 31) & 3 判斷最高位元判斷正負來決定是否加 3 。

int mod3(unsigned n)
{
    n = popcount(n ^ 0xAAAAAAAA) + 23;  // Now 23 <= n <= 55.
    n = popcount(n ^ 0x2A) - 3;  // Now -3 <= n <= 2
    return n + ((n >> 31) & 3);
}

另一種變形是利用 lookup table ，其中 n 介於 0 到 32 透過查表直接得知。

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

然而在 HackerDelight 中提到不使用指令的方法

int remu3(unsigned n) {
 n = (n >> 16) + (n & 0xFFFF); // Max 0x1FFFE.
 n = (n >> 8) + (n & 0x00FF); // Max 0x2FD.
 n = (n >> 4) + (n & 0x000F); // Max 0x3D.
 n = (n >> 2) + (n & 0x0003); // Max 0x11.
 n = (n >> 2) + (n & 0x0003); // Max 0x6.
 return (0x0924 >> (n << 1)) & 3;
}

這段程式碼主要是利用 \(4^k \equiv 1 (mod \ \ 3)\) 簡化求 \(n \ (mod \ \ 3)\)
首先假設 n 前的 16 位元是 \(c\) ，後 16 位元是 \(d\) ，那麼我們就可以將 n 的值寫成 \(n=c\cdot2^{16}+d\) ，再透過 \(4^k \equiv 1 (mod \ \ 3)\) 可以將 \(n (mod \ \ 3)\)表示為

\[ n (mod \ \ 3) = c\cdot2^{16}+d(mod \ \ 3) = c+d(mod \ \ 3) \]

上述的數學式可轉換成 n = (n >> 16) + (n & 0xFFFF) ，以此類推重複 5 次將 n 縮小至最大 0x6
最後一行類似查表，將餘數值存在 0x0924 ，透過位元操作取值，其中因為 3 的餘數是 {0,1,2} 可以用 2 個位元即可表示，所以 0x0924 是每隔 2 位元存一個餘數值，並透過 (n << 1) 與 & 3 搭配取得。

0x0924   =  00 00 10 01 00 10 01 00
(n << 1) =     12 10  8  6  4  2  0   
 n       =      6  5  4  3  2  1  0

tictactoe.c 中使用到的 mod7 方法說明
上述提到利用 \(4^k \equiv 2^{2k}\equiv 1 (mod \ \ 3)\) 簡化求 \(n \ (mod \ \ 3)\) ，同理利用 \(8^k \equiv 2^{3k} \equiv 1 (mod \ \ 7)\) 簡化求 \(n \ (mod \ \ 7)\) ，這裡 (x >> 15) 使用 15 是因為 \(2^{3k}\) ，所以是每隔 3 位元有同餘的結果，而我們要將原 32 位元的數縮小，找 16 附近的 3 的倍數 15 。
在這個方法中不像 HackerDelight 提到的例子，繼續縮小輸入的值後查表求餘數，而是透過 \(n (mod \ 7) \equiv \lfloor\frac{8}{7}n\rfloor (mod \ 8)\) 求 7 的餘數，先將數字乘以 \(\lceil\frac{2^{32}}{7}\rceil\) 也就是 0x2492492 再右移 29 得到剩餘的三位數即為解。

不太明白為何 \(n (mod \ 7) \equiv \lfloor\frac{8}{7}n\rfloor (mod \ 8)\)

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

tictactoe.c 是一個井字遊戲，他的特別之處在於儲存棋子的方式，以往可能是用一個矩陣儲存井字中每個位置的棋子狀態，這個程式是用棋子可以連成的線來儲存，這個方法的優點在於可以更快速的判斷贏家，也可以更方便取得連線的資訊，

2024q1 Homework4 (quiz3+4)

第 3 週測驗題

測驗 1

ffs / fls 取代 __builtin_clz

測驗 2

測驗 3

測驗 4

測驗 5

改進程式碼，使其得以處理 x = 0 的狀況，並仍是 branchless

第 4 週測驗題

測驗 1

針對 LeetCode 477. Total Hamming Distance 撰寫出更高效的程式碼。

測驗 2

測驗 3

Read more

嵌入式作業系統分析與實做 Final Project-自動餵食器

2024q1 Homework5 (assessment)

2024q1 Homework2 (quiz1+2)

2024q1 Homework6 (integration)

測驗 `1`

測驗 `2`

測驗 `3`

測驗 `4`

測驗 `5`

測驗 `1`

測驗 `2`

測驗 `3`