2024q1 Homework4 (quiz3+4)

contributed by < gawei1206 >

第三周測驗題

測驗 1

以 Digit-by-digit calculation 方式求得

$x$ 的平方根
　　　　

$x$ =

$N^2$ =

$(a_n+a_{n-1}+a_{n-2}+...+a_0{)}^2$,

$a_m$ =

$2^n$

$or$

$0$

$P_m$ =

$a_n+a_{n-1}+a_{n-2}+...+a_{m+1}+a_m$，而

$P_0$ 即為

$x$ 的平方根

$N$

再來可以觀察到：
　　　　

$P_m$ =

$P_{m+1}$ +

$a_m$

因為現在的目的是要求出

$P_0$，也就是要從

$a_n$ 累加到

$a_0$，在每回合判斷是否加入

$a_m$

$$\{\begin{array}{ll}{P}_m=P_{m+1}+2^m,& \text{if }{P}_m^2\le N^2\\ P_m=P_{m+1},& \text{otherwise}\end{array}$$

同先前版本實作的效果

if ((result + a) * (result + a) <= N)

由於每輪計算

$N^2-P_m^2$ 的成本過高，因此改為利用上輪的差值

$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}=2P_{m+1}{a}_m+a_m^2$

也就是紀錄上一輪的

$P_{m+1}$ 來計算

$Y_m$。

再進一步將

$Y_m$ 拆成兩個變數相加，

$2P_{m+1}{a}_m+a_m^2=P_{m+1}{2}^{m+1}+(2^m{)}^2$

$$c_m=P_{m+1}{2}^{m+1},d_m=(2^m{)}^2,Y_m=\{\begin{array}{ll}{c}_m+d_m& if{a}_m=2^m\\ 0& 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-1}=2^{2m-2}={\displaystyle \frac{d_m}{4}}$

最後以上述的推導求出

$P_0$ 即從

$a_n$ 累加到

$a_0$ 的和，以下為參數的對應：
z :

$c_m$
m :

$d_m$
b :

$Y_m$
x :

$X_{m+1}$

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

其中 if (x >= b) 就是以

$X_{m+1}$ 和

$Y_m$ 判斷

$X_m$，即

$N^2-P_m^2$

延伸問題

以 ffs 取代 __builtin_clz:
先前的方法是以 (31 - __builtin_clz(x)) 來計算 MSB，ffs()的回傳是目前 LSB 的位置，以每次求得的 LSB 位置累加並右移，也可以計算出 MSB。

但不知道不依賴 GNU extension 會有什麼優點

int tmp = x, msb = 0;
while (tmp) {
    int s = ffs(tmp);
    msb += s;
    tmp >>= s;
}
msb--;
for (int m = 1UL << (msb & ~1UL); m; m >>= 2)

提供分支和無分支 (branchless) 的實作:
參考 LindaTing0106 同學的作法，將迴圈中的判別式改寫，達到無分支的實作。

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

測驗 2

採用 bitwise operation 來實作除法，但是除數無法直接用 2 的冪來表示的話結果會是不精確的，因此我們的目標是計算新的除數，使 tmp / 10 直到小數點後第一位都是精確的。

題目的除數為 10，再來因為a[i], b[i] 的範圍是 0 ~ 9，carry 則是 0 ~ 1，，所以 tmp 的範圍是 0 ~ 19。

int tmp = (b[i] - '0') + (a[i] - '0') + carry;
carry = tmp / 10;
tmp = tmp - carry * 10;

$x=\frac{2^N}{a}$ 為這邊想要求得的新除數，可以用 2 的冪來表示且

$\frac{tmp}{x}$ 仍然會在精確度內，其中

$N$ 為非負整數

$a$ 為正整數。

證明猜想成立後，接下來只需要針對 19 來做討論即可。

$a.b\le \frac{n}{x}\le a.b9\Rightarrow 1.9\le \frac{19}{x}\le 1.99\Rightarrow 9.55\le x\le 10$

接下來需要找出可以符合條件的

$N$ 和

$a$，題目以

$2^N=128，a=13$ 為其中一解，使得

$\frac{128}{13}\approx 9.84$ 符合條件。

嘗試透過 bitwise operation 來配對數值，首先要得到

$\frac{13tmp}{8}=\frac{tmp}{8}+\frac{tmp}{2}+tmp$

(tmp >> 3) + (tmp >> 1) + tmp

接下來再左移 3 bit 後就可以得到

$13tmp$，但在右移時會有部分位元被捨棄，因此要另行處理。

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

但在這部分我有些疑惑，因為上述程式碼所計算出來的結果有時可能不是

$13tmp$，而只是接近

$13tmp$

雖然最後要再右移 7 bit 結果應該是不影響，但如果真的要將被捨棄的位元加回去，應該如下較為正確：

q = ((((x >> 3) + (x >> 1) + x) << 3) + (d0 << 2) + d2) >> 7;
r = tmp - (((q << 2) + q) << 1);

再透過

$r=f-g\cdot Q$ 算出餘數。

最後包裝成題目的程式碼如下

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 = (in | 1) - (in >> 2) 算出

$x=\frac{3}{4}\cdot in$，其中 (in | 1) 表示當 in 為偶數時會將它加一，再來 q = (x >> 4) + x 的意思是

$\frac{1}{16}\cdot \frac{3}{4}\cdot in+\frac{3}{4}\cdot in=\frac{51}{64}\cdot in\approx 0.79\cdot in$，後續四行重複的操作是為了將結果逼近

$0.8\cdot in$，最後右移 3 bit 得到

$\frac{1}{10}\cdot in$。
一樣透過

$r=f-g\cdot Q$ 算出餘數，q & ~0x7 與 *div << 3 是一樣意思的，再加上 *div << 1 求得

$10\cdot div$，即可計算出餘數。

測驗 3

實作 ilog2 的目標是找出二進位表示的 msb，在第二種實作方法中的想法是以二元搜尋快速找到 msb，相較於方法一逐位元的判斷來說更有效率。

運用 GNU extension __builtin_clz 來改寫，為了避免 0 的例外狀況需將輸入與 1 做 or。

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

測驗 4

在對 ewma 結構體初始化時，會確保 weight 及 factor 的值是 2 的冪，目的是為了保證後續位元操作的正確性，代替乘除法的使用。

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$ 的值，在等號右邊的 internal 代表的是

$S_{t-1}$，而 val 代表的是

$Y_t$，由下方的程式碼可以發現，在 internal 是 0 的時候程式會執行 (val << avg->factor)，所以這邊推測上方的 FFFF 也是做相同的操作，再來觀察
(((avg->internal << EEEE) - avg->internal) + (val << avg->factor)) >> avg->weight
這段程式碼與

$\alpha Y_t+(1-\alpha )\cdot S_{t-1}$ 的關係，可以轉換成

$((S_{t-1}\cdot 2^E-S_{t-1})+Y_t)\times 2^{-weight}$，可以發現到

$2^{-weight}$ 就是

$\alpha$ ，再將

$2^{-weight}$ 乘進去想要算出

$(1-\alpha )\cdot S_{t-1}$，所以 EEEE 就是
avg->weight。

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

測驗 5

與第三題類似，目的是要算出 MSB，但這邊要求出

$\lceil lo{g}_2(x)\rceil$，所以當 x 為二的冪時需要先將 x-- ，使得最後結果正確。

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

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

後續步驟與第三題找出 MSB 的部分相同，不過是以 branchless 的方式實作，以變數 r 紀錄位移量，因為位移量都是二的冪，所以 r |= shift 可以看做 r += shift，
最後 return (r | shift | x > 1) + 1;可以看做兩部分，r | shift 是將前一次的 shift 加上，而在最後一次位移後 x 的值可能會是 0x3 或是 0x2，所以後半部會判斷 x 是否大於 1 並將它加上，最後再加上 1 達到取 ceil 的效果。

延伸問題

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

我們需要判斷 x 是否大於 0 再將其減一，在之前的作業中有看過 likely 還有 unlikely，其中這兩個巨集的實作中有用到 !!(x) 的用法，目的是透過兩次 NOT 來確保值一定是 0 或 1，所以可將程式中 x-- 的部分修改如下：

+ x -= !!(x);
- x--;