2024q1 Homework4 (quiz3+4)

contributed by < MathewSu-001 >

第一周測驗1

程式碼理解

若要求得平方根 N ，以 2的冪相加的表示式為

N^{2} = (a_{n} + a_{n - 1} + a_{n - 2} + . . . . + a_{0})^{2}

根據測驗一以及 Digit-by-digit calculation 所提及的定義，

P_{m} = \sum_{i = 0}^{n - m} a_{n - i}

是迄今為止找到的近似平方根，則

P_{0} = a_{n} + a_{n - 1} + . . . + a_{0}

即為所求平方根

N

。

所以我們就是要去計算出所有

a_{n - i}

的總和，從

a_{n}

開始加到

a_{0}

，而求得

a_{m}

只要試驗

P_{m}^{2} \leq N^{2}

是否成立

{\begin{cases} P_{m} = P_{m + 1} + 2^{m}, & if P_{m}^{2} \leq N^{2} \\ P_{m} = P_{m + 1}, & otherwise \end{cases}

考量到運算成本問題，一種想法是從上一輪的差值

X_{m + 1}

減去

Y_{m}

得到

X_{m}

，

$X_{m} = N^{2} - P_{m}^{2} = (N^{2} - P_{m + 1}^{2}) - (P_{m}^{2} - P_{m + 1}^{2}) = X_{m + 1} - Y_{m}$
$Y_{m} = P_{m}^{2} - P_{m + 1}^{2} = (P_{m + 1} + a_{m})^{2} - P_{m + 1}^{2} = 2 P_{m + 1} a_{m} + a_{m}^{2}$

當

X_{n} = 0

時,找到精確的平方根；如果不是，則

P_{m}

給出平方根的合適近似值，其中

X_{n}

是近似誤差。

進一步調整，將

Y_{m}

拆成

c_{m}, d_{m}

$c_{m} = 2 P_{m + 1} a_{m} = P_{m + 1} 2^{m + 1} (if a_{m} = 2^{m})$
$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}, d_{m}

推出下一輪

c_{m - 1}, d_{m - 1}

$c_{m - 1} = P_{m} 2^{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} = a_{m - 1}^{2} = (2^{m - 1})^{2} = \frac{d_{m}}{4}$

所以可以知道當

c_{- 1}

時，得到值

P_{0}

就是平方根 N 所求。

結合上述方法撰寫演算法來尋求

P_{0} = a_{n} + a_{n - 1} + . . . + a_{0}

，初始化條件

$X_{n + 1} = N^{2}$
$c_{n} = 2 P_{n + 1} a_{n}$ 且
$a_{n}$ 已是已知最大，所以
$P_{n + 1} = 0 \to c_{n} = 0$
$d_{n} = a_{n}^{2} = (2^{n})^{2} = 2^{2 n}$ ，所以可以設
$d_{n} = m$ 得到以下程式碼:

int m = 1UL << ((31 - __builtin_clz(x)) & ~1UL)

上述程式碼可以得到 x 最近的偶數符合

2^{2 n} = 4^{n}

。

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;

在以上程式碼中，

z = c_{n}, m = d_{n}, b = Y_{n}, x = X_{n}

。所以依照上述的演算法，不斷的迴圈直到

m = d_{n} = 0

的時候，那所求得的 z 即為

c_{- 1} = P_{0}

。

取代 `__builtin_clz`

使用 linux 提供的 __fls 來進行修改。因為他回傳的 return 值是找到 x 最高位數，所以引用方法如下:

int z = 0;
-for (int m = 1UL << ( (31 - __builtin_clz(x)) & ~1UL); m; m >>= 2) {
+for (int m = 1UL << ( __fls(x) & ~1UL); m; m >>= 2) {
    int b = z + m;
    z >>= 1;
    if (x >= b)
        x -= b, z += m;               
}

58557c9

第一周測驗2

程式碼解析

第一周測驗3

程式碼解析

從版本一開始看起，

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

可以知道找尋

l o g_{2}

的做法就是取得最高位元的次方是多少。
所以在版本二中，

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

首先就是先看 i 的值有沒有超過

2^{16}

(32 位元的一半)，假如有的話將 result 加上 16，並且數值右移 16 位，然後再依次砍半到

2^{8} 、 2^{4} 、 2^{1}

看數值有無超過，這相當於找 fls。

所以版本三中，

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

因為 __builtin_clz的功用就是回傳 v 中從最高有效位開始的前導零位數，且如果 v 為 0，則結果未定義。所以加上 | 1 就是定義如果 v 為零時，回傳 0。

從測驗一中，我學到 fls 的用途就是找到 v 最高位數，與 31 - __builtin_clz(v | 1) 作用相同，因此我有進行嘗試並且證實最後結果相同。

a1de600

第一周測驗4

程式碼解析

透過Exponentially Weighted Moving Average(EWMA; 指數加權移動平均)的數學定義:

S_{t} = {\begin{cases} Y_{0} & t = 0 \\ α Y_{t} + (1 - α) \cdot S_{t - 1} & t > 0 \end{cases}

可以知道程式碼的 ewma_add 的計算方式

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

在變數上，val 為新的一筆資料；avg 中的 internal 為最後取的平均值，對應到公式中的

S_{t}

; weight 為權重值，對應到公式中的

α

；最後 factor 與

Y_{t}

的計算有關係，為縮放因子。

當 avg->internal 為零時，

S_{t} = Y_{0}, t = 0

當 avg->internal 不為零時，

S_{t} = α Y_{t} + (1 - α) \cdot S_{t - 1} = [Y_{t} + (\frac{1}{α} - 1) \cdot S_{t - 1}] \cdot α

這邊需要注意的是，

α

介在 0 ～ 1 之間，但這邊的 weight 為大於零的 2 的冪，所以在作法上，原本的除法要變成乘法( bitewise operation 做右移)，乘法變成除法( bitewise operation 做左移)。

第一周測驗5

程式碼解析

在一開始我先將程式做改寫，想要知道 x-- 的用途:

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;   
+   return (r | shift | x > 1) ; 
}

就發現如果做出這樣的改動的話，程式碼功能就會跟測驗三的 ilog2 相同。因此先做 x-- ，然後 return 值時再加一的用途，就是利用 ilog2 會以 2 的冪做為一個區段，原本

2^{k} \sim 2^{k + 1} - 1

return 的值為 k ，但改動後會變成

2^{k} + 1 \sim 2^{k + 1}

return 的值為 k+1，達到 ceil 效果。

處理 `x = 0` 的狀況

當 x = 0 時，做 x-- 會取得 x = 0xFFFFFFFF，而導致最後回傳的值是不正確的，所以在做法上參考了 vax-r 同學的做法，只有在 x > 0 才進行減法。

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;

-    x--;
+    x -= !!x;

20f56ef

2024q1 Homework4 (quiz3+4)

第一周測驗1

程式碼理解

取代 __builtin_clz

第一周測驗2

程式碼解析

第一周測驗3

程式碼解析

第一周測驗4

程式碼解析

第一周測驗5

程式碼解析

處理 x = 0 的狀況

第二周測驗1

程式碼解析

Read more

2024q1 Homework5 (assessment)

2024q1 Homework6 (integration)

homework2

2024q1 Homework1 (lab0)

取代 `__builtin_clz`

處理 `x = 0` 的狀況