Try   HackMD

2024 Homework4 (quiz3+4)

第三周測驗

2024q1 第 3 週測驗題

測驗 1 : Using bitwise operations to compute square root

先看原程式碼

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

編譯方式

$ gcc -o i_sqrt i_sqrt.c -lm

因為呼叫 log2,所以需要連結 libm (math library),即 -lm 選項

假設 N =

3610 =
1001002
msb = 5 => a =
1000002

接著進入迴圈中

迴圈次數 a result 是否進入迴圈
1
1000002
0 y
2
100002
0 y
3
10002
0 y
4
1002
1002
y
5
102
1102
y
6
12
0 n

在此 bitwise 操作中, (result + a) 的平方逐一比較 N 值,當值較小時就疊加當前的值,即 result += a ,直到為 N 的開方根值。

此題一開始用十進位去看比較容易明白, a >>= 1 看作 a / 2
編譯方式:$ gcc -o i_sqrt i_sqrt.c -lm

版本二使用右移操作直到 n 值為 0 ,藉此計算 msb

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

版本三使用 Digit-by-digit calculation 方法進行開方根。將平方數拆成 2 的冪相加

x=N2=(an+an1+...+a0)2 ,其中 N 即為所求的開方根

Pm=an+an1+...+am
P0=an+an1+...+a0=N
即為所求

推論出:

Pm=Pm+1+am

後續將求出所有

am=2m or
0
,從
m=n
一直到
m=0
,而求得
am
只須比較
Pm2N2
,有兩種情況:
{Pm=Pm+1+2m  ; Pm2N2Pm=Pm+1            ; otherwise

由於每輪逐一比較的運算成本過高,提出了另個想法
N2Pm2
的值表示位
Xm

{Xm=N2Pm2(1)Xm+1=N2Pm+12(2)
將(1)式減(2)式可得
XmXm+1=Pm2+Pm+12(3)

其中

Pm 平方為

Pm2=Pm+12+2amPm+1+am2=Pm+12+(2Pm+1+am)am

Ym=Pm2Pm+12=(2Pm+1+am)am

最後將

Ym 代回到(3)式可得:

Xm=Xm+1Ym

Ym 拆解
cm,dm

cm=Pm+12m+1
dm=(2m)2

得:

{Ym=cm+dm ;if am=2mYm=0              ;if am=0
藉由位元運算從
cm,dm
推出下一輪

cm1=Pm2m=(Pm+1+am)2m=Pm+12m+am2m
{cm/2+dm  if am=2mcm/2            if am=0

dm1=dm/4

Linux 核心原始程式碼 linux/block/blk-wbt.c

static void rwb_arm_timer(struct rq_wb *rwb)
{
	struct rq_depth *rqd = &rwb->rq_depth;

	if (rqd->scale_step > 0) {
		/*
		 * We should speed this up, using some variant of a fast
		 * integer inverse square root calculation. Since we only do
		 * this for every window expiration, it's not a huge deal,
		 * though.
		 */
		rwb->cur_win_nsec = div_u64(rwb->win_nsec << 4,
					int_sqrt((rqd->scale_step + 1) << 8));
	} else {
		/*
		 * For step < 0, we don't want to increase/decrease the
		 * window size.
		 */
		rwb->cur_win_nsec = rwb->win_nsec;
	}

	blk_stat_activate_nsecs(rwb->cb, rwb->cur_win_nsec);
}

尚未完成

測驗 2 : Using bitwise operations to perform mod 10 and div 10

參考《Hacker's Delight》,採用 bitwise operation 來實作除法

n a2N 此題選用 a=13 N=7 =>
n 1327
,其中
1327
0.10156

完整程式碼

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×1327 改寫成
tmp×13÷16÷8
,其中
tmp×138×8=13tmp
可從下方程式碼得到

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

由於 tmp 右移三位會遺失原本最低三位的 bit ,因此將分別用 d0d1d2 來儲存原先的 bit ,再將遺失的 bit 補回來,最後在除以128,即為

13128tmp,得到一個近似
tmp÷10
的結果 ( q )。

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

上述程式碼中不太理解為何補足遺失的最低三位 bit 是加上 + d0 + d1 + d2 ,在我實作中發現就算缺少這段程式碼也不會影響其結果,其大致的觀念我能理解是為了增加其精度,但又為何不直接加上 + d2 就好。

下方程式碼根據餘數定理:

fg×Q=r,可以得到
r
(餘數)

  • f
    : 被除式
  • g
    : 除式
  • Q
    : 商
  • r
    : 餘數
r = tmp - (((q << 2) + q) << 1);

其中(((q << 2) + q) << 1) 等於

(q×4+q)×2=q×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));   
}

uint32_t x = (in | 1) - (in >> 2) 得到

34in
uint32_t q = (x >> 4) + x; 得到
q=34in÷16+34in=5164in
,其中
5164in
約為
0.796875in
近似
810

再經由 bitwise 操作使 0.796875 更接近 0.8

    q = (q >> 8) + x; // q = 0.799987793
    q = (q >> 8) + x; // q = 0.799999952
    q = (q >> 8) + x; // q = 0.799999998
    q = (q >> 8) + x; // q = 0.8

最後 *div = (q >> 3); 即為

div=810in÷8=110in
((q & ~0x7) + (*div << 1));
810in+110in×2=in10×10
,其中
in10
為商。

在由餘式定理:

fg×Q=r 求得 *mod 值為何

練習撰寫不依賴任何除法指令的 % 9 (modulo 9) 和 % 5 (modulo 5) 程式碼

尚未進行

測驗 3 : Using bitwise operations to calculate the logarithm base 2

ilog2 計算以 2 為底的對數,且其輸入和輸出皆為整數。

稍微改寫一下原程式碼方便識讀,舉例

i=410=1002

while ( i >= 2 ) 原 i 後 i log
第 1 次迴圈 100 10 1
第 2 次迴圈 10 1 2

輸出 log = 2

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

改寫後的程式碼,將 i 用一定數量的位元去進行判斷,有效提升數字較大時的計算速度,其原因為原程式碼每個位元逐一去判斷,反覆進行多次的迴圈,改後後判斷分成 4 個階段,若 i 大於

216282421則直接進行位元操作。

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

在 Linux 核心原始程式碼找出 log2 的相關程式碼 (或類似的形式),並解說應用案例

看到 log2 為底的形式就想到 clzffsfls 等相關的操作都有類似的手法。

linux/arch/arc/include/asm/bitops.h 中找到 fls 的實作,其作用為返回最高有效 bit 位(由最低位往左數),當判斷式為 0 進行位移操作,其手法相似 log2 ,不同於位元操作向左進行。

該順序不同於 log2constant_fls 從 1 起始,當沒有有效位返回 0 。

static inline int constant_fls(unsigned int x)
{
	int r = 32;

	if (!x)
		return 0;
	if (!(x & 0xffff0000u)) {
		x <<= 16;
		r -= 16;
	}
	if (!(x & 0xff000000u)) {
		x <<= 8;
		r -= 8;
	}
	if (!(x & 0xf0000000u)) {
		x <<= 4;
		r -= 4;
	}
	if (!(x & 0xc0000000u)) {
		x <<= 2;
		r -= 2;
	}
	if (!(x & 0x80000000u))
		r -= 1;
	return r;
}

測驗 4 : Using bitwise to calculate EWMA

參考 stevendd5432024q1 第 3 週測驗題

Exponentially Weighted Moving Average
賦予資料指數衰減的權重,主要目的使越舊的資料權重越低,反之越新的資料對整體平均值的影響越大,以下為此統計手法的公式:

St={Y0t=0αYt+(1α)St1t>0

  • α
    表示歷史資料加權降低的程度,介在 0 ~ 1 之間,越高的
    α
    會使歷史資料減少的越快
  • Yt
    表示時間 t 時的資料點
  • St
    表示時間 t 時的 EWMA
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;
}

@internal : 看作

t 當前的位置
@factor : 縮放因子,必須是 2 的冪次
@weight : 資料的權重,在這
α=2avg>weight
透過右移來達到
α=2(avg>weight)=12avg>weight
,必須是 2 的冪次

struct ewma *ewma_add(struct ewma *avg, unsigned long val)
{
    avg->internal = avg->internal
                        ? (((avg->internal << weight) - avg->internal) +
                           (val << avg->factor)) >> avg->weight
                        : (val << avg->factor);
    return avg;
}

一開始

t=0 ,因此
St=Y0
對應到 avg->internal = (val << avg->factor);
,接著 avg->internal = (((avg->internal << avg->weight) - avg->internal) + (val << avg->factor)) >> avg->weight 的理解較為困難,看起來與題意的公式不同,在 stevendd543 筆記中的舉例有非常清楚解釋:

αα1=1,將公式改寫成

α1St=α1(αYt+(1α)St1)
          =Yt+(α11)St1

其中

(α11)St1 等效於 ((avg->internal << avg->weight) - avg->internal)
Yt
等效於 (val << avg->factor) ,最後將整個公式 >> avg->weight 即為
α1St×α=St
(當時
t
的EWMA)。

unsigned long ewma_read(const struct ewma *avg)
{
    return avg->internal >> avg->factor;
}

將計算的平均值除以縮放因子而得到真正的平均值。

在 Linux 核心原始程式碼找出 EWMA 的相關程式碼 (或類似的形式),至少該涵蓋無線網路裝置驅動程式

測驗 5 : Test 3 Extension

延伸測驗 3,已知輸入必為大於 0 的數值 (即 x > 0),以下程式碼可計算 [log2] ,也就是 ceil 和 log2 的組合並轉型為整數。

x-- 目的為避免 x 為 2 的冪,接著使用類似測驗 3 的作法,在指定的位元寬度中比較,若條件成立就會進行位元操作,下方程式碼寫的更加複雜因此實際帶入數值去進一步解釋。

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

case 1
x = 0000 0000 0000 0000 0000 1000 0000 0000

x-- = 0000 0000 0000 0000 0000 0111 1111 1111
r = (x > 0xFFFF) << 4 = 0 << 4 = 0
x = x >> r = 0000 0000 0000 0000 0000 0111 1111 1111
shift = (x > 0xFF) << 3 = 1 << 3 = 8
x >>= shift = 0000 0000 0000 0000 0000 0000 0000 0111
r |= shift = 8
shift = (x > 0xF) << 2 = 0 << 2 = 0
x >>= shift = 0000 0000 0000 0000 0000 0000 0000 0111
r |= shift = 8
shift = (x > 0x3) << 1 = 1 << 1 = 2
x >>= shift = 0000 0000 0000 0000 0000 0000 0000 0001
return (r | shift | x > 1) + 1 = return (11)
  • 透過 shift = (x > 0xFF) << 3 決定 x 要右移的大小
  • r |= shift 等效 r += shift ,用來記錄由第 0 bit 往左數到 msb 的有效位
  • 最後結果把先前 x-- 的值加回來

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

int ceil_ilog2(uint32_t x)
{
    uint32_t r, shift;
+   bool mark = x;
+   x -= !!x 
-   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) + mark;
}

x = 0ceil_ilog2 應要回傳 0 ,因此設立一個 mark 來紀錄最後是否 +1 。

Linux 核心原始程式碼
尚未完成

第四周測驗

2024q1 第 4 週測驗題

測驗 1 : Total Hamming Distance

閱讀 population count 了解計算二進位表示式中有多少個位元為 1 、 了解 Hamming distance 所代表的意思。

更進階的操作參照測驗1,大致內容為每4個位元為一個單位,求出其中的 set bits 個數,32位元中總共有 8 組,將每組算得的 set bits 個數加總在一起,裡面的推導過程需詳讀就能得知最後為何 v >> 24 即為 popcount 。

考題針對 LeetCode 477. 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;
}

Leetcode 範例描述到

Input: nums = [4,14,2]
The answer will be:
HammingDistance(4, 14) + HammingDistance(4, 2) + HammingDistance(14, 2) = 2 + 2 + 2 = 6.

在考題程式碼中為兩兩成對的去計算 popcount ,因此會出現HammingDistance(4, 14)與HammingDistance(14, 4)這種情況,因此最後回傳要除以 2 (即 return total >> 1)。

不依賴 popcount,嘗試以上述歸納的規則,針對 LeetCode 477. Total Hamming Distance 撰寫出更高效的程式碼

參考 solution 中的解答

int totalHammingDistance(int* nums, int numsSize){
    int ans = 0;
    for (int i = 0; i < 32; ++i) {
        int c = 0;
        for (int j = 0; j < numsSize; ++j) {
            c += (nums[j] >> i) & 1;
        }
        ans += c * (numsSize - c);
    }
    return ans;
}

使用 bitwise 逐一比較二進制中每個位元是否為 1 ,其判斷式 (nums[j] >> i) & 1
每當比較完一位元,便計算其 HammingDistance 的值 ( ans += c * (numsSize - c) )

n 'th bit 3 2 1 0
input 4 0 1 0 0
input 14 1 1 1 0
input 2 0 0 1 0
  1. 第 0 個 bit 觀察:
    全部都是 0 => Hamming Distance = 0
  2. 第 1 個 bit 觀察:
    1 有兩個 => Hamming Distance = (2)* (numsSize - 2)
  3. 第 3 個 bit 觀察:
    1 有一個 => Hamming Distance = (1)* (numsSize - 1)

測驗2 : Remainder by Summing digits

了解 [ 離散數學 ] 同餘(Mod)的性質及定義
詳讀 Hacker's Delight 中提及 Remainder by Summing digits

以除 3 為例,

11(mod3) ,
21(mod3)

2k
{11(mod3),if k is even21(mod3),if k is odd

n 以二進制表示成
n=bn1bn2...b1b0

n=i=0n1bi2i ni=0n1bi(1)i(mod3)

將上式結合 popcount 寫成程式

n = popcount(n & 0x55555555) - popcount(n & 0xAAAAAAAA)

0x55 = 0b01010101
0xAA = 0b10101010

再由定理

popcount(x&m)popcount(x&m)=popcount(xm)popcount(m)

定理證明在 Hacker's Delight 裡有提到,頁碼:P13

化簡成

n = popcount(n ^ 0xAAAAAAAA) - 16
int mod3(unsigned n)
{
    n = popcount(n ^ 0xAAAAAAAA) + 23; //n = popcount(n ^ 0xAAAAAAAA) - 16 + 39;
    n = popcount(n ^ 0x2A) - 3;
    return n + ((n >> 31) & 3);        //避免 n 為負數值
}

為了避免 n 變成負數,因此加上一個足夠大的 3 的倍數,範例中為加 39 ,在第二輪應用使 n 介於 -3~2 之間而非 -3~3 之間。
return n + ((n >> 31) & 3)主要用作避免 n 為負數值,若為 n 為正則不受影響。

舉例 若 n = -1
-1 -> 0b1111 1111 1111 1111 1111 1111 1111 1111  # n
0b11                                             # (n >> 31) & 3
0b10                                             # n + ((n >> 31) & 3)

lookup table

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

使用查表的方式去尋找對應的餘數。
tictactoe.c

將上述手法應用於第三次作業中,追求更有效的賽局判定機制。

測驗 3 : XTree

AVL Tree 是一種Binary Search Tree (BST) 實作方法,透過高度計算及 balance factor 判斷是否需要旋轉來平衡。

red-black tree 遵循以下規則來達成平衡:

  1. 節點只有黑或紅
  2. 根節點必為黑
  3. 若節點是紅色的,則其子節點必為黑色,反之不必為真(黑色節點的下個節點可為黑或紅)
  4. 所有從該節點走到其子樹下的末端節點 (leaf) 的上之黑色節點數量必定

此題 XTree 的平衡機制與 AVL tree 相似,利用 hint 來評估是否需要做 rotate,然而 xt_node 的 hint 的計算方式是,其左右節點 hint 最大值 +1 並且只有在被呼叫到 xt_update 時才會被更新。

#define xt_root(r) (r->root)
#define xt_left(n) (n->left)
#define xt_right(n) (n->right)
#define xt_rparent(n) (xt_right(n)->parent)
#define xt_lparent(n) (xt_left(n)->parent)
#define xt_parent(n) (n->parent)

xt_create : 初始化 XTree
xt_destroy

利用遞迴

    if (xt_root(tree))
        __xt_destroy(tree, xt_root(tree));

將樹裡的每個節點釋放,在 __xt_destroy()

    if (xt_left(n))
        __xt_destroy(tree, xt_left(n));

    if (xt_right(n))
        __xt_destroy(tree, xt_right(n));

xt_left (左子樹)或 xt_right (右子樹)存在則繼續往下走直到盡頭 tree->destroy_node(n) 來釋放。

xt_first

用遞迴的方式尋找最左邊的節點

    if (!xt_left(n))
        return n;

    return xt_first(xt_left(n));

xt_last

xt_first 一樣,但變成尋找最右變的節點

xt_rotate_left

    struct xt_node *l = xt_left(n), *p = xt_parent(n);

    xt_parent(l) = xt_parent(n);
    xt_left(n) = xt_right(l);
    xt_parent(n) = l;
    xt_right(l) = n;

這段程式碼為下圖操作







Tree

before


p

p



n

n



p->n





l

l



n->l





r

r



n->r





A

A



l->A





B

B



l->B











Tree

after


p

p



l

l



p->l





A

A



l->A





n

n



l->n





B

B



n->B





r

r



n->r





    if (p && xt_left(p) == n)
        xt_left(p) = l;
    else if (p)
        xt_right(p) = l;

假設 p 的右子樹不存在







Tree

before


null



null1



p

p



p->null





l

l



p->l





l->null1





n

n



l->n











Tree

after


null



null1



p

p



p->null





l

l



p->l





l->null1





n

n



l->n





    if (xt_left(n))
        xt_lparent(n) = n;

重新定義 n 節點的 parent

xt_rotate_right

操作手法與 xt_rotate_left 相似,因此僅用圖示







Tree

before


p

p



n

n



p->n





l

l



n->l





r

r



n->r





A

A



r->A





B

B



r->B











Tree

after


p

p



r

r



p->r





n

n



r->n





B

B



r->B





l

l



n->l





A

A



n->A





xt_balance

與 AVL Tree 中 Balance Factor of T, BF(T) 的平衡判斷相似,這裡的 hint 為左子樹高 - 右子樹高

xt_max_hint

選擇左子樹與右子樹的 hint 的最大值 +1 即為節點 n 的 hint 值(更新 hint)

xt_update

若樹不平衡 xt_balance < -1 進行 xt_rotate_right , xt_balance > 1 進行 xt_rotate_left

xt_insert

插入新節點在樹中

xt_remove

當移除一節點 del 時,將該右子樹的最左來填補,即 least 來替換掉原本的 del ,並檢查是否平衡 xt_update,反之則用左子樹的最右節點來點補,若被刪除的節點為 leaf 將不用進行交換。

static void __xt_remove(struct xt_node **root, struct xt_node *del)
{
    if (xt_right(del)) {
        struct xt_node *least = xt_first(xt_right(del));
        if (del == *root)
            *root = least;

        xt_replace_right(del, least);
        xt_update(root, xt_right(least));
        return;
    }

    if (xt_left(del)) {
        struct xt_node *most = xt_last(xt_left(del));
        if (del == *root)
            *root = most;

        xt_replace_left(del, most);
        xt_update(root, xt_left(most));
        return;
    }

    if (del == *root) {
        *root = 0;
        return;
    }

    /* empty node */
    struct xt_node *parent = xt_parent(del);

    if (xt_left(parent) == del)
        xt_left(parent) = 0;
    else
        xt_right(parent) = 0;

    xt_update(root, parent);
}