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}^n{a}_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}or0$。從

$m=n$ 一路試到

$m=0$ ，試驗

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

$a_m$

$$\{\begin{array}{l}{P}_m=P_{m+1}+2^m,if{P}_m^2\le N^2\\ P_m=P_{m+1},if{a}_m=0\end{array}$$

然而如果每次都試驗

$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_0{P}_1+P_1^2$$
並可推展出一般式

$$P_m^2=a_m^2+2a_m{P}_{m+1}+P_{m+1}^2$$
經過移項

$$P_m^2-P_{m+1}^2=2a_m{P}_{m+1}+a_m^2$$
假設

$Y_m=P_m^2-P_{m+1}^2=2a_m{P}_{m+1}+a_m^2$

以及假設

$X_m=N^2-P_m^2$ ，

$X_m\ge 0$ 即為每次判斷

$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\ge 0$ 的問題轉換為

$(X_{m+1}-Y_m)\ge 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{array}{l}{c}_m+d_m ,if a_m=2^m\\ 0 ,if a_m=0\end{array}$

以此推出下一輪

• $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_m{2}^m=\{\begin{array}{l}{c}_m/2+d_m ,if a_m=2^m\\ c_m/2 ,if a_m=0\end{array}$
• $d_{m-1}=(2^{m-1}{)}^2={\displaystyle \frac{d_m}{4}}$

這裡說明初始與中止條件，因為是判斷

$(X_{m+1}-Y_m)\ge 0$ ，所以我們需要準備兩者的初始條件，在程式進行時，會由

$a_n$ 往下試

• $X_{n+1}=N$
• $n$ 是最大的位數，使的
  $d_n^2=a_n^2=(2^n{)}^2=4^n\le 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_0{2}^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}={\displaystyle \frac{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{array}{l}{c}_m/2+d_m ,if a_m=2^m\\ c_m/2 ,if a_m=0\end{array}$ ，不管

$a_m=2^m$ 是否成立

$c_m$ 都需要除以二，這裡以 z >>= 1 實作

且若 則減去

$Y_m$ 得

$X_{m-1}$。

前面提到我們將原本

$X_m\ge 0$ 的問題轉換為

$(X_{m+1}-Y_m)\ge 0$ 同時也等價

$X_{m+1}\ge 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}$ 可以表示如下 :

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

透過上面的轉換，我們將問題轉換為確保

$x$ 的範圍，來達到

$\frac{n}{x}$ 在精度以內，下面是對

$x$ 範圍的探討

$x$ 的右不等式

$x\le \frac{n}{a.b}=\frac{ab}{a.b}=10$ ，因此一定在精確度內。

$x$ 的左不等式

$\frac{n}{a.b9}\le x$ 無法直接證明，透過代入

$l$ 來證明，

$l$的不等式如下:

$c.d\le \frac{l}{x}\le c.d9$ ，透過將

$\frac{n}{a.b9}\le x$ 代入可表示為

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

左不等式：

$$c.d\le cd\times \frac{a.b9}{ab}=cd\times (\frac{a.b}{ab}+\frac{0.09}{ab})=c.d+\frac{0.09cd}{ab}$$
最後推得

$c.d\le c.d+\frac{0.09cd}{ab}$ 左式顯然成立。

右不等式：

$$cd\times \frac{a.b9}{ab}\le c.d9\Rightarrow c.d+\frac{0.09cd}{ab}\le c.d9\Rightarrow c.d+\frac{0.09cd}{ab}\le c.d+0.09$$
其中

$\frac{0.09cd}{ab}\le 0.09$ ，因為

$ab>cd$ 所以右不等式也成立

於是前述猜想也成立，接下來只需要針對

$n$ 來做討論即可，而因為 int tmp = (b[i] - '0') + (a[i] - '0') + carry 所以 tmp 最大是 9+9+1=19，

$n$ 也即為 19 ，下面針對 19 做討論即可

$$\begin{gathered} 1.9\le \frac{19}{x}\le 1.99 \\ \Rightarrow 9.55\le x\le 10 \end{gathered}$$

上述不等式可得知除數介於

$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\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 >> 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) 等於是

$\mathrm{商}\times 8+\mathrm{商}\times 2=\mathrm{商}\times (8+2)=\mathrm{商}\times 10$ ，根據餘數定理

$\mathrm{餘}\mathrm{數}=\mathrm{被}\mathrm{除}\mathrm{數}(in)-\mathrm{商}(div)\times 10\mathrm{除}\mathrm{數}(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{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

首先看到 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{array}{l}{Y}_0,t=0\\ \alpha Y_t+(1-\alpha )\cdot S_{t-1},t>0\end{array}$$

最後看到 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\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 減 1，這樣可以確保當 x 是 2 的冪時的正確性，因為

$x\le 2^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(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=\sum _{i=0}^{n-1}{b}_i\cdot 2^i$ ，由以上的推論可得

$n\equiv \sum _{i=0}^{n-1}{b}_i\cdot (-1{)}^i(mod3)$ 將奇偶數分開

$n\equiv \sum _{i=even}{b}_i-\sum _{i=odd}{b}_i(mod3)$

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

$$popcount(x\mathrm{\&}\bar{m})-popcount(x\mathrm{\&}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(mod3)$ 簡化求

$n(mod3)$
首先假設 n 前的 16 位元是

$c$ ，後 16 位元是

$d$ ，那麼我們就可以將 n 的值寫成

$n=c\cdot 2^{16}+d$ ，再透過

$4^k\equiv 1(mod3)$ 可以將

$n(mod3)$表示為

$$n(mod3)=c\cdot 2^{16}+d(mod3)=c+d(mod3)$$

上述的數學式可轉換成 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(mod3)$ 簡化求

$n(mod3)$ ，同理利用

$8^k\equiv 2^{3k}\equiv 1(mod7)$ 簡化求

$n(mod7)$ ，這裡 (x >> 15) 使用 15 是因為

$2^{3k}$ ，所以是每隔 3 位元有同餘的結果，而我們要將原 32 位元的數縮小，找 16 附近的 3 的倍數 15 。
在這個方法中不像 HackerDelight 提到的例子，繼續縮小輸入的值後查表求餘數，而是透過

$n(mod7)\equiv \lfloor \frac{8}{7}n\rfloor (mod8)$ 求 7 的餘數，先將數字乘以

$\lceil \frac{2^{32}}{7}\rceil$ 也就是 0x2492492 再右移 29 得到剩餘的三位數即為解。

不太明白為何

$n(mod7)\equiv \lfloor \frac{8}{7}n\rfloor (mod8)$

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

測驗 3