2024q1 Homework4 (quiz3+4)

contributed by < padaray >

quiz3

測驗 1

計算平方根

版本一：
對輸入的參數 N 以 2 為底數取 log，意義是當我們將 N 取 log2 時，可以得到 most significant bit (msb)。
對 1 左移 msb 個位元後存入變數 a ，意義是找出一個必定大於 N 的平方根的數，利用此數來逼近找到 N 的平方根。
最後使用 while 迴圈，計算 a + result 的平方是否大於 N，若小於 N 則將此數存入 result 變數，對 a 右移 1 bit，繼續計算，直到 a = 0。

版本二：
版本二的差別是不使用 math.h 的 log2 函式，直接對參數 N 進行右移，右移至參數為 0 時，則代表找到 most significant bit。此方法將程式的可移植性增強，並且不影響原本的函式原型 (prototype) 和精度 (precision)

- int msb = (int) log2(N);
+ int msb = 0;
+ int n = N;
+ while (n > 1) {
+     n >>= 1;
+     msb++;
+ }

版本三：

完整程式碼(無遮蔽處)

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

__builtin_clz(x) ：是 gcc 內建處理二進位的其中一個函式巨集，此函式會回傳 counting leading zero，幫助我們找到 most significant bit。

在 for 迴圈中左移 31 - __builtin_clz(x) bit 數得到 msb，& ~1UL 可以確保最後一個 bit 為 0，理由是變數 m 必須是偶數對應到

$d_n=(2^n{)}^2$。而 m 右移兩個 bit 是對應到

$d_{m-1}=d_m/4$

測驗 2

針對正整數在相鄰敘述進行 mod 10 和 div 10 操作，如何減少運算成本？

除法原理：

$f=g\times Q+r$

• $f$: 被除數
• $g$: 除數
• $Q$: 商
• $r$: 餘數

程式碼如下：

carry = tmp / 10;  // carry 代表商(Q)
tmp = tmp - carry * 10;  // tmp 代表餘數(r)

如果我們想要使用 bitwise operation 來實作除法，如果被除數包含了除了 2 以外的因數，就會造成結果不精確，所以使用以下方式。
假設除數為 10，我們必須找到

${\displaystyle \frac{1}{10}}\approx {\displaystyle \frac{a}{2^N}}$，最後發現是

$\approx {\displaystyle \frac{13}{128}}$
程式碼如下：

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

(((tmp >> 3) + (tmp >> 1) + tmp) << 3) 使用 bitwise operation 取得一個分子有 13 的分數，因此用

${\displaystyle \frac{tmp}{8}}+{\displaystyle \frac{tmp}{2}}+tmp={\displaystyle \frac{13tmp}{8}}$，之後再左移三位消掉分母 8。
d0、d1、d2 是為了將 tmp 右移時損失的 bit 數加回，加回後再右移 7 位也就是除以 128 得到商數。
(((q << 2) + q) << 1) 相當於

$(q\times 4+q)\times 2$，也就是最一開始程式碼的 carry * 10。

最終程式碼可包裝為以下：

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

測驗 3

考慮 ilog2 可計算以 2 為底的對數，且其輸入和輸出皆為整數。

版本一：
初始化變數 log 為負一，直接對參數 i 進行右移，每次右移將 log++，右移至參數 i 為 0 時

版本二：
完整程式碼(無遮蔽處)

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

為了避免數字過大每次位移一個 bit 太慢，例如：65536 需要位移 16 個 bit，因此加入條件判斷，若大於特定 2 的冪次方，直接移動該次方的位元數。

版本三：

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

測驗 4

解釋 EWMA 的實作

Exponentially Weighted Moving Average (EWMA) 是一種統計方法，透過對每筆資料給予不同的權重計算平均，通常越新的資料權重就越高，越舊的資料權重越小，好處是可以適應新的資料變化，又保留一定的歷史資料影響。

數學定義：

$S_t=\{\begin{array}{ll}{Y}_0,& t=0\\ \alpha Y_t+(1-\alpha )\cdot S_{t-1},& t>0\end{array}$

$\alpha$ : 舊資料加權降低的程度，越高代表降低越快。介於 0 ~ 1 之間

$Y_t$ : 時間 t 時的資料

$S_t$ : 時間 t 時算出的 EWMA

以下為程式碼實作：

 * ewma_add() - Exponentially weighted moving average (EWMA)
 * @avg: Average structure
 * @val: Current value
 *
 * Add a sample to the average.
 */
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;
}

internal 代表目前計算的 ewma (數學式中的

$S_t$) ， (avg->internal << avg->weight) 代表數學式中的

$\alpha ={\displaystyle \frac{1}{2^{weight}}}$ (但不知道原因為何)，因此

$1-\alpha ={\displaystyle \frac{2^{weight}-1}{2^{weight}}}$。
當 internal 為 0 時，internal 代入 (val << avg->factor) 寫成公式為

$2^{factor}\cdot val$。
當 internal 不為 0 時，internal 代入 (((avg->internal << avg->weight) - avg->internal) + (val << avg->factor)) >> avg->weight 寫成公式為

${\displaystyle \frac{2^{weight}-1}{2^{weight}}}\cdot internal+{\displaystyle \frac{2^{factor}}{2^{weight}}}\cdot val$，轉換一下可以得到

$(1-\alpha )\cdot internal+2^{factor}\cdot \alpha \cdot val$

測驗 5

延伸測驗 3 使用下方程式碼計算

$\lceil {\log }_{2}x\rceil$

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 是否大於 0xFFFF，意義是確定 x 的 msb 是否位於高位的 16 bit，若為真，將 x 右移 16 bit。
第七、八、九行判斷參數 x 在這剩下的 16 bit 中，msb 是否位於高位的 8 bit，若為真，將 x 右移 8 bit，將右移的位數紀錄在 r 中。
第十行到第十四行依此類推。
第四、十五行將 x 先 -1 後 + 1 回傳是為了處理 2 的冪次方，若不處理會讓 2 的冪次方多計算一次，例：

$2^7$ 經過無條件進位法，會變成 8，因此要先減 1

quiz4

測驗 1

popcount(population count)