HackMD
2022q2 Homework (quiz2)
contributed by <ray90514>
測驗1
#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 。
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
EWMA的公式為
可化簡為
若 weight 為 2 的冪次方的話，可以使用 shift 加速計算。
也就是 include/linux/average.h 的作法。
WRITE_ONCE(e->internal, internal ? \
(((internal << weight_rcp) - internal) + \
(val << precision)) >> weight_rcp : \
(val << precision));
為了保留小數點以下的計算，先將值做 precision 的位移，而讀取EWMA時要補回來。
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 建立一個結構體儲存每次 EWMA 的值。然後宣告了對計算 EWMA 初始化、讀取、加入的函式。
struct ewma_##name { \
unsigned long internal; \
};
測驗2
#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位元有號數的改寫如下
#include <stdint.h>
int32_t max(int32_t a, int32_t b)
{
return a ^ ((a ^ b) & -(a < b));
}
測驗3
#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 。
延伸問題
- 使用
__builtin_ctz改寫後如下。
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% 。
- The image file may be corrupted
- The server hosting the image is unavailable
- The image path is incorrect
- The image format is not supported
- lib/math/gcd.c 實作兩種 gcd ，依照可不可以快速取得第一位非 0 的位數來分。
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 的冪的因數提出。
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)
a >>= 1;
if (a & r)
a += b;
a >>= 1;
測驗4
#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 二進位表示如下
而 -b = ~b + 1 ， -b 二進位表示如下
b 和 -b 做 and 後， EXP6 為 bitset & -bitset 或是 -bitset & bitset。
延伸問題
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 。
struct rem_node {
int key;
int index;
struct list_head link;
};
size = 1333;
struct list_head *heads = malloc(size * sizeof(*heads));
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;
}
先是計算兩數相除的商與餘數，接著處理小數點以下的部份。
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
/*
* 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
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 ，對照以下的示意
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 符合要求。
