# 2025q1 Homework2 (quiz1+2) contributed by < `rota1001` > ## Q1 間接指標、鏈結串列與記憶體配置 ### 測驗 1 ```cpp static inline void list_insert_before(list_t *l, list_item_t *before, list_item_t *item) { list_item_t **p; for(p = &l->head; *p != before; p = &(*p)->next) ; *p = item; item->next = before; } ``` 有一個單向鏈結串列 `l`,我們要將 `item` 這個節點接到 `before` 節點之前,於是使用了一個間接指標 `p`,指向一個指向下一個節點的指標的指標(一開始指向 head),好處是我可以透過 `p` 修改它指向的值。就算 `head` 一開始是 `NULL` 還是可以通過 `p` 在不做特判時修改,以下以圖例解釋: 如果我今天要 `list_insert_before(l, &node2, &node3)`,首先我讓 `p` 指向 `head`。 ```graphviz digraph linked_list { node [shape=record]; rankdir=LR; p [label = "p"] head [label = "head"] node0 [label = "<name>node0|{<n>&node1}"] node1 [label = "<name>node1|{<n>&node2}"] node2 [label = "<name>node2|{<n>NULL}"] p->head head->node0:name node0:n->node1:name node1:n->node2:name } ``` 發現這個時候對 `p` 取值是 `&node0`,於是 `p = &(*p)->next`。 ```graphviz digraph linked_list { node [shape=record]; rankdir=LR; p [label = "p"] head [label = "head"] node0 [label = "<name>node0|{<n>&node1}"] node1 [label = "<name>node1|{<n>&node2}"] node2 [label = "<name>node2|{<n>NULL}"] p->node0:n head->node0:name node0:n->node1:name node1:n->node2:name } ``` 發現這個時候對 `p` 取值是 `&node1`,於是 `p = &(*p)->next`。 ```graphviz digraph linked_list { node [shape=record]; rankdir=LR; p [label = "p"] head [label = "head"] node0 [label = "<name>node0|{<n>&node1}"] node1 [label = "<name>node1|{<n>&node2}"] node2 [label = "<name>node2|{<n>NULL}"] p->node1:n head->node0:name node0:n->node1:name node1:n->node2:name } ``` 這個時候發現 `p` 取值是 `&node2`,於是跳出迴圈,將 `p` 指向的值改成 `&node3`: ```graphviz digraph linked_list { node [shape=record]; rankdir=LR; p [label = "p"] head [label = "head"] node0 [label = "<name>node0|{<n>&node1}"] node1 [label = "<name>node1|{<n>&node3}"] node2 [label = "<name>node2|{<n>NULL}"] node3 [label = "<name>node3|{<n>uninitialized}"] p->node1:n head->node0:name node0:n->node1:name node1:n->node3:name } ``` 然後把 `node3` 的 `next` 改成 `&node2`: ```graphviz digraph linked_list { node [shape=record]; rankdir=LR; p [label = "p"] head [label = "head"] node0 [label = "<name>node0|{<n>&node1}"] node1 [label = "<name>node1|{<n>&node3}"] node2 [label = "<name>node2|{<n>NULL}"] node3 [label = "<name>node3|{<n>&node2}"] p->node1:n head->node0:name node0:n->node1:name node1:n->node3:name node3:n->node2:name } ``` 接下來看一下如果 `l` 原本是空的,進行 `list_insert_before(l, NULL, &node0)` 操作會發生什麼吧,首先初始狀態是這樣,`p` 指向 `head`: ```graphviz digraph linked_list { node [shape=record]; rankdir=LR; p [label = "p"] head [label = "head(NULL)"] p->head } ``` 首先發現 `p` 的取值就是 `NULL` 了,所以直接跳出迴圈,將 `p` 指向的值改為 `&node0`: ```graphviz digraph linked_list { node [shape=record]; rankdir=LR; p [label = "p"] head [label = "head"] node0 [label = "<name>node0|<n>uninitialized"] p->head head->node0:name } ``` 接下來把 `node0` 的 `next` 改成 `NULL`: ```graphviz digraph linked_list { node [shape=record]; rankdir=LR; p [label = "p"] head [label = "head"] node0 [label = "<name>node0|<n>NULL"] p->head head->node0:name } ``` :::success 延伸問題: 在現有的程式碼基礎上,撰寫合併排序操作 ::: 這裡無法利用雙向鏈結串列的特性將子串列連接起來,所以使用陣列來實作堆疊。而 Linux 中對於 `list_sort` 的實作運用了鏈結串列可以在中間做修改的特性,所以不在堆疊頂部做合併並不會有什麼負面效果,而如果使用陣列的話不在頂部合併就要把頂部的東西往下面搬,成本較高,所以這裡使用的實作方式是我未引入 `list_sort` 方法時的[實現方式](https://hackmd.io/@rota1001/linux2025-homework1#q_sort)。也就是建立一個堆疊,用 `count` 的低位連續的 `1` 數量決定這次要合併幾次。 ```cpp static inline void merge_sort(list_t *l) { if (!l->head) return; list_item_t *stack[32]; int count = 0, sp = 0; for (list_item_t *head = l->head, *safe; head && ((safe = head->next) || 1); head = safe) { head->next = NULL; stack[sp++] = head; for (int tmp = count; tmp & 1; tmp >>= 1, sp--) stack[sp - 2] = merge_two(stack[sp - 1], stack[sp - 2]); count++; } for (;sp > 1 ;sp--) stack[sp - 2] = merge_two(stack[sp - 1], stack[sp - 2]); l->head = stack[0]; } ``` 另外,`merge_two` 的部份模仿 Linux ,不使用間接指標,而是使用結構單純的分支,以降低編譯後的指令數。這個部份在我的 [lab0](https://hackmd.io/@rota1001/linux2025-homework1#%E6%95%88%E8%83%BD%E6%AF%94%E8%BC%83) 有做過底層的分析。 ```cpp static inline list_item_t *merge_two(list_item_t *l, list_item_t *r) { list_item_t *head = NULL, **indir = &head; for (;;) { if (!l) { *indir = r; break; } if (!r) { *indir = l; break; } if (l->value <= r->value) { *indir = l; indir = &l->next; l = l->next; } else { *indir = r; indir = &r->next; r = r->next; } } return head; } ``` ### 測驗 2 `remove_free_tree` 是要找到特定的 `target`,並且將它從這棵樹中移除。 方法是如果它左右子樹都有的話,那麼拿左子樹中最右邊的節點(一直往右走到最底)來代替 `target`,如果只有一邊有子樹的話就拿那個子樹的跟節點來代替 `target`,否則 `target` 就會設成 `NULL`(是 `target` 所在位置的值被設成 `NULL`) :::success 解釋上方程式碼運作的原理,並嘗試補完程式碼,使其可獨立執行,應提供相關的測試程式碼 ::: 這裡只要去實作 `find_free_tree` 這個函式就好,它在二元搜尋樹上找到大於等於 `target` 的 `size` 中最小的節點所在的位置(指向它的指標)。 ```cpp block_t **find_free_tree(block_t **root, block_t *target) { if (!(*root)) return root; if ((*root)->l && ((*root)->l->size >= target->size)) return find_free_tree(&(*root)->l, target); if ((*root)->size >= target->size) return root; return find_free_tree(&(*root)->r, target); } ``` 這裡如果發現所有的節點都比 `target` 的 `size` 還要小的話,就會回傳一個指向值為 `NULL` 的指標的指標,至於這裡我沒有在 `remove_free_tree` 中做意外處理,因為他在 `remove_free_tree` 中沒有任何取得被移除的節點的途徑,所以我判斷這個函式的使用場景是已經確定有這個合法節點的情況下了。由於介面已經被定出來了,我沒有去修改,然而我認為界面應該要這樣定比較好: ```cpp block_t **find_free_tree(block_t **root, size_t target_size); void remove_free_tree(block_t **target) ``` `find_free_tree` 可以傳入 `size`,找到指向一個指向 `size` 大於等於 `target_size` 的節點的指標的指標。因為在 `find_free_tree` 中已經找到目標節點位於什麼位置了,所以 `remove_free_tree` 可以不用每次都呼叫 `find_free_tree` 從根節點開始找,並且外部也能透過 `find_free_tree` 獲得目標節點(它只是被從樹中移除,可以做後續操作)。 這裡為了測試寫了 `insert_tree` 和 `print_tree` 兩個函式,分別能做二元搜尋樹的插入和中序輸出節點的 `size`: ```cpp void print_tree(block_t *node) { if (!node) return; print_tree(node->l); printf("%ld ", node->size); print_tree(node->r); } void insert_tree(block_t **root, size_t size) { if (!(*root)) { *root = malloc(sizeof(block_t)); if(!(*root)) return; (*root)->size = size; (*root)->l = (*root)->r = NULL; return; } if ((*root)->size <= size) insert_tree(&(*root)->r, size); else insert_tree(&(*root)->l, size); } ``` 然後我用某種順序插入 0 到 9,然後一個一個刪除: ```cpp int main() { block_t *root = NULL; for (int i = 0; i < 10; i++) { insert_tree(&root, (i * 3) % 10); } print_tree(root); puts(""); block_t new_block; for (int i = 0; i < 10; i++) { new_block.size = (i + 4) % 10; remove_free_tree(&root, &new_block); print_tree(root); puts(""); } } ``` 結果如下: ``` 0 1 2 3 4 5 6 7 8 9 0 1 2 3 5 6 7 8 9 0 1 2 3 6 7 8 9 0 1 2 3 7 8 9 0 1 2 3 8 9 0 1 2 3 9 0 1 2 3 1 2 3 2 3 3 ``` 可以發現呈現預期行為。 :::success 參閱 memory-allocators,針對你補完 (及後續改進) 的記憶體配置器,進行效能評比,並解讀其效能表現 ::: 我在 github 上面開了一個 respository > [custom-memory-allocator](https://github.com/rota1001/custom-memory-allocator) 首先我的 `allocator` 的結構體是這樣定義的: ```cpp typedef struct allocator { size_t size; block_t *root; char mem[0]; } allocator_t; ``` `size` 指的單純給 `mem` 的空間有多少(不計算 header),`root` 是二元搜尋樹的根,`mem` 是長度為 0 的陣列,但在這裡這樣的好處是能讓 `mem` 的空間取決於我分配給他的空間,而且它本身是不耗空間的。 然後我的 `block` 在原始程式的基礎做修改: ```cpp typedef struct block { size_t size; size_t prev_size; struct block *l, *r; char mem[0]; } block_t; ``` 這裡參考了一些 glibc 中的 `malloc` 的結構體的概念,紀錄 `prev_size` 可以讓我們計算出上一個區塊在哪裡,size 可以讓我們計算下一個區塊在哪裡,而我讓 `prev_size` 和 `size` 都向 8 位元組對齊,所以可以把他的最低位拿來紀錄現在有沒有在使用。`l` 和 `r` 依然是代表二元搜尋樹上的左右子節點,而 `mem` 和上面是一樣的概念。不一樣的地方是,這裡的 `size` 是紀錄包含 header 的大小,因為多數的使用場景是要拿來計算上一個和下一個區塊的地址,這樣可以少一些運算。 我的界面是這樣定的: ```cpp /** * alloc_create - Create a new allocator with a * maximum memory size, which will be aligned up to 8 bytes * @size: maximum memory size */ allocator_t *alloc_create(size_t size); /** * alloc_alloc - Allocate memory with allocator and size * If there is not any invalid block, it will return NULL * @allocator: the allocator * @size: the size of memory needed */ void *alloc_alloc(allocator_t *allocator, size_t size); /** * alloc_free - Free an allocated memory * @allocator: the allocator * @ptr: pointer to the memory region */ void alloc_free(allocator_t *allocator, void *ptr); ``` - `alloc_create` 創建一個指定總共記憶體大小的 `allocator`。 - `alloc_alloc` 使用 `allocator` 來分配一塊指定大小的記憶體區塊。 - `alloc_free` 釋放一塊記憶體區塊。 `alloc_create` 首先會把 `size` 做 `ALLIGN_UP`,然後分配空間。我會在最後面塞一塊一直 `inuse` 的 `block_t`,讓每個區塊去存取他的 `NEXT_BLOCK` 都不會出問題。然後把整個記憶體變成一個區塊,放進二元搜尋樹中: ```cpp allocator_t *alloc_create(size_t size) { size = ALIGN_UP(size); allocator_t *allocator = malloc(sizeof(allocator_t) + sizeof(block_t) * 2 + size); allocator->size = size; allocator->root = NULL; block_t *first_block = (block_t *) (allocator + 1); block_t *last_block = (block_t *) (first_block->mem + size); first_block->size = size + sizeof(block_t); first_block->prev_size = 1; insert_node(&allocator->root, first_block); last_block->size = 1; return allocator; } ``` 接下來是 `alloc_alloc`,我首先會 `ALIGN_UP`,然後確保我申請到的可用空間可以至少塞下兩個 header 加上輸入的 `size`,讓我總是可以把它分成兩塊。接下來把原本的區塊從二元搜尋樹中移除,然後將切完剩下的插入: ```cpp void *alloc_alloc(allocator_t *allocator, size_t size) { size = ALIGN_UP(size); block_t **pblock = find_node_by_size(&allocator->root, size + sizeof(block_t) * 2); if (!pblock) return NULL; block_t *block = *pblock, *left_block = (block_t *) ((char *) block + size + sizeof(block_t)); remove_node(pblock); left_block->size = block->size - size - sizeof(block_t); block->size = (size + sizeof(block_t)) | 1; left_block->prev_size = block->size; insert_node(&allocator->root, left_block); return block->mem; } ``` 這裡來介紹一下這裡用到的巨集: ```cpp #define ALIGN_UP(size) (((size) + 7) & ~0x7) #define PREV_INUSE(block) (((block_t *) (block))->prev_size & 1) #define INUSE(block) (((block_t *) (block))->size & 1) #define CLEAR_USE_BIT(size) ((size) & ~0x1) #define NEXT_BLOCK(block) \ ((block_t *) ((char *) (block) + CLEAR_USE_BIT(((block_t *) (block))->size))) #define PREV_BLOCK(block) \ ((block_t *) ((char *) (block) - CLEAR_USE_BIT(((block_t *) (block))->prev_size))) #define NEXT_INUSE(block) (INUSE(NEXT_BLOCK(block))) #ifndef container_of #define container_of(ptr, type, member) \ ((type *) ((char *) (ptr) - offsetof(type, member))) #endif ``` - `ALIGN_UP` 向上對齊到 8 位元組 - `PREV_INUSE` 上一個區塊是否被使用中 - `NEXT_INUSE` 下一個區塊是否被使用中 - `INUSE` 現在區塊是否被使用中 - `CLEAR_USE_BIT` 將大小去除 `USE_BIT` - `NEXT_BLOCK` 下一個區塊 - `PREV_BLOCK` 前一個區塊 - `container_of` 將 `list.h` 中使用的巨集引入,可以從成員得到結構體 接下來是 `alloc_free`,首先會去看這個 `block` 的前後是否在使用中,如果不是在使用中就把他們從樹中移出,並且和 `block` 合併,最後把它插入樹中。 ```cpp void alloc_free(allocator_t *allocator, void *ptr) { block_t *block = container_of(ptr, block_t, mem); block_t **pblock; if (!PREV_INUSE(block)) { pblock = find_node(&allocator->root, PREV_BLOCK(block)); remove_node(pblock); block = merge(PREV_BLOCK(block), block); } if (!NEXT_INUSE(block)) { pblock = find_node(&allocator->root, NEXT_BLOCK(block)); remove_node(pblock); block = merge(block, NEXT_BLOCK(block)); } block->size = CLEAR_USE_BIT(block->size); NEXT_BLOCK(block)->prev_size = block->size; insert_node(&allocator->root, block); } ``` 這裡也對最初的在樹上找點與刪除節點的函數做了改良,首先是界面的部份: ```cpp static block_t **find_node(block_t **root, block_t *target); static block_t **find_node_by_size(block_t **root, size_t size); static void remove_node(block_t **node_ptr); ``` 使用 `find_node` 可以找到特定指標的節點在哪(指向他的指標在哪),`find_node_by_size` 可以找到大於等於 `size` 中最小的節點的位置,`remove_node` 是直接輸入指向指向節點的指標的指標。因為 `find_node` 必須要找到指標一樣的節點,所以二元搜尋樹做插入時會把指標也當作大小評判的標準,讓每個元素在大小關西上都是獨一無二的。另外,原本 `remove_node` 中多次呼叫 `find_node`,而且是從跟節點開始找,而當我直接找到節點去做刪除的話,就只要找一次就好,詳細實作可以看 Commit [f8279c3](https://github.com/rota1001/custom-memory-allocator/commit/f8279c34ae1418937eda5414de12258a96cccf4c) 測試的時候發現有些操作會造成 segmentation fault,原因是因為在 `alloc_alloc` 中,我讓區塊大小改變後沒有去改變下一個區塊的 `prev_size`,在 Commit [c056a09](https://github.com/rota1001/custom-memory-allocator/commit/c056a0963a410380421be641771f025dacdc5f35) 做了修正。 在 Commit [69326d7](https://github.com/rota1001/custom-memory-allocator/commit/69326d7caf1550623a2d73b027988a846a9e9f3a),我增加了比對效率的測試,會執行 1e6 次,隨機大小和操作的 alloc 或 free,以下是結果: ![custom-alloc](https://hackmd.io/_uploads/Bkflbqlnkl.png) 會發現我的 custom-alloc 比 malloc 慢,我將 alloc 和 free 拆開看時,會發現我的 alloc 是比 malloc 快的,比較慢的是在 free 的部份: ``` $ ./test.elf malloc: 278674908 custom-alloc: 207225006 0 ``` 在 Commit [20c760b](https://github.com/rota1001/custom-memory-allocator/commit/20c760b45e6c933109a7a9459494895badfd7e31),我在 `block` 結構體中加入了 `pparent`,紀錄他的父親指向它的那個指標的位置,讓我在 `free` 中可以不用去二元搜尋樹中搜尋,而且在 `insert` 的時候也不去比較指標大小,它帶來了 15% 左右的效率提昇,以下是它與 `malloc` 的比較: ![custom-alloc](https://hackmd.io/_uploads/HJA5U5l3yg.png) 雖然還是比較慢,但差距縮小了。 ### 測驗 3 首先簡單的填個空: - `GGGG`: `head->prev = prev` - `HHHH`: `list_entry(pivot, node_t, list)->value` - `IIII`: `list_entry(n, node_t, list)->value` - `JJJJ`: `pivot` - `KKKK`: `right` :::success 解釋上述程式碼運作原理 ::: 會有個由 `list_head *` 組成的陣列 `begin`,代表著是還沒被處理完的鏈結串列,並且是以一個堆疊的形式做插入和取出。每次會創建兩個空的鏈結串列 `left` 和 `right`,接下來從 `begin` 取出一個鏈結串列。如果這個鏈結串列只有一個元素就把它放進 `result` 中,否則將鏈結串列中的第 0 個元素當作 `pivot`,將此鏈結串列中大於 `pivot` 的元素放到 `right` 中,將小於等於 `pivot` 的元素放到 `left` 中,並依序將 `left`、`pivot`、`right` 放進 `begin` 中。這個時候可以發現,`begin` 中保證堆疊下面的元素一定小於等於堆疊上面的元素,所以最後會得到排序完的結果。 可以發現,這段程式碼並沒有充分發揮雙向鏈結串列的優勢,首先因為在排序中不需要維護雙向鏈結串列的性質,最後一個元素可以指向 `NULL`,而在裡面 `R` 這個指向最後一個元素的只有用到一個地方,而且實際上可以將 `L != R` 改成 `L->next != NULL` 而達到一樣的效果(如果讓最後一個元素總是指向 `NULL` 的話)。另外,參考 `list_sort` 的實作,可以利用雙向鏈結串列的 `prev` 來將子串列連接起來,以實現堆疊,這樣就不用在堆疊上面開空間了。 :::success 研讀〈[A comparative study of linked list sorting algorithms](https://dl.acm.org/doi/pdf/10.1145/225890.225893)〉,分析針對雙向鏈結串列的排序演算法考量,並運用 Linux 核心風格的 List API 予以實作 ::: 雙向鏈結串列可以將 `prev` 和 `next` 當作術的 `l` 和 `r` 來實作 tree sort,或者在排序時改成單向鏈結串列以減少維護的成本。至於快速排序法的部份,為了維持穩定排序,所以多用了一個鏈結串列 `equal` 來紀錄和 `pivot` 一樣大的元素。在以下,承襲這個概念,並且用迭代的方式實作快速排序法: ```cpp #define MAGIC (void *)1 void quick_sort(struct list_head *list) { int value; struct list_head *result = NULL, *left = NULL, *right = NULL, *equal = NULL, **presult = &result, **pleft = &left, **pright = &right, **pequal = &equal; list->prev->next = NULL; struct list_head *begin = NULL; list->next->prev = begin; begin = list->next; while (begin) { if (!begin->next || begin->next->prev == MAGIC) { *presult = begin; while (*presult) presult = &(*presult)->next; begin = begin->prev; continue; } value = list_entry(begin, node_t, list)->value; struct list_head *p = begin; begin = begin->prev; for (; p; p = p->next) { int n_value = list_entry(p, node_t, list)->value; if (n_value > value) { *pright = p; pright = &(p->next); } else if (n_value < value) { *pleft = p; pleft = &(p->next); } else { *pequal = p; pequal = &(p->next); } } *pright = *pequal = *pleft = NULL; if (equal && equal->next) equal->next->prev = MAGIC; if (right) { right->prev = begin; begin = right; } if (equal) { equal->prev = begin; begin = equal; } if (left) { left->prev = begin; begin = left; } left = right = equal = NULL; pright = &right; pleft = &left; pequal = &equal; } list->next = result; rebuild_list_link(list); } ``` 我將一個鏈結串列以 `pivot` 分為 `left`、`equal`、`right` 三個鏈結串列,對於每次從 `begin` 中取出元素時,會先去檢查他是不是 `equal`。如果它只有一個元素的話就是 `equal`,如果他有多個元素但是他的 `next->prev` 是 `MAGIC` 的話,那麼它也是 `equal`。如果他不是 `equal` 的話就取第 0 個元素為 `pivot`,並且開始分類。因為在 `equal` 中的元素是和原本順序一樣的,所以他是一種穩定排序法,而它也通過了原本題目中的 `list_is_ordered` 測試。 ## Q2 鏈結串列、雜湊表與位元運算 ### 測驗 1 先簡單的填個空: - `AAAA`: `list_entry` - `BBBB`: `list_del` - `CCCC`: `list_move_tail` - `DDDD`: `list_add` - `EEEE`: `list_splice` - `FFFF`: `list_splice_tail` :::success 解釋上方程式碼運作原理,改進 `random_shuffle_array` 也探討 `list_for_each_entry_safe` 建構的迴圈中,若將 `list_move_tail` 換為 `list_move`,為何會無法滿足 `stable sorting` ::: 首先回答為什麼用 `list_move` 就不是 `stable sorting`。因為 `list_move_tail` 是把一個節點移到目標鏈結串列的尾巴,所以依序的插入後,先後順序是不會變的,而 `list_move` 是從頭插入,所以在鏈結串列中的順序會和插入的順序相反,所以相同大小的元素的順序不保證和一開始一樣。 上面的那段程式碼,基本上就是〈[A comparative study of linked list sorting algorithms](https://dl.acm.org/doi/pdf/10.1145/225890.225893)〉提到的關於快速排序法的遞迴實作,這個在 [第一週 測驗 3](https://hackmd.io/@rota1001/linux2025-homework2#測驗-3) 已經予以討論,並且以迭代的方式改寫。 至於 `random_shuffle_array` 的改進可以使用 Fisher-Yates shuffle algorithm 去實作: ```cpp static inline void random_shuffle_array(uint16_t *operations, uint16_t len) { for (; len > 1; len--) { uint16_t choose = get_unsigned16() % len; if (choose == len - 1) continue; operations[choose] ^= operations[len - 1]; operations[len - 1] ^= operations[choose]; operations[choose] ^= operations[len - 1]; } } ``` ### 測驗 2 #### clz2 ```cpp static const int mask[] = {0, 8, 12, 14}; static const int magic[] = {2, 1, 0, 0}; unsigned clz2(uint32_t x, int c) { if (!x && !c) return 32; uint32_t upper = (x >> (16 >> c)); uint32_t lower = (x & (0xFFFF >> mask[c])); if (c == 3) return upper ? magic[upper] : 2 + magic[lower]; return upper ? clz2(upper, c + 1) : (16 >> (c)) + clz2(lower, c + 1); } ``` 對於輸入的 `x` 和 `c`,代表這個函式要計算 `x` 這個數字的最低 $2^{5-c}$ 個位元中的前綴 0。首先計算高位和低位各 $2^{4-c}$ 個位元,並且 4 個位元以下建表查詢。否則如果高位出現 1 就查詢高位中的前綴 0,否則就把高位的位數加上去,並且查詢低位中的前綴 0。 #### 整數開平方根 要思考範例程式的正確性,我先不一次到位,從改進程式的觀點來看這個程式。首先看一下以下程式碼: ```cpp uint64_t sqrti_simple(uint64_t x) { uint64_t m, y = 0; if (x <= 1) return x; int total_bits = 64; int shift = (total_bits - 1 - clz64(x)) & ~1; uint64_t p = 0; for(int i = (shift / 2); i >= 0; i--) { uint64_t b = (p << (i + 1)) + (1 << (2 * i)); if (x >= b) { x -= b; p += (1 << i); } } return y; } ``` 這個程式的目的很簡單,就是 [Digit-by-digit-Method](https://hackmd.io/@sysprog/sqrt#Digit-by-digit-Method) 的實作。定義 $P_i=\sum^{\frac{shift}{2}}_{j=i}a_j$,其中 $a_j$ 是 $2^j$ 或 $0$,取決於第 $j$ 個位元是不是 $1$。每次計算 $P$ 都會去檢查 $(P_{i+1}+2^i)^2\leq N$ 有沒有成立,其中 $(P_{i+1}+2^i)^2=P_{i+1}^2+2^{i+1}P_{i+1}+2^{2i}$,而我們會令 $X_i=N-P_i^2$,所以等同於去檢查 $2^{i+1}P_{i+1}+2^{2i}\leq X_{i+1}$ 是否成立。而關於 $X$ 的更新方式,很明顯的是: $$X_i=N-P_i^2=N-(P_{i+1}^2+2a_iP_{i+1}+a_i^2)=X_{i+1}-(2aiP_{i+1}+a_i^2)$$ 而這個程式有改進空間,就是位移的運算用太多了,於是我們可以用 $y_i = p_i \times 2^i$ 來替換,而當 `i` 為 0 的時候可以發現兩者相等。可以發現上面的 `(p << (i + 1))` 實際上是 $p_{i+1}\times2^{i+1}$,所以可以變成 `y`,而 `p += (1 << i)` 實際上是 $p_i = p_{i+1}+2^i$ ,同乘 $2^i$ 會變成 $2^ip_i=2^ip_{i+1}+2^{2i}$,又可以變成 $y_i=y_{i+1}2^{-1}+2^{2i}$,所以 `p += (1 << i)` 可以變成 `y = (y >> 1) + (1 << (2 * i))`,我貼近範例程式的實作會變成: ```cpp uint64_t sqrti_simple(uint64_t x) { uint64_t m, y = 0; if (x <= 1) return x; int total_bits = 64; int shift = (total_bits - 1 - clz64(x)) & ~1; for(int i = (shift / 2); i >= 0; i--) { uint64_t b = y + (1 << (2 * i)); y >>= 1; if (x >= b) { x -= b; y += (1 << (2 * i)); } } return y; } ``` 可以發現 $2^{2i}$ 可以被替換成 `m`,而且可以用簡單的位移搞定,並且 `i >= 0` 可以用 `m > 0` 來代替,所以可以得到原本的程式碼: ```cpp uint64_t sqrti(uint64_t x) { uint64_t m, y = 0; if (x <= 1) return x; int total_bits = 64; int shift = (total_bits - 1 - clz64(x)) & ~1; m = 1ULL << shift; while(m) { uint64_t b = y + m; y >>= 1; if (x >= b) { x -= b; y += m; } m >>= 2; } return y; } ``` :::success 解釋上述程式碼,並探討擴充為 $\lceil \sqrt{x}) \rceil$ (向上取整數) 實作,該如何調整並保持只用整數加減及位元操作 ::: 向上取整的整數平方根其實只是從找到最大的 $y$ 使得 $y^2\leq N$ 變成找到最大的 $y$ 使得 $(y-1)^2<N$ 而已,所以只要把 `x >= b` 改成 `x > b`,並且把結果加 1 就好。 ```cpp uint64_t sqrt_ceil(uint64_t x) { uint64_t m, y = 0; if (x <= 1) return x; int total_bits = 64; int shift = (total_bits - 1 - clz64(x)) & ~1; m = 1ULL << shift; while (m) { uint64_t b = y + m; y >>= 1; if (x > b) { x -= b; y += m; } m >>= 2; } return y + 1; } ``` #### 浮點數開平方根 :::success 參照計算機結構測驗題 C 及其注釋,以 C 語言 (搭配 GNU extension) 撰寫不依賴除法指令的 sqrtf,確保符合 IEEE 754 規格。對照 sqrtf, rsqrt, [reciprocal implementations in Arm assembly](https://github.com/picolibc/picolibc/issues/956) ::: :::warning 這裡最大的突破是我用建表的方式去改良 [reciprocal implementations in Arm assembly](https://github.com/picolibc/picolibc/issues/956) 中的 `mantissa_sqrt` ::: 我首先會說明它原始的程式的原理,然後再說明我怎麼做出改良的(我有對原始程式以盡量簡潔的方式先整理過,所以和他的會有一點不同)。 首先是 `custom_sqrtf`,它會將輸入 `x` 分為 `sign`、`exponent` 和 `mantissa` 三個部份。如果 `exponent` 是 0 的話就依照 `sign` 回傳 `+0` 或 `-0`,如果 `x` 是負數或 `NaN` 的話,那就回傳 `qNaN`。如果 `x` 是 `INF` 的話就回傳 `INF`。否則,就進入以下平方根的計算: ```cpp int32_t new_exp = ((int32_t)exponent - (127 << 23)) >> 1; new_exp = (new_exp + (126 << 23)) & (0xFF << 23); uint32_t new_mantissa = mantissa_sqrt(mantissa); if (((uint64_t) new_mantissa * new_mantissa + new_mantissa) < ((uint64_t) mantissa << 23)) new_mantissa++; u.i = new_exp + new_mantissa; return u.f; ``` 先前情提要一下,`mantissa` 指的是本來是 `1.XXXX` 的數字乘上 $2^{23}$,並且因為指數部份會捨去最低位,所以如果它的最低位是 1 的話就會把 `mantissa` 乘以 2,然後 `1X.XXX` 的平方根一樣會在 1 和 2 之間,所以還是可以當作一般的 `mantissa` 去處理。而這裡所有的數字都是和它原本在 `uint32_t` 中的位置對齊的(這裡的對齊是指 fix point 的那種對齊,而不是 align)。 它會先把 `exponent` 的部份減掉 127 ,除以 2 完後再加上 126,這裡有兩個有趣的點要解釋。首先是最高的 9 個位元在加減法一樣會符合群的性質,並且因為轉型成 `int32_t`,所以如果是在 GNU C 的規範下,`int32_t` 是 Arithmetic shift,也就是有號整數左移時會在左邊補 sign bit,所以在這些運算下可以把它當作一個 9 位元的有號整數來使用。另一個有趣的點是為什麼是加 126 而不是 127,這是因為 `mantissa` 在開完根號後,最高位一定有一個 1,代表 `1.XXXX`,所以剛好把這個 1 補上去就變 127 了。 接下來是用 `mantissa_sqrt` 計算出 `new_mantissa`,這部份等下再回來講,它會將 `mantissa` 開根號,然後一樣是以 $2^{23}$ 作為 fix point 的對齊。這個回傳結果是一個最大的 `x` 使得 $(x\times2^{-23})^2\leq mantissa\times 2^{-23}$ (因為這裡的 fix point 之後會影響計算,以下計算都以實際值來表示),而 `x + 1` 有可能讓結果更接近,所以可以列出這個不等式: \begin{aligned} &((x + 1)\times 2^{-23}) ^ 2 - mantissa \times 2^{-23}<mantissa \times 2^{-23} - (x \times 2^{-23})^2 \\ &\iff 2x^2+2x+1 < 2\times mantissa \times 2^{23} \\ &\iff x^2+x < mantissa \times 2^{23} \end{aligned} 所以可以說明這個不等式成立時,`new_mantissa` 應該要加 1。 接下來先要說明的是它原始的 `mantissa_sqrt` 實現: ```cpp uint32_t mantissa_sqrt(uint32_t x) { // sqrt of 9.23 fixed point number in the range [1,4) x <<= 6; uint32_t y = x << 1; uint32_t t = x + (x >> 1); if (t < (1u << 31)) { t += t >> 1; if (t < (1u << 31)) { x = t; y += y >> 1; } } for (uint32_t i = 2; i < 14; i += 1) { uint32_t t = x + (x >> i); t += t >> i; if (t < (1u << 31)) { x = t; y += y >> i; } } uint32_t hi = ((uint64_t) (x) * (y)) >> 32; y += y >> 1; y -= hi; return y >> 8; } ``` 這裡要先講一下他的核心概念,並且在推導的過程中,有別於原始推導,加入了 fix point 的計算以貼近程式碼,並且補充一些數學細節。 對於 $x > 1$,我都可以找到一個 $S\subseteq\mathbb{Z}$ 使得 $\prod_{i\in S}(1+2^{-i})$ 足夠接近 $x$,這裡用到了在 $k$ 很小的時候 $\ln(1+k) \approx k$ (這直接用泰勒級數的一次項來逼近),所以: $$\ln x \approx \ln(\prod_{i\in S}(1+2^{-i}))=\sum_{i\in S}\ln(1+2^{-i}) \approx \sum_{i \in S}2^{-i}$$ 可以發現選 $S$ 就是設法將 $\ln x$ 用二進位制的方式表示出來,所以它很顯然的可以找到一個夠近的估計值。並且如果 $1 \leq x < 2$,那 $0 \leq \ln x < \ln2 < 1$,所以只要考慮 $i>0$ 的部份就行了。 為了溝通方便,以下以數學符號表示的 $x$、$y$、$k$ 指的都是沒有經過 fix point 調整的值,譬如說 `x` 的初始值是 $x\times2^{23}$。首先,它會先把 `x` 左移 6,也就是把它變成 $x\times 2^{29}$,讓 `y` 變成 $x\times 2^{30}$,這些數字的考量下面會講。這個程式做的事情是找到 $S\subseteq \{1, 2, ...,14\}$ 使得 $\prod_{i\in S}(1+2^{-i})$ 逼近於 $\sqrt{\frac{2^{31}}{x\times 2^{29}}}$,而這件事情從用二進位制逼近 $\ln x$ 的角度來看就是從高位開始選擇這位是 1 還是 0。 所以這個時候,我們可以找到 $k$ 使得 $x \times 2^{29}\times k^2=2^{31}(1-\epsilon)$,其中 $\epsilon$ 是接近 0 的誤差值。而會發現當每次選取位元,也就是 `t < (1u << 31)` 的時候,`x` 都會乘兩次 $1+2^{-i}$,而 `y` 會乘一次,所以 `y` 會變成 $x \times2^{30}\times k$,簡單的做一下整理: \begin{aligned} &k = 2\sqrt{\frac{1-\epsilon}{x}} \\ &y = 2^{30}xk = 2^{31} \sqrt{x} \times\sqrt{1-\epsilon} \\ &2^{31}\sqrt{x} = \frac{y}{\sqrt{1-\epsilon}} \end{aligned} 然而,雖然 $1-\epsilon$ 已經很接近 1 了,但是因為他在分母,我們可以用差不多的成本進行更精確的估計,這時候可以使用二項式定理: \begin{aligned} &\frac{y}{\sqrt{1-\epsilon}}=y\sum_{i=0}^\infty (-1)^i\tbinom{-\frac{1}{2}}{i}\epsilon^i \approx y(1+\frac{\epsilon}{2}) & & \end{aligned} 最後 `x` 的值,我們叫它 $x'$,它等於 $2^{31}(1-\epsilon)$,所以 $\epsilon = 1-2^{-31}x$,代入上面式子: $$y(1+\frac{\epsilon}{2}) = y(1+\frac{1-2^{-31}x}{2}) = y(1+\frac{1}{2})-2^{-32}xy$$ 到這裡,可以發現下面這段程式碼就是將 $2^{31}\sqrt{x}$ 算出來: ```cpp uint32_t hi = ((uint64_t) (x) * (y)) >> 32; y += y >> 1; y -= hi; ``` 然後左移 8 就可以得到最後的 `mantissa`。 接下來講一下我的改進: ```cpp static const uint32_t table[] = {0, 954437177, 1374389535, 1696777203, 1902269252, 2019305101, 2081915509, 2114318375, 2130804226, 2139119553, 2143295481, 2145388031, 2146435455, 2146959455}; uint32_t mantissa_sqrt(uint32_t x) { x <<= 6; uint32_t y = x << 1; for (uint32_t i = 1; i < 14; i++) if (x <= table[i]) { x += x >> i; x += x >> i; y += y >> i; } uint32_t hi = ((uint64_t) (x) * (y)) >> 32; y += y >> 1; return (y - hi) >> 8; } ``` 可以發現,在不考慮溢位的情況下,當我的 `t` 滿足 `(t + (t >> 1)) + ((t + (t >> 1)) >> 1) >= (1u << 31)`,任何一個比 `t` 大的數字一定都滿足,所以可以對於最大的 `t` 使得 `(t + (t >> 1)) + ((t + (t >> 1)) >> 1) < (1u << 31)` 去建一個表,這樣只要檢查 t 有沒有比那個數小就好了。而且,因為第一次的 `t` 在計算中可能會溢位,所以原始程式中需要特判,這個程式就不用特判。然後會發現 `t` 這個變數不必要,全部用 `x` 來判斷就好。 #### 藉由查表的整數平方根 :::success 參照從 $\sqrt 2$ 的存在談開平方根的快速運算提到的「藉由查表的平方根運算」,以 C 語言實作,並比較不同手法的整數開平方根效能 ::: 以下為 [Integer square roots](https://github.com/mdickinson/snippets/blob/main/notebooks/Integer%20square%20roots.ipynb) 中的做法,首先證明一個定理(原文中先用其他部份限縮 $k$ 的範圍再予以討論,這裡以 $a \geq 2^k$ 的前提下討論,並且再後續程式碼的解讀中再去關注 $a$ 的範圍): $\forall a, n \in \mathbb{Z}^+,\ k\in \mathbb{Z}^+\cup\{0\},\ a \geq 2^k,$ $\text{if}\ (a-1)^2 < \lfloor \frac{n}{4^{k+1}} \rfloor < (a+1)^2,$ $\text{then}\ \lvert 2^k a + \lfloor \frac{n}{2^{k+2}a}\rfloor - \sqrt{n} \rvert <1$ 以下證明: $a-1 < \frac{\sqrt{n}}{2^{k + 1}} < a + 1$ $\Rightarrow \lvert \sqrt{n} - 2^{k+1}a \rvert < 2^{k + 1}$ $\Rightarrow 0 \leq n + 4^{k+1}a^2 - 2^{k+2}a\sqrt{n} < 4^{k+1}$ $\Rightarrow 0 \leq 2^ka+\frac{n}{2^{k+2}a} - \sqrt{n} < \frac{2^k}{a} \leq 1$ $\Rightarrow \lvert 2^k a + \lfloor \frac{n}{2^{k+2}a}\rfloor - \sqrt{n} \rvert <1$ 也就是說我們可以找合適的 $k$ 從一個和 $\sqrt{\frac{n}{4^{k + 1}}}$ 距離小於 1 的數推出一個和 $\sqrt{n}$ 距離小於 1 的數,接下來來解讀程式碼: ```python def isqrt_small(n): assert 0 < n < 2**64 # Shift left to make sure we're in range 2**62 <= n < 2**64 c = (n.bit_length() - 1) // 2 n <<= 62 - 2*c a = 1 + (n >> 62) # 2 <= a <= 4, a**2 ≈ n >> 60 a = (a << 1) + (n >> 59) // a # 8 <= a <= 16, a**2 ≈ n >> 56 a = (a << 3) + (n >> 53) // a # 128 <= a <= 256, a**2 ≈ n >> 48 a = (a << 7) + (n >> 41) // a # 2^15 <= a <= 2^16, a**2 ≈ n >> 32 a = (a << 15) + (n >> 17) // a # 2^31 <= a <= 2^32, a**2 ≈ n a -= a * a - 1 >= n return a >> 31 - c; ``` 首先它會將 `n` 乘上 $2^{62-2c}$,開完根號後再乘上 $2^{c-31}$,保證初始的 `n` 滿足 $2^{62}\leq n < 2^{64}$,這裡以這個前提開始討論。 一開始,因為 $1 \leq\frac{n}{2^{62}} < 2^2$,所以可以視為有個初始的 $a$ 為 1,這時候選一個最大的 $k=0$ 讓 $a\geq2^k$。這個時候滿足 $(a-1)^2<\lfloor \frac{2^{-60}n}{4^1}\rfloor < (a+1)^2$,所以可以找到一個和 $\sqrt{2^{-60}n}$ 距離小於 1 的 $2^0\times1+\lfloor \frac{2^{-60}n}{2^2\times 1}\rfloor$ ```python a = 1 + (n >> 62) ``` 接下來,同樣的我們知道 $(a - 1)^2<2^{-60}n<(a+1)^2$,其中,我們從上面產生出 $a$ 的地方可以產生一個算幾不等式導出 $a \geq \sqrt{2^{-60}n} \geq2$,所以可以選 $k = 1$。這個時候滿足 $(a-1)^2<\lfloor \frac{2^{-56}n}{4^2}\rfloor < (a+1)^2$,所以可以找到一個和 $\sqrt{2^{-56}n}$ 距離小於 1 的 $2^1a+\lfloor \frac{2^{-56}n}{2^3a}\rfloor$ ```python a = (a << 1) + (n >> 59) ``` 選 $k=3$,滿足 $(a-1)^2<\lfloor \frac{2^{-48}n}{4^4}\rfloor < (a+1)^2$,所以可以找到一個和 $\sqrt{2^{-48}n}$ 距離小於 1 的 $2^3a+\lfloor \frac{2^{-48}n}{2^5a}\rfloor$ ```python a = (a << 3) + (n >> 53) // a ``` 接下來分別選 $k=7, k=15, k=31$ 就可以得出一個和 $\sqrt{n}$ 距離小於 1 的 $a$。而且因為算幾不等式,$a \geq \sqrt{n}$,所以只要多考慮 $a > \sqrt{n}$ 的狀況就好。 然後接下來,它建了一個表,可以快速的找出一個 $a$ 使得 $(a-1)^2<\lfloor 2^{-48}n\rfloor < (a+1)^2$ 讓它可以少做兩次的 $a$ 的運算: ```cpp // isqrt64_tab[k] = isqrt(256*(k+65)-1) for 0 <= k < 192 static const uint8_t isqrt64_tab[192] = { 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, ... }; // integer square root of a 64-bit unsigned integer uint32_t isqrt64(uint64_t n) { if (n == 0) return 0; int lz = clz64(n) & 62; n <<= lz; uint32_t a = isqrt64_tab[(n >> 56) - 64]; a = (a << 7) + (n >> 41) / a; a = (a << 15) + (n >> 17) / a; a -= n < (uint64_t)a * a; return a >> (lz >> 1); } ``` 一樣把 $n$ 變成一個在 $[2^{62}, 2^{64})$ 中的整數再開始處理,然後會去查 `isqrt64_tab[(n >> 56) - 64]` 這一格。它代表的是 $\lfloor \sqrt{256\times(2^{56}n-64+65)-1} \rfloor = \lfloor \sqrt{256\times2^{56}n+255} \rfloor$。而這一個表需要估計出的是 $\sqrt{2^{-48}n}$,$\lfloor \sqrt{256\times2^{56}n+255} \rfloor$ 顯然較大,我們用數學計算一下兩者的差距(前情提要的一點是 $2^{-56}n \geq 64$): 顯然我們可以把前後雙方代換成 $\sqrt{256k+255}$ 和 $\sqrt{256k+x}$,其中 $k \geq 64$ 且 $0 \leq x < 256$,然後雙方差距為: $\sqrt{256k+255}-\sqrt{256k+x} = \frac{255-x}{\sqrt{256k+255}+\sqrt{256k+x}}\leq\frac{255-0}{\sqrt{245\times 64 + 255} + \sqrt{256\times 64 + 0}}\approx 0.992 < 1$ 所以 $\lfloor\sqrt{256k+255}\rfloor - \sqrt{256k+x} < 1$。令 $a = \lfloor\sqrt{256k+255}\rfloor$,因為 $a+1>\sqrt{256k+255}\geq \sqrt{256k+x}$,所以 $(a+1)^2>256k+x$。而另一邊移項可得 $a-1<\sqrt{256k+x}$,所以 $(a-1)^2<256k+x$。 綜上所述,可得 $(a-1)^2<\lfloor 2^{-48}n\rfloor < (a+1)^2$ ### 測驗 3 先簡單的填個空: - `AAAA`: `map->bits` - `BBBB`: `map->bits` - `CCCC`: `first->pprev` - `DDDD`: `n->pprev` - `EEEE`: `n->pprev` :::success 解釋上述程式碼運作原理,應提供測試程式碼 >針對 Linux 核心的 hash table 實作,提供更多圖例並揣摩 Linux 核心開發者 ::: 以下是測試程式碼,我將 1000000 以下所有經過雜湊函數後是某個隨機數的倍數的值都加進雜湊表,並且觀察行為是否正確。最後釋放空間,發現結果是沒錯的。 ```cpp int main() { srand(time(NULL)); int mod = rand() % 10000; map_t *map = map_init(10); for (int i = 0; i < 1000000; i++) { if (!(hash(i, 32) % mod)) map_add(map, i, malloc(10)); } for (int i = 0; i < 1000000; i++) { void *ret = map_get(map, i); assert((ret && !(hash(i, 32) % mod)) || (!ret && (hash(i, 32) % mod))); } map_deinit(map); } ``` 首先看一下結構體,`map_t` 中有一個陣列 `ht`,裡面放著一堆 `hlist_head`。`hlist_head` 中放的是一個節點 `first`。接下來是節點的結構體 `hlist_node`,裡面有一個指向下一個節點的指標和一個指向上一個節點中指向下一個節點的指標的指標。然後再把這個 `hlist_node` 內嵌到 `hash_key` 結構體中。 ```cpp struct hlist_head { struct hlist_node *first; }; struct hlist_node { struct hlist_node *next, **pprev; }; typedef struct { int bits; struct hlist_head *ht; } map_t; struct hash_key { int key; void *data; struct hlist_node node; }; ``` 以下圖解: ```graphviz digraph hash_table { node [shape=record]; rankdir=LR; map [label="{<ht0>ht[0]|<first0>first}|{ht[1]|<first1>first}|{ht[2]|first}"] node0 [label = "<name>node0|{<pprev>pprev|<next>next}"] node1 [label = "<name>node1|{<pprev>pprev|<next>next}"] node2 [label = "<name>node2|{<pprev>pprev|<next>next}"] map:first0->node0:name node0:next->node1:name map:first1->node2:name node1:pprev->node0:next node0:pprev->map:first0 node2:pprev->map:first1 } ``` 如果我要插一個東西到 `ht[0]` 的話,如果發現這個 `key` 不存在的話,它會插在開頭,像這樣: ```graphviz digraph hash_table { node [shape=record]; rankdir=LR; map [label="{<ht0>ht[0]|<first0>first}|{ht[1]|<first1>first}|{ht[2]|first}"] node0 [label = "<name>node0|{<pprev>pprev|<next>next}"] node1 [label = "<name>node1|{<pprev>pprev|<next>next}"] node2 [label = "<name>node2|{<pprev>pprev|<next>next}"] node3 [label = "<name>node3|{<pprev>pprev|<next>next}"] map:first0->node3:name node3:next->node0:name node0:next->node1:name map:first1->node2:name node1:pprev->node0:next node0:pprev->node3:next node3:pprev->map:first0 node2:pprev->map:first1 } ``` 使用 `pprev` 的原因可能是這樣,不是每個節點的前面都有 `prev` 和 `next`,譬如開頭的前面就只有 `first`,但是每個節點都有一個那來指向他的 `list_head *`,加了它之後刪除節點時就不用特判他是不是第一個節點。這裡的程式碼中,會用到這個功能的只有 `map_deinit`,這裡圖示一下如果移除 `node2` 並釋放會發生什麼: ```graphviz digraph hash_table { node [shape=record]; rankdir=LR; map [label="{<ht0>ht[0]|<first0>first}|{ht[1]|<first1>first}|{ht[2]|first}"] node0 [label = "<name>node0|{<pprev>pprev|<next>next}"] node1 [label = "<name>node1|{<pprev>pprev|<next>next}"] node3 [label = "<name>node3|{<pprev>pprev|<next>next}"] map:first0->node3:name node3:next->node0:name node0:next->node1:name node1:pprev->node0:next node0:pprev->node3:next node3:pprev->map:first0 } ``` 至於 `hash` 這個函數是雜湊函數: ```cpp #define GOLDEN_RATIO_32 0x61C88647 static inline unsigned int hash(unsigned int val, unsigned int bits) { /* High bits are more random, so use them. */ return (val * GOLDEN_RATIO_32) >> (32 - bits); } ``` `0x61C88647` 是 $\lceil\phi^{-2}\times2^{32}\rceil$,其中 $\phi$ 是黃金比例。會取這個數字的原因是因為隨機的數字乘上 $\phi^{-1}$ 後取小數點後的位數,他的分布會是均勻的,而乘兩次同樣也是均勻的。`hash` 做的事情就是取二進位制中的小數點後 `bits` 位。 :::success 進行《[The Art of Computer Programming](https://www-cs-faculty.stanford.edu/~knuth/taocp.html)》(vol 3) 一書section 6.4 的 exercise 8,也就是證明 Theorem S ::: 這本書是要錢的,但是在 [Linux 核心的 hash table 實作](https://hackmd.io/@sysprog/linux-hashtable#%E6%AF%94%E8%BC%83-golden-ratio-%E8%88%87%E9%9D%9E-golden-ratio)有貼出此定理內容的截圖。 首先,這裡用了一種連分數的遞迴表示法,對於連分數 $[a_1,a_2,...]$,我們定義: $\begin{cases} p_0 = 1 \\ p_1 = 0 \\ p_{k+1} = a_kp_k + p_{k-1} \text{, if}\ k \ge 1 \end{cases}$ $\begin{cases} q_0 = 0 \\ q_1 = 1 \\ q_{k+1} = a_kq_k + q_{k-1} \text{, if}\ k \ge 1 \end{cases}$ 這個時候 $[a_1, a_2, ..., a_{k-1}]$ 可以表示成 $\frac{p_k}{q_k}$ 對於 $k \ge 2$,以下用數學歸納法證明一下: 首先,$\frac{p_2}{q_2}=\frac{a_1}{1}=a_1=[a_1]$ 且 $\frac{p_3}{q_3}=\frac{a_2a_1+1}{a_2}=a_1+\frac{1}{a_2}=[a_1, a_2]$。 假設所有小於等於 $k, k \ge3$ 的情況都成立 $[a_1, a_2, ..., a_k] = [a_1,a_2,...,a_{k-1}+\frac{1}{a_k}]=\frac{p_{k-1}(a_{k-1}+\frac{1}{a_k})+p_{k-2}}{q_{k-1}(a_{k-1}+\frac{1}{a_k})+q_{k-2}}$ $=\frac{p_k+\frac{p_{k-1}}{a_k}}{q_k+\frac{q_{k-1}}{a_k}}=\frac{a_kp_k+p_{k-1}}{a_kq_k+q_{k-1}}=\frac{p_{k+1}}{q_{k+1}}$ 這樣就證明完了。 另外定義一下 $\{x\}$ 代表 $x$ 的小數部份, $\{x\}^+$ 代表 $\{x\}+1$ 接下來說明一下定理的內容: 對於無理數 $\theta$,我們有他的連分數表示 $[a_1, a_2, ...]$ 和對應的 $\{q\}_0^\infty$、$\{p\}_0^\infty$。對於每個正整數 $n$ 而言,有唯一一組 $k \ge 1,\ 1 \le r \le a_k,\ 0 \le s < q_k$ 使得 $n = rq_k + q_{k-1} + s$,而在插入 $\{n\theta\}$ 之前的狀態正好依序有以下 3 種區間: - $s$ 個長度為 $\{(-1)^k(rq_k+q_{k-1})\theta\}$ 的區間 - $n-q_k$ 個長度為 $\{(-1)^kq_k\theta\}$ 的區間 - $q_k-s$ 個長度為 $\{(-1)^k((r-1)q_k+q_{k-1})\theta\}$ 的區間 看到這裡,去想了一下後發現,如果不是我文字理解錯誤就是他有問題: 首先很顯然的,$\{\{a\}+\{b\}\}=\{a+b\}$,然後這 3 種區間總長度的總和應該要是 1,所以: $\{s\times\{(-1)^k(rq_k+q_{k-1})\theta\}+(n-q_k)\{(-1)^kq_k\theta\}+(q_k-s)\{(-1)^k((r-1)q_k+q_{k-1})\theta\}\}=0$ 也等同於是: $\{s(rq_k+q_{k-1})\theta+(n-q_k)q_k\theta+(q_k-s)((r-1)q_k+q_{k-1})\theta\}=0$ 等同於 $s(rq_k+q_{k-1})\theta+(n-q_k)q_k\theta+(q_k-s)((r-1)q_k+q_{k-1})\theta$ 是整數,於是我們把它展開: $s(rq_k+q_{k-1})\theta+(n-q_k)q_k\theta+(q_k-s)((r-1)q_k+q_{k-1})\theta$ $=sq_k\theta + (n-q_k)q_k\theta + ((r-1)q_k + q_{k-1})q_k \theta$ $=q_k\theta(s + n-q_k+(r-1)q_k+q_{k-1})$ $=q_k\theta(s+(rq_k+q_{k-1}+s)-q_k+(r-1)q_k+q_{k-1})$ $=q_k\theta((2r-2)q_k+2q_{k-1}+2s)$ $=2q_k\theta(n-q_k)$ 而因為 $2q_k(n-q_k)$ 是整數且 $\theta$ 不是整數,所以它唯一是整數的可能只有 $0$,而 $n-q_k=0$ 很明顯就不是恆真的,所以他們長度和不為 1,因此我認為這裡的敘述有問題。 這裡我想他在構造一種方式去證明任一個無理數 $\theta$,我們用 $\{\theta\},\{2\theta\},\{3\theta\}...$ 去切割 $[0, 1]$ 區間時,不會有超出 3 個種類的長度,目前還沒有什麼想法。