2024q1 Homework4 (quiz3+4)

contributed by < jimmylu890303>

第三週測驗題

題目

測驗一

版本一

#include <math.h>
int i_sqrt(int N)
{
    int msb = (int) log2(N);
    int a = 1 << msb;
    int result = 0;
    while (a != 0) {
        if ((result + a) * (result + a) <= N)
            result += a;
        a >>= 1;
    }
    return result;
}

在版本一中，會透過 log2(N) 去尋找 N 的最高位元，再從高位逐步去尋找滿足 (result + a) * (result + a) <= N 的位元。

以輸入 N = 36 為例， msb = 5 且 a = 100000

a	result	(result + a) * (result + a) <= N
100000(32)	0	X
10000(16)	0	X
1000(8)	0	x
100(4)	0	O
10(2)	4	O
1(1)	6	X

所以 result 的值為 6 ，在此版本中使用到 log2 ，所以需要在編譯時多新增編譯選項 -lm

版本二

在版本二中，去改善版本一使其不依賴 log2 並保持原本的函式原型 (prototype) 和精度 (precision)，而改善的作法就是使用 shift 的方式去找出最高位的位元。

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

版本三

在版本三中使用 __builtin_clz(x) 的巨集來找出最高位在哪個位元算出 leading zero 個數。

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	m	z	b	z(after shift)	x >= b	x-=b	z+=m
36	10000(16)	0	16	0	O	20	16
20	100(4)	16	20	8	O	0	12
0	1(1)	12	13	6	X	0	6

測驗二

針對正整數在相鄰敘述進行 mod 10 和 div 10 操作，如何減少運算成本？
採用 bitwise operation 來實作除法

先考慮除以 10 的情況，除以 10 相當於乘以

${\displaystyle \frac{1}{10}}$
但由於

${\displaystyle \frac{1}{10}}$ 無法以二進位的方式去表示，所以只能透過一個近似的值去替代
而

${\displaystyle \frac{1}{10}}\approx {\displaystyle \frac{13}{128}}$ ，因為

${\displaystyle \frac{128}{13}}\approx 9.84$
所以 x 除以 10 相當於

$x\times {\displaystyle \frac{13}{128}}$
更進一步拆解上述的式子

$x\times 13÷8÷16$

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

考慮以下式子，相當於

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

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

由於往右作 shift 3個位元會造成最右邊三個位元資訊損失，所以在往左 shift 後要將最右邊三個位元加回來。所以以下操作值會等於

$13\times tmp$

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

所以最後往右 shift 7個位元相當於除上128，這樣就會得到

$tmp÷10$ 的結果

q = ((((tmp >> 3) + (tmp >> 1) + tmp) << 3) + d0 + d1 + d2) >> 7;

而 tmp 取 10 的餘數則利用餘數定理

$r=f-g\cdot Q$

r = tmp - (((q << 2) + q) << 1);

q 為商數，則以下操作為

$q\times 10$，即可算出餘數 r

(((q << 2) + q) << 1)

以上的概念可以包裝成以下程式碼

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

在以下操作中，

$x={\displaystyle \frac{3}{4}}\times in$

uint32_t x = (in | 1) - (in >> 2);

$q={\displaystyle \frac{3}{4}}\times in+{\displaystyle \frac{3}{64}}\times in={\displaystyle \frac{51}{64}}\times in\approx 0.79\times in$

uint32_t q = (x >> 4) + x;

而這邊的操作是將 0.79 逼近至 0.8

x = q;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;
q = (q >> 8) + x;

$div={\displaystyle \frac{8}{10}}\times in\times {\displaystyle \frac{1}{8}}={\displaystyle \frac{1}{10}}\times in$

*div = (q >> 3);

考慮以下操作
q >> 3 為商數
而 q & ~0x7 將最低三位清零，相當於商數往左 shift 3位，即為商數的 8 倍
*div << 1 為商數的兩倍
((q & ~0x7) + (*div << 1)) 所以這裡相當於商數乘上 10

*mod = in - ((q & ~0x7) + (*div << 1));   

測驗三

版本一

在以下程式碼中，使用右移操作來計算 log 值

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

以 8 為例

round	log	i	i >>= 1	log++
1	-1	1000	100	0
2	0	100	10	1
3	1	10	1	2
4	2	1	0	3

以 7 為例

round	log	i	i >>= 1	log++
1	-1	0111	011	0
2	0	011	01	1
3	1	01	0	2

由以上觀察可發現 shift 操作次數是由 msb 在第幾個位元而定，所以這個實作也可以藉由 counting leading zero 來完成

版本二

以下程式碼是將版本一改善，假設輸入資料為 65536 ，在版本一中的 while 迴圈次數需要 65536+1 次，而在版本二中只需要 1 次的迴圈，對於輸入資料很大時，可以減少所需的時間。

- while (i) {
-    i >>= 1;
-    log++;
- }
+ 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;
+    }

版本三

使用 GNU extension __builtin_clz 來計算最高位前方有幾個 0 ，在透過 31 - 0 個個數即為最高位元的位置，此值即為 log 的值。

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

測驗四

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

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

以下為 ewma 的結構， internal 為

$S_{t-1}$ ，factor 為縮放因子， weight 為權重

/* Exponentially weighted moving average (EWMA) */
struct ewma {
    unsigned long internal;
    unsigned long factor;
    unsigned long weight;
};

考慮以下實作，

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

$\alpha ={\displaystyle \frac{1}{2^{weight}}}$

$1-\alpha ={\displaystyle \frac{2^{weight}-1}{2^{weight}}}$
若 avg->internal 為 0 ，則為

$2^{factor}\cdot val$
若 avg->internal 不為 0 ，則為

${\displaystyle \frac{2^{weight}-1}{2^{weight}}}\cdot internal-{\displaystyle \frac{2^{factor}}{2^{weight}}}\cdot val$

測驗五

以下程式碼計算

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

以下用於判別 msb 在左邊 16 位元或右 16 位元，若在左 16 位元則 x > 0xFFFF 條件成立，則將 x shift 16位元。

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

以下用於判別 msb 在這 16 位元中的左 8 位元或右 8 位元，r 用於紀錄 shift 幾個 bits

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

以此類推即可求出

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

改進程式碼

commit - 8d405ff

第一個為 x=0 時，做 x– 會導致 x 值變成 0xffffffff，所以作了以下改進
x -= (!!x) 會使 x > 0 時才會作減法的動作。
第二個為 x=1 時，應該要回傳 0 的值，但是在原本的程式中會回傳 1 ，所以去檢查 !(r | shift | x > 0) 有無成立，若該條件成立代表 x = 0 或者 x = 1 的狀態，再將回傳值清零。

-    x --;
+    x -= (!!x);
-    return (r | shift | x > 1)+1;
+    x = !(r | shift | x > 0) ? 0 : (r | shift | x > 1)+1;
+    return x;       

ceil_log2 in Linux

linux/tools/testing/selftests/mm/thuge-gen.c

int ilog2(unsigned long v)
{
	int l = 0;
	while ((1UL << l) < v)
		l++;
	return l;
}

第四週測驗題

題目

測驗一

popcount

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

在最初的 popcount實作中，是使用迴圈將 LSB 的值清 0 ，
這邊使用特別操作的手法為 v &= (v-1) ， 此外 v &= (v-1) 的操作也可以用於判斷是否為 2 的冪次方。

0001 0100  # 20          ; LSB in bit position 2
0001 0011  # 20 - 1
0001 0000  # 20 & (20 - 1)

n = -(~n) 為 n++ 的操作，

$-N=\sim N+1$ 為二的補數操作，則

$-(\sim N)=N+1$

branchless popcount

unsigned popcount_branchless(unsigned v)
{
    unsigned n;
    n = (v >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;
    n = (n >> 1) & 0x77777777;
    v -= n;

    v = (v + (v >> 4)) & 0x0F0F0F0F;
    v *= 0x01010101;                                     

    return v >> 24;
}

根據以下公式，我們可以計算 popcount(x)

$popcount(x)=x-\lfloor {\displaystyle \frac{x}{2}}\rfloor -\lfloor {\displaystyle \frac{x}{4}}\rfloor -\lfloor {\displaystyle \frac{x}{8}}\rfloor -...-\lfloor {\displaystyle \frac{x}{2^{31}}}\rfloor =\sum _{n=0}^{31}(2^n-\sum _{i=0}^{n-1}{2}^i)b_n=\sum _{n=0}^{31}{b}_n$

以下操作為計算

$x-\lfloor {\displaystyle \frac{x}{2}}\rfloor -\lfloor {\displaystyle \frac{x}{4}}\rfloor -\lfloor {\displaystyle \frac{x}{8}}\rfloor$
0x77777777 的設定是因為這邊是以 4 個 bits 為單位 (nibble) ，為了防止另外一組的 nibble 值 shift 過來。

n = (v >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;
n = (n >> 1) & 0x77777777;
v -= n;

所以操作完後的值會結果如下(只考慮最右邊的 nibble)，這個 nibble 即代表這 4 個 bits 的 popcount 之值 :

$(2^3-2^2-2^1-2^0)b_3+(2^2-2^1-2^0)b_2+(2^1-2^0)b_1+2^0{b}_0=b_3+b_1+b_1+b_0$

 3       2          1             0
b3 (b2-b3) (b1-b2-b3) (b0-b1-b2-b3)

假設

$B_n$ 為 nibble ，所以以下操作可以看出它的結果

v = (v + (v >> 4)) & 0x0F0F0F0F;

B7 B6 B5 B4 B3 B2 B1 B0  // v
 0 B7 B6 B5 B4 B3 B2 B1  // (v >> 4)
// (v + (v >> 4))
B7 (B7+B6) (B6+B5) (B5+B4) (B4+B3) (B3+B2) (B2+B1) (B1+B0)
// (v + (v >> 4)) & 0x0F0F0F0F
0 (B7+B6) 0 (B5+B4) 0 (B3+B2) 0 (B1+B0)

最後在透過 v *= 0x01010101 將各個 nibble 相加

A6 = B7+B6 A4 = B5+B4 A2 = B3+B2 A0 = B1+B0

                         0 A6  0 A4  0 A2  0 A0
                      x  0  1  0  1  0  1  0  1
---------------------------------------------------
                         0 A6  0 A4  0 A2  0 A0
                   0 A6  0 A4  0 A2  0 A0  0
             0 A6  0 A4  0 A2  0 A0  0
       0 A6  0 A4  0 A2  0 A0  0
---------------------------------------------------
                            ↑_______________________A6+A4+A2+A0

可以發現 A6+A4+A2+A0 是我們期望的結果，所以透過 v >> 24 取得。

利用 popcount 計算 Total Hamming distance

int totalHammingDistance(int* nums, int numsSize)
{
    int total = 0;;
    for (int i = 0;i < numsSize;i++)
        for (int j = 0; j < numsSize;j++)
            total += __builtin_popcount(nums[i] ^ nums[j]); 
    return total >> 1;
}

使用 XOR 運算 可以得出兩數字的差異，在利用巨集 __builtin_popcount 即可以算出兩者的 Hamming Distance ，而 Total Hamming distance 就是每個數字都跟其他數字相比，所以時間複雜度會被限制在

$N^2$ 。

改善程式碼

int totalHammingDistance(int* nums, int numsSize)
{
    int total = 0;
    for(int i=0 ;i<32;i++){
        int sum=0;
        for(int n=0;n<numsSize;n++){
            int lsb = (nums[n] & 0x1);
            sum += lsb;
            nums[n]>>=1;
        }
        total += sum * (numsSize-sum);
    }
    return total;
}

假設輸入資料為 4, 14, 4 三筆資料，我們從最低位元到最高位元將每個數字的該位元加總，可以發現以下規律 : 該位元的加總 * (數字個數 - 該位元的加總) = 該位元的 Total Hamming distance ，並且時間複雜度為

$N$

 numsize = 3
---------------------------------------------
 4 = 0100
14 = 1110
 4 = 0100
 
 0th bit : 0 + 0 + 0 = 0 // 0 * (numsize - 0) = 0
 1th bit : 0 + 1 + 0 = 1 // 1 * (numsize - 1) = 2
 2th bit : 1 + 1 + 1 = 3 // 3 * (numsize - 3) = 0
 3th bit : 0 + 1 + 0 = 1 // 1 * (numsize - 1) = 2
 
 total = 0 + 2 + 0 + 2 = 4

測驗二