# 2022q2 Homework (quiz2)
contributed by <`ray90514`>
## 測驗1
```c
#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
return (a >> 1) + (b >> 1) + (EXP1);
}
```
`(a >> 1) + (b >> 1)` 可以看成 `a / 2 + b / 2` ，但因為最低位會被捨去，所以若 a 和 b 最低位都為 1 ，相加會進位要補上 1。 `EXP1` 為 `a & b & 1` 。
```c
uint32_t average(uint32_t a, uint32_t b)
{
return (EXP2) + ((EXP3) >> 1);
}
```
`a + b` 可以寫成另外一種形式 `a ^ b + (a & b) << 1` ， `(a & b) << 1` 是進位， `a ^ b` 是相加後原位的結果。 `(a + b) / 2` 則可寫成 `(a + b) >> 1` 再套用前面的形式， `(a ^ b) >> 1 + (a & b) << 1 >> 1` ， 所以 `EXP2` 就是 `a & b` ， `EXP3` 為 `a ^ b` 。
### 延伸問題
以下為開 -O3 的輸出。
第一種實作
```
average:
.LFB23:
.cfi_startproc
endbr64
movl %edi, %eax
movl %esi, %edx
andl %esi, %edi
shrl %eax
shrl %edx
andl $1, %edi
addl %edx, %eax
addl %edi, %eax
ret
.cfi_endproc
```
第二種實作
```
average:
.LFB23:
.cfi_startproc
endbr64
movl %edi, %eax
xorl %esi, %eax
shrl %eax
andl %esi, %edi
addl %edi, %eax
ret
.cfi_endproc
```
第一個比第一個多了一個 `movl` ，一個 `shrl` ，一個 `andl` ，少一個 `xorl` 。
### EWMA
$s_n$ 為第 n 個 EWMA 的值， $a_n$ 為第 n 個測量到的值。
EWMA的公式為 $s_n = a_n / weight + (1 - 1 / weight) * s_{n-1}$
可化簡為 $s_n = (a_n + weight * s_{n - 1} - s_{n-1}) / weight$
若 weight 為 2 的冪次方的話，可以使用 shift 加速計算。
$s_n = (a_n + s_{n - 1} << (log_2(weight)) - s_{n - 1}) >> log_2(weight)$
也就是 [include/linux/average.h](https://github.com/torvalds/linux/blob/master/include/linux/average.h) 的作法。
```c
WRITE_ONCE(e->internal, internal ? \
(((internal << weight_rcp) - internal) + \
(val << precision)) >> weight_rcp : \
(val << precision));
```
為了保留小數點以下的計算，先將值做 `precision` 的位移，而讀取EWMA時要補回來。
```c
static inline unsigned long \
ewma_##name##_read(struct ewma_##name *e) \
{ \
BUILD_BUG_ON(!__builtin_constant_p(_precision)); \
BUILD_BUG_ON(!__builtin_constant_p(_weight_rcp)); \
BUILD_BUG_ON((_precision) > 30); \
BUILD_BUG_ON_NOT_POWER_OF_2(_weight_rcp); \
return e->internal >> (_precision); \
}
```
[include/linux/average.h](https://github.com/torvalds/linux/blob/master/include/linux/average.h) 建立一個結構體儲存每次 EWMA 的值。然後宣告了對計算 EWMA 初始化、讀取、加入的函式。
```c
struct ewma_##name { \
unsigned long internal; \
};
```
## 測驗2
```c
#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b)
{
return a ^ ((EXP4) & -(EXP5));
}
```
`max` 回傳的值只會是 a 或 b ，而已知 `a ^ (a ^ b) = b`, `a ^ 0 = a` ，所以 ``(EXP4) & -(EXP5)`` 的結果不是 `a ^ b` 就是 `0` ，我們可以讓 `a ^ b` & bitmask 獲得這兩種值，而 `EXP5` 為 a 和 b 的比較操作，所以 `EXP4` 是 `a ^ b` 。 `EXP5` 的結果只有 1 和 0 加上負號為 -1 和 0 ， -1 和 0 的二進位分別表示為全 1 和全 0 ，可以用作 bitmask 。所以如果當 a "較大"時要得到 a 的值，`(EXP4) & -(EXP5)` 的值須為 0，則 `EXP5` 比較的結果為 0，所以 `EXP5` 是 `a < b` 。
### 延伸問題
32位元有號數的改寫如下
```c
#include <stdint.h>
int32_t max(int32_t a, int32_t b)
{
return a ^ ((a ^ b) & -(a < b));
}
```
## 測驗3
```c
#include <stdint.h>
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
u /= 2, v /= 2;
}
while (!(u & 1))
u /= 2;
do {
while (!(v & 1))
v /= 2;
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (COND);
return RET;
}
```
先提取 2 的冪的最大公因數，接著把除數剩餘的 2 的冪的因數除掉，接下來重複做將被除數減掉除數，如果除數比被除數大就交換兩數，直到除數為 0 ，所以 `COND` 為 `v` 。因為 2 不再是公因數，可以把被除數有 2 的冪的因數除掉。`RET` 為剩下的被除數補上 2 的冪的最大公因數所以是 `u << shift` 。
### 延伸問題
1. 使用 `__builtin_ctz` 改寫後如下。
```c
uint64_t gcd64(uint64_t u, uint64_t v)
{
if (!u || !v) return u | v;
int shift;
shift = __builtin_ctz(u & v);
u >>= __builtin_ctz(u);
do {
v >>= __builtin_ctz(v);
if (u < v) {
v -= u;
} else {
uint64_t t = u - v;
u = v;
v = t;
}
} while (v);
return u << shift;
}
```
每次對 1000000 組隨機數字做 gcd ，總共測試十次得到以下結果，改寫後的版本比原本的快 80% 。
2. lib/math/gcd.c 實作兩種 gcd ，依照可不可以快速取得第一位非 0 的位數來分。
```c
unsigned long gcd(unsigned long a, unsigned long b)
{
unsigned long r = a | b;
if (!a || !b)
return r;
b >>= __ffs(b);
if (b == 1)
return r & -r;
for (;;) {
a >>= __ffs(a);
if (a == 1)
return r & -r;
if (a == b)
return a << __ffs(r);
if (a < b)
swap(a, b);
a -= b;
}
}
```
一樣是先提取 2 的冪的公因數，重複兩者相減，將被除數的 2 的冪的因數提出。
```c
unsigned long gcd(unsigned long a, unsigned long b)
{
unsigned long r = a | b;
if (!a || !b)
return r;
/* Isolate lsbit of r */
r &= -r;
while (!(b & r))
b >>= 1;
if (b == r)
return r;
for (;;) {
while (!(a & r))
a >>= 1;
if (a == r)
return r;
if (a == b)
return a;
if (a < b)
swap(a, b);
a -= b;
a >>= 1;
if (a & r)
a += b;
a >>= 1;
}
}
```
跟上面不一樣的是多了以下幾行。 `r` 是 2 的冪的公因數， `a & r` 相當於只看約掉後的值的最低位。
> If a and b are all odds, then gcd(a,b) = gcd((a-b)/2, b) = gcd((a+b)/2, b)
```c
a >>= 1;
if (a & r)
a += b;
a >>= 1;
```
## 測驗4
```c
#include <stddef.h>
size_t improved(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
size_t pos = 0;
uint64_t bitset;
for (size_t k = 0; k < bitmapsize; ++k) {
bitset = bitmap[k];
while (bitset != 0) {
uint64_t t = EXP6;
int r = __builtin_ctzll(bitset);
out[pos++] = k * 64 + r;
bitset ^= t;
}
}
return pos;
}
```
以下用 `b` 代替 `bitset` ，假設最低位的 1 是第 n 個， `b` 的二進位表示如下 $b_{64}, b_{63}, ..., 1(b_n), 0, 0, ..., 0$
`~b` 二進位表示如下
$\lnot b_{64}, \lnot b_{63}, ..., 0(b_n), 1, 1, ..., 1$
而 `-b = ~b + 1` ， `-b` 二進位表示如下
$\lnot b_{64}, \lnot b_{63}, ..., 1(b_n), 0, 0, ..., 0$
`b` 和 `-b` 做 and 後， $b_{64}$ 到 $b_{n + 1}$ 互為相反的值，所以結果都為 0 ， $b_{n}$ 兩者都為 1 所以結果為 1， $b_{n-1}$ 到 $b_{0}$ 兩者皆為 0 所以結果為 0 ， 這樣就只剩下 $b_{n}$ 為 1 ，也就是最低位元的 1 。 `EXP6` 為 `bitset & -bitset` 或是 `-bitset & bitset`。
### 延伸問題
```c
uint64_t get_bitset(int count)
{
uint64_t result = 0;
for (int i = 0; i < count; i++) {
result |= 1 << (rand() % 64);
}
return result;
}
```
使用以上程式碼產生不同分佈的數，對這些分佈比較兩種實作的時間，每個分佈有 1000000 個數。從結果可以看出在數據較稀疏的時候，使用 `__builtin_ctzll__` 的效果較好。
## 測驗5
概念是這樣的，做除法的時候，如果出現相同的餘數代表有循環小數，因此使用 hashtable 來儲存餘數。搭配 `find` 來尋找 hashtable 內的 value 。這裡的 hashtable 採用 array 搭配 linked list 。
```c
struct rem_node {
int key;
int index;
struct list_head link;
};
size = 1333;
struct list_head *heads = malloc(size * sizeof(*heads));
```
```c
static int find(struct list_head *heads, int size, int key)
{
struct rem_node *node;
int hash = key % size;
list_for_each_entry (node, &heads[hash], link) {
if (key == node->key)
return node->index;
}
return -1;
}
```
先是計算兩數相除的商與餘數，接著處理小數點以下的部份。
```c
for (int i = 0; remainder; i++) {
int pos = find(heads, size, remainder);
if (pos >= 0) {
while (PPP > 0)
*p++ = *decimal++;
*p++ = '(';
while (*decimal != '\0')
*p++ = *decimal++;
*p++ = ')';
*p = '\0';
return result;
}
struct rem_node *node = malloc(sizeof(*node));
node->key = remainder;
node->index = i;
MMM(&node->link, EEE);
*q++ = (remainder * 10) / d + '0';
remainder = (remainder * 10) % d;
}
```
一樣是兩數相除，將商加進字串中，以及將餘數 * 10 成為新的被除數。判斷是否有同樣的餘數出現，如果有代表是循環小數，將結果的字串依照儲存的位置補上 `'('` 和 `')'` ，如果沒有則將餘數跟出現的位置加進 hashtable ，直到除盡或是有循環小數。 `pos` 是循環小數開始出現的位置，所以如果有循環小數，要先將非循環小數的結果加到原本的字串， `PPP` 是 `pos--` 。
`MMM(&node->link, EEE)` 是要將節點加入至 hashtable ，所以 `MMM` 是 `list_add` 或 `list_add_tail` ，而 `EEE` 是被 array 儲存的 linked list 開頭，而根據 `find` 函式， hash 的值為 `key % size` ，所以是 `heads[remainder % size]` 。
## 測驗6
```c
/*
* ALIGNOF - get the alignment of a type
* @t: the type to test
*
* This returns a safe alignment for the given type.
*/
#define ALIGNOF(t) \
((char *)(&((struct { char c; t _h; } *)0)->M) - (char *)X)
```
`struct { char c; t _h; }` 這個結構體內的 `_h` 會因為 `c` 要做 alignment ，所以將 _h 的位址與結構體開頭的位址相減，可以得知 alignment 的狀況。 所以 `M` 為 `_h` ， `&((struct { char c; t _h; } *)0)->M` 是 `_h` 的位址， `X` 為結構體的起始位址所以是 0 ，因為是以 byte 為單位，所以轉型成指向 `char` 的指標，相減之後會獲得 t 是以幾個 byte 來對齊的。
## 測驗7
```c
static inline bool is_divisible(uint32_t n, uint64_t M)
{
return n * M <= M - 1;
}
static uint64_t M3 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 3 + 1;
static uint64_t M5 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 5 + 1;
int main(int argc, char **argv)
{
for (size_t i = 1; i <= 100; i++) {
uint8_t div3 = is_divisible(i, M3);
uint8_t div5 = is_divisible(i, M5);
unsigned int length = (2 << KK1) << KK2;
char fmt[9];
strncpy(fmt, &"FizzBuzz%u"[(9 >> div5) >> (KK3)], length);
fmt[length] = '\0';
printf(fmt, i);
printf("\n");
}
return 0;
}
```
`length` 是要從 `"FizzBuzz%u"` 複製的長度，如果是 3 的倍數 length 為 4，如果是 5 的倍數 length 為 4 ， 如果同時是 3 和 5 的倍數 length 為 8，都不是的話 length 為 2 。 所以 `KK1` 和 `KK2` 為 `div3` 和 `div5` 。 `(9 >> div5) >> (KK3)` 是 offset ，對照以下的示意
```c
string literal: "fizzbuzz%u"
offset: 0 4 8
```
若是 3 的倍數要輸出 `"fizz"` ， offset 為 0，若是 5 的倍數要輸出 `"buzz"` ， offset 為 4 ，若同時是 3 和 5 的倍數要輸出 `"fizzbuzz"` ， offset 為 0 ，若甚麼都不是 offset 為 8 。 `(9 >> 1) >> (KK3) = 4, KK3 = 0` ， `(9 >> 0) >> (kk3) = 0, kk3 >= 4` ， `(9 >> 1) >> (KK3) = 0, KK3 >= 4` ， `KK3` 為 `div3 << 2` 和 `div3 * 4` 符合要求。
###### tags: `linux2022`
Sign in with Wallet
Connect another wallet