Try   HackMD

2022q1 Homework2 (quiz2)

contributed by < qwe661234 >

測驗題目

測驗一

1-1

  1. a is even, b is odd => average = a / 2 + b / 2
  2. a is odd, b is even => average = a / 2 + b / 2
  3. a is even, b is even => average = a / 2 + b / 2
  4. a is odd, b is odd => average = a / 2 + b / 2 + 1

因此, 我們需要加上一個修正值, 而這個修正值在 a is odd 且 b is odd 時是 1, 而其他情況下是 0。
我們可以觀察到當數字是奇數時, 最後一個 bit 是 1, 而偶數是 0, 因此我們可以利用這個特性, 當 a, b 都是奇數時, a & b 結果的最後一個 bit 是 1, 而其他情況是 0, 接著再 &1 來取出最後一個 bit。

EXP1 = a & b & 1

#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a >> 1) + (b >> 1) + (a & b & 1);
}

1 - 2

我們知道 a + b 可以寫成 a ^ b + ((a & b) << 1)

  • a ^ b : a + b 未加上 carry 的結果
  • (a & b) << 1 : carry 的結果, 因為 carry 要加到下一位, 所以要左移

所以 (a + b) / 2 => (a + b ) >> 1 => (a ^ b + ((a & b) << 1)) >> 1 => (a & b) + ((a ^ b) >> 1)

EXP2 = a & b
EXP3 = a ^ b

#include <stdint.h>
uint32_t average(uint32_t a, uint32_t b)
{
    return (a & b) + ((a ^ b) >> 1);
}

比較上述實作在編譯器最佳化開啟的狀況,對應的組合語言輸出,並嘗試解讀

  • mov dest src: 將資料從 src 傳送一份到 dest
  • and src dest: dest = dest & src
  • xor src dest: dest = dest ^ src
  • shr dest: dest = dest >> 1
  • retq: 回傳
0x0000000000001180 <+0>:	endbr64 
0x0000000000001184 <+4>:	mov    %edi,%eax
0x0000000000001186 <+6>:	mov    %esi,%edx
0x0000000000001188 <+8>:	and    %esi,%edi
0x000000000000118a <+10>:	shr    %eax
0x000000000000118c <+12>:	shr    %edx
0x000000000000118e <+14>:	and    $0x1,%edi
0x0000000000001191 <+17>:	add    %edx,%eax
0x0000000000001193 <+19>:	add    %edi,%eax
0x0000000000001195 <+21>:	retq  

兩個 mov 分別傳送一份 a, b 到另外兩個暫存器
兩個 and 對應到 a & b & 1, 先將 a & b 的結果存下來, 再將 a & b 後的結果 & 1
兩個 shr 分別對應 a >> 1 和 b >> 1
第一個 add 是作 a >> 1 + b >> 1
第二個 add 把前一個 add 的結果與 (a & b & 1) 的結果相加

0x0000000000001180 <+0>:	endbr64 
0x0000000000001184 <+4>:	mov    %edi,%eax
0x0000000000001186 <+6>:	and    %esi,%edi
0x0000000000001188 <+8>:	xor    %esi,%eax
0x000000000000118a <+10>:	shr    %eax
0x000000000000118c <+12>:	add    %edi,%eax
0x000000000000118e <+14>:	retq 

mov 先傳送一份 a 或 b 到另一個暫存器, 者邊先假設 a
and 對應到 a & b
xor 對應到 a ^ b
shr 對應到 (a ^ b) >> 1, 將上一個 xor 後的結果做右移
add 是把 a & b 的結果與 (a ^ b) >> 1 的結果相加

測驗二

利用 xor 的特性:

  • x ^ x = 0
  • x ^ 0 = x

我們可以利用這兩個特性, 若 ((EXP4) & -(EXP5)) 結果是 a ^ b, 我們可以得到 b (因為 a ^ a ^ b = b), 若 ((EXP4) & -(EXP5)) 結果是 0, 我們可以得到 a (因為 a ^ 0 = a)。

上述例子我們可以看出, EXP4 = a ^ b, 接著分析 EXP5, 當 EXP4 = a ^ b 時, 若 a > b, 則 EXP5 的結果要是 0, 這樣 ((a ^ b) & 0) = 0, 反之, a < b 時, EXP5 的結果要是 1, 這樣 ((a ^ b) & -1) = a ^ b , 所以得 EXP5 = a < b

EXP4 = a ^ b
EXP5 = a < b

#include <stdint.h>
uint32_t max(uint32_t a, uint32_t b){
    return a ^ ((a ^ b) & -(a < b));
}

針對 32 位元無號/有號整數,撰寫同樣 branchless 的實作

原本認為有號正數無法使用此 function, 因為在做 a < b 比較時, 比較的部份是使用減法, 應該會發生 overflow, 但設法輸入會發生 overflow 的值結果仍正確, 於是分析 a < b 的 assembly code 來找出原因。

cmpl	-8(%rbp), %eax
setl	%al
movzbl	%al, %eax

cmpl %eax(src2), %edx(src1)

cmpl 指令是用來比較 a (存於 register %edx) 和 b (存於 register %eax), 先做 a - b, 然後用 a - b 的結果來設定 SF(sign flag)OF(overflow flag), a - b 的果會在設定為 flag 後捨棄。

  • SF set if (a - b) < 0 (as signed)
  • OF set if two’s-complement (signed) overflow
    (a > 0 && b < 0 && (a - b) < 0) || (a < 0 && b > 0 && (a - b) > 0)

setl %al

setl 指令是將 SF ^ OF 的結果存於 register %al 中

思考 overflow 情況

  1. a = INT32_MAX, b = - 1

    a - b < 0 -> SF = 1
    a > 0, b < 0, t < 0(因為 sign bit 變成 1) -> OF = 1
    SF ^ OF = 0 (a < b = 0)

  2. a = -2, b = INT32_MAX

    a - b > 0 -> SF = 0
    a < 0, b > 0, t > 0(因為 sign bit 變成 0) -> OF = 1
    SF ^ OF = 1 (a < b = 1)

  3. a = 1, b = INT32_MIN

    a - b < 0 -> SF = 1
    a > 0, b < 0, t < 0(因為 sign bit 變成 1) -> OF = 1
    SF ^ OF = 10 (a < b = 0)

  4. a = INT32_MIN, b = 1

    a - b > 0 -> SF = 0
    a < 0, b > 0, t > 0(因為 sign bit 變成 0) -> OF = 1
    SF ^ OF = 1 (a < b = 1)

因為 a < b 的結果由 signed flag 和 overflow flag 共同控制, 因此即使 signed integer 發生 overflow, 結果仍然正確, 故可以和 unsigned integer 用同一個 max fucniton

int32_t max(int32_t a, int32_t b)
{
    return a ^ ((a ^ b) & -(a < b));
}

reference:

http://www.cs.cmu.edu/afs/cs/academic/class/15213-s14/www/lectures/06-machine-control.pdf

測驗三

if (!u || !v) return u | v;

用來判斷以下三種情形

  1. u = 0 && v = 0 -> gcd(0, 0) = u | v = 0
  2. u = 0 && v != 0 -> gcd(0, v) = u | v = v
  3. u != 0 && v = 0 -> gcd(u, 0) = u | v = u
int shift;
for (shift = 0; !((u | v) & 1); shift++) {
    u /= 2, v /= 2;
}

因為 u 和 v 都是偶數, 則 gcd(u, v) = 2 * gcd(

u/2,
v/2
)
shift 是用來記錄除 2 多少次, 等等 return value 要乘回來
for loop 中 !((u | v) & 1) 是用來判斷 u 和 v 是不是都為偶數, 因為偶數的最後一個 bit 會是 0, 所以可以看成 !(u & 1) & !(v & 1), 合併後可得 !((u | v) & 1)

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

如果 u 是偶數, v 是奇數 gcd(u, v) = gcd(

u/2, v), 第一個 while loop 用來將 u 持續除以 2, 直到 u 為奇數。

如果 u 是奇數, v 是偶數 gcd(u, v) = gcd(u,

v/2), 在 do-while loop 中的第一個 while loop 用來將 v 持續除以 2, 直到 v 為奇數, 接著做輾轉相除法, 重複此步驟直到 v = 0 時, u 為最大公因數, 所以 do-while 的中止條件為 v = 0。

最後 return 時, 要記得把剛剛 shift 乘回來, 每除一次 2, shift++, 所以

u(2shift) = u << shift

COND = v
RET = u << shift

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 (v);
    return u << shift;
}

在 x86_64 上透過 __builtin_ctz 改寫 GCD,分析對效能的提升

__builtin_ctzl

此 funciton 會從最尾端的 bit 開始算起, 直到 bit 不為 0 停止, 並回傳 0 的個數, 透過 __builtin_ctzl 的結果來 shift, 可以做到和用 while loop 不斷除 2 相同的效果, 因為從最尾端的 bit 算起, 每遇到一個 0, 就會除以 2, 也就是 >> 1, 直到最尾端的 bit 不為 0。

while (!(u & 1))
    u /= 2;
u >> __builtin_ctzl(u)

上述這兩段程式碼效果是等價的。
而 shift 只須紀錄 __builtin_ctzl(u)__builtin_ctzl(v) 較小者, 者代表 u 和 v 都同時除以 2 的情況, 等同於 gcd64 的 shift。

uint64_t gcd64_ctz(uint64_t u, uint64_t v)
{
    if (!u || !v) return u | v;
    int shiftU, shiftV, shift;
    shiftU = __builtin_ctzl(u);
    shiftV = __builtin_ctzl(v);
    u = u >> shiftU;
    v = v >> shiftV;
    shift = shiftU > shiftV ? shiftV : shiftU;
    while (v){
        v = v >> __builtin_ctzl(v);
        if (u < v)
        {
            v -= u;
        }
        else
        {
            uint64_t t = u - v;
            u = v;
            v = t;
        }
    };
    return u << shift;
}

測試程式

隨機產生 100000 對 64bit 的 unsigned int 作為 function gcd64 的 input, 重複 10 次取平均。

uint64_t rand_64() {
    uint64_t r = rand(); 
    r = r << 30 | rand();
    return r;
}

int main() {
    srand(time(0));
    for(int j = 0; j < 10; j++) {
        clock_t start, end;
        start = clock();
        for(int i = 0; i < 100000; i ++){
            gcd64(rand_64(), rand_64());
        }
        end = clock();
        double diff = end - start;
        printf("%f\n" , diff / CLOCKS_PER_SEC);
    }
}

測試結果

效能提升約 50 %
去除第一次的結果取平均, 第一次結果接近其他 9 次結果的 2 倍, 故去除。
未使用 __builtin_ctzl: 0.069411 sec
使用 __builtin_ctzl: 0.035895 sec

測驗四

原理是 iterate 整個 bitmap array, 找出 array 中每個 element 的 bit 為 1 的位置回傳, 外層 for loop 是用來 iterate 整個 bitmap array, 內層 for loop 是用來 iterate element 中的 64 個 bit, 並紀錄 bit 為 1 的位置。

size_t naive(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
    size_t pos = 0;
    for (size_t k = 0; k < bitmapsize; ++k) {
        uint64_t bitset = bitmap[k];
        size_t p = k * 64;
        for (int i = 0; i < 64; i++) {
            if ((bitset >> i) & 0x1)
                out[pos++] = p + i;
        }
    }
    return pos;
}

改進的部份改用__builtin_ctzll 來實做, 因為 ctz 會計算最尾端 1 後面有幾個 0, 因此再做完一次後, 要消除最尾端的 1, 這樣才能再對 element 做 __builtin_ctzll

bitset & -bitset 可以使得最後一個為 1 的 bit 留下, 而其他是 0, 例如 010(2) & 110(-2) = 010, 將此結果與 bitset 做 xor 可以消除最尾端的 1。

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 = bitset & -bitset;
            int r = __builtin_ctzll(bitset);
            out[pos++] = k * 64 + r;
            bitset ^= t;                           
        }
    }
    return pos;
}

真實使用案例 Bloom Filter

Bloom Filter(布隆過濾器)由 Burton Howard Bloom 在 1970 構思出來,用來測試一個元素是否「可能存在」或是「絕對不存在」在容器中特定集合中, "可能存在"是因為有機會出現假陽性(false positive), 但絕不會有假陰性(false negative)。

Bloom Filter 的功能其實和 hash table 很相似, 但 Bloom Filter 可在 O(1) 時間複雜度驗證成員是否存在,且只需要相對少的儲存空間, 不過缺點是會有假陽性(false positive) 的可能存在。

Bloom Filter 實作:

  1. 大小為 n 的 bit array
  2. m 個 hash fucniton
  • 新增: 假設要新增 key 時, 透過 m 個 hash function 轉換會得到 m 個 index, 將 bit array 中的 m 個 index 都設成 1 。
  • 查詢: 要查尋某個 key 時, 一樣先透過 m 個 hash function 轉換會得到 m 個 index, 如果 bit array 中的這 m 個 index 都是 1, 則代表 key 「可能存在」, 如果非 m 個 index 都為 1, 則表示 key 絕對不存在」。

reference:

https://medium.com/@Kadai/資料結構大便當-bloom-filter-58b0320a346d
https://rust-algo.club/collections/bloom_filter/

設計實驗,檢驗 ctz/clz 改寫的程式碼相較原本的實作有多少改進?應考慮到不同的 bitmap density;

參考圖片以 64bit 大小為 16 的 bit array 模仿圖片中的 bit 排列, 比較原本 function 以及以 __builtin_ctzll 改進後的 funciton 的效能 。

naive: 0.001063 (second)
imporved_ctz: 0.0004879 (second)

naive: 0.001940 (second)
imporved_ctz: 0.000510 (second)

naive: 0.002162534 (second)
imporved_ctz: 0.001259721 (second)

naive: 0.002204634 (second)
imporved_ctz: 0.001884666 (second)

由實驗可知, 在 1 越頻繁出現的圖片中, 改進後的 function 效能會越接近原本的 function, 但改進後的 function 效能仍較佳。

思考進一步的改進空間

對原本的版本進行改進, 原本一次只看 1 個 bit, 改進成一次可以看 2 個 bit, 並以 switch 來對應不同結果

size_t naive_v2(uint64_t *bitmap, size_t bitmapsize, uint32_t *out)
{
    size_t pos = 0;
    for (size_t k = 0; k < bitmapsize; ++k) {
        uint64_t bitset = bitmap[k];
        uint64_t res;
        size_t p = k * 64;
        for (int i = 0; i < 64; i+=2) {
            switch ((bitset >> i) & 0x3){
            case 3:
                out[pos++] = p + i + 1;
                out[pos++] = p + i;
                break;
            case 2:
                out[pos++] = p + i + 1;
                break;
            case 1:
                out[pos++] = p + i;
                break;
            }
        }  
    }
    return pos;
}

效能比對

naive: 0.001063 (second)
imporved_ctz: 0.0004879 (second)
naive_v2: 0.0004869 (second)

naive: 0.001940 (second)
imporved_ctz: 0.000510 (second)
naive_v2: 0.0005161 (second)

naive: 0.002162534 (second)
imporved_ctz: 0.001259721 (second)
naive_v2: 0.0009134 (second)

naive: 0.002204634 (second)
imporved_ctz: 0.001884666 (second)
naive_v2: 0.001647 (second)

由實驗可知, 在 1 越頻繁出現的圖片中, 以一次看2個 bit 改進的 function 效能較以 __builtin_ctzll 改進的 function 好。

測驗五

LeetCode 166. Fraction to Recurring Decimal
此演算法的步驟是以字串的方式來存數值, 首先先處理整數的部份再處理小數部份

  1. 整數部份:
    將 numerator 除以 denominator, 取計算結果的整數並以 sprintf 將整數部份存成字串, 算出整數結果後要先判斷結果是正或負來決定是否加上負號, 以 numerator mod denominator 計算餘數, 如果餘數不為 0 要加上小數點。
  2. 小數部份:
    以 numerator mod denominator 計算餘數, 並以餘數去計算小數部份, 步驟是先去 hash table 找是否有相同的餘數, 如我存在代表接下來的小數不分將會循環, 所以依照題目將循環部份以 ( 循環部份 ) 的方式儲存並回傳, 如果 hash table 中沒有將餘數存進 hash table 中, 接著將餘數乘以 10 再除以 denominator 存於小數部份, 新的餘數則變為原本的餘數乘以 10 mod denominator。

PPP = pos - -
因為這邊得到的 pos 是開始重複的位置, 所以要先將 decimal 中, 再重複之前的部份接在 p 上
MMM = list_add_tail
將當前的餘數存於 hash table 中
EEE = &heads[remainder % size]
先經過 hash function 得到 index, 將 node 接在 hash table 的 heads[index] 之後

#include <stdbool.h>                       
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "list.h"
    
struct rem_node {
    int key;
    int index;
    struct list_head link;
};  
    
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;
}

char *fractionToDecimal(int numerator, int denominator)
{
    int size = 1024;
    char *result = malloc(size);
    char *p = result;

    if (denominator == 0) {
        result[0] = '\0';
        return result;
    }

    if (numerator == 0) {
        result[0] = '0';
        result[1] = '\0';
        return result;
    }

    /* using long long type make sure there has no integer overflow */
    long long n = numerator;
    long long d = denominator;

    /* deal with negtive cases */
    if (n < 0)
        n = -n;
    if (d < 0)
        d = -d;

    bool sign = (float) numerator / denominator >= 0;
    if (!sign)
        *p++ = '-';

    long long remainder = n % d;
    long long division = n / d;

    sprintf(p, "%ld", division > 0 ? (long) division : (long) -division);
    if (remainder == 0)
        return result;

    p = result + strlen(result);
    *p++ = '.';

    /* Using a map to record all of reminders and their position.
     * if the reminder appeared before, which means the repeated loop begin,
     */
    char *decimal = malloc(size);
    memset(decimal, 0, size);
    char *q = decimal;

    size = 1333;
    struct list_head *heads = malloc(size * sizeof(*heads));
    for (int i = 0; i < size; i++)
        INIT_LIST_HEAD(&heads[i]);
    for (int i = 0; remainder; i++) {
        int pos = find(heads, size, remainder);
        if (pos >= 0) {
            while (pos-- > 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;

        list_add_tail(&node->link, &heads[remainder % size]);

        *q++ = (remainder * 10) / d + '0';
        remainder = (remainder * 10) % d;
    }

    strcpy(p, decimal);
    return result;
}

程式碼改進

sign 判斷, 除法改為 bitwise operation

bool sign = (float) numerator / denominator >= 0;
if (!sign)
    *p++ = '-';

其實在 sign 的判斷部份只需要 numerator > 0 ^ denominator > 0, 可以避免用除法就該避免, 因為除法對於硬體的開銷較大

if (numerator > 0 ^ denominator > 0)
    *p++ = '-';

無須判斷 division 正負號

sprintf(p, "%ld", division > 0 ? (long) division : (long) -division);

因為在上面已經把 n 和 d 都轉為正數, 因此 division 必為正數, 不須判斷正負號, 直接以改為 %lld, 無須再額外轉型。

sprintf(p, "%lld", division);

strcpy -> strncpy

strcpy(p, decimal);

strcpy 有濳在危險, 可能發生 overflow, 覆蓋掉其他資料, 因此改為 strncpy

strncpy(p, decimal, strlen(decimal));
result[strlen(result)] = '\0';

測驗六

ALIGNOF 是用來算某個型態在 struct 中會以幾個 byte 來對齊。

((size_t)&((TYPE *)0)->MEMBER) 的原理是先將 TYPE 型態結構體的地址變更為 0,然後再加上成員 MEMBER 的偏移量, ALIGNOF 是先算出 type t 的偏移量, 因為成員 char 會用 type t 的 alignment 來對齊, 所以算 type t 的偏移量 - 0 等同於算 type t 的 alignment。

至於這邊使用 char 而非其他 type 的原因是 struct 會以成員中 alignemnt 最大的來做對齊, 所以才選 alignemnt 最小的 char。

M = _h
X = 0

/*
 * 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)->_h) - (char *)0

在 Linux 核心原始程式碼中找出 alignof 的使用案例 2 則,並針對其場景進行解說

1. linux/include/linux/kstrtox.h

linux/include/linux/kstrtox.h 中, kstrtol 這個 function 是用來將字串轉為 long, 這邊利用 sizeof 加上 __alignof__ 來檢查 data model, 檢查 long 和 long long 是否皆為 64 bit, 如果 long 和 long long 皆為 64 bit, 則將字串轉成 long long, 否則轉成 long。

為什麼除了 sizeof 檢查以外還需要 __alignof__ 檢查, 目前查到的可能原因為部份硬體架構中, sizeof 的結果並不等於 __alignof__, 因此需要額外檢查 __alignof__, 才能安全的把 long 視為 long long。
例如: 在AMD64架構上,物件是採用自然對齊,多大的物件就對齊多大的記憶體,但在i386架構上,int 佔 4 byte,且為 4-byte aligned,long 也佔 4 bytes,亦為 4-byte aligned。但 long long 和 double 都佔了 8bytes,卻也為 4-byte aligned。

Reference: 關於記憶體對齊(Alignment)

static inline int __must_check kstrtol(const char *s, unsigned int base, long *res)
{
	/*
	 * We want to shortcut function call, but
	 * __builtin_types_compatible_p(long, long long) = 0.
	 */
	if (sizeof(long) == sizeof(long long) &&
	    __alignof__(long) == __alignof__(long long))
		return kstrtoll(s, base, (long long *)res);
	else
		return _kstrtol(s, base, res);
}

在 Linux 核心源使程式碼找出 ALIGN, ALIGN_DOWN, ALIGN_UP 等巨集,探討其實作機制和用途,並舉例探討

#define ALIGN_UP(x, align_to)	(((x) + ((align_to)-1)) & ~((align_to)-1))
#define ALIGN_DOWN(x, align_to) ((x) & ~((align_to)-1))

align_to 須為 power of two, ALIGN_UPALIGN_DOWN 用來得到一個以 x 為上下界且是 align_to 的 m 倍的值, ALIGN_UP 表示 x 是下界, ALIGN_DOWN 則表示 x 是上界。

ALIGN_UP 的實做機制是 align_to

2n, align_to -1 + x 做完 & ~(align_to -1) 後 0 ~ (n - 1) 個 bit 都會為 0, 所以會得到一個 >= x 且為 align_to * m 的數值, 所以ALIGN_UP 回傳的結果會是一個 >= x 且最接近 x 的 align_to * m。

ALIGN_DOWN 的實做機制則是 align_to

2n 的數值, x 做完 & ~(align_to -1) 後 0 ~ (n - 1) 個 bit 都會為 0, 所以會得到一個 <= x 且為 align_to * m 的數值, 所以ALIGN_DOWN 回傳的結果會是一個 <= x 且最接近 x 的 align_to * m。

x 如果為 0 則結果為 0。

測驗七

LeetCode 412. Fizz Buzz

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

        char fmt[9];
        strncpy(fmt, &"FizzBuzz%u"[(8 >> div5) >> (div3 << 2)], length);
        fmt[length] = '\0';

        printf(fmt, i);
        printf("\n");
    }
    return 0;
}

KK1 = div3
KK2 = div5
KK3 = div3 << 2

觀察 div3, div5 與 length

  • div3 = 0, div5 = 0, length = 2 (2 << 0)
  • div3 = 1, div5 = 0, length = 4 (2 << 1)
  • div3 = 0, div5 = 1, length = 4 (2 << 1)
  • div3 = 1, div5 = 1, length = 8 (2 << 2) == (2 << 1 << 1)

由觀察可得答案為 (2 << div3) << div5, div3 跟 div5 位置可互換

接著我們觀察起始位置與 div3, div5 之間的關係

  • div3 = 0, div5 = 0, 起始位置 = 8
  • div3 = 1, div5 = 0, 起始位置 = 0
  • div3 = 0, div5 = 1, 起始位置 = 4
  • div3 = 1, div5 = 1, 起始位置 = 0

如果 div3 為 1 時要右移 4 位, 因此 KK3 為 div3 << 2, 由於 div3 和 div5 為 0 時, 起始位置應為 8, 所以題目原本的 9 應該是錯誤的, 當 div3 和 div5 都為 1 時雖然會右移 5 位, 但 8 >> 4 或 8 >> 5 都為 0, 所以不影響。

評估 naive.c 和 bitwise.c 效能落差

以 perf 來分析 naive.c 和 bitwise.c 的效能

perf stat -r 10 ./xxx
  • naive: 0.002719 sec
  • bitwise: 0.0011213 sec

由實驗結果可見, naive 的執行時間約為 bitwise 的 2.4 倍

分析 Faster remainders when the divisor is a constant: beating compilers and libdivide 一文的想法, 設計另一種 bitmask,然後再改寫 bitwise.c 程式碼,試圖運用更少的指令來實作出 branchless

首先將 bitmask 改為用來判斷是否可除以 15 的, 接著利用 fastmod 的實做手法將 除以 15 後的餘數算出, 接著利用餘數來判斷, 如果餘數 mod 3 == 0 代表可被 3 整除並將最後一個 bit 設為 1, 如果 mod 5 == 0, 代表可被 5 整除並將倒數第二個 bit 設為 1, 如果 mod 3 和 5 都為 0 就是可被 15 整除, 這邊使用到 !! 技巧, 不管餘數是多少只會得到 0 or 1, 得到回傳結果後 div3 = result & 1, 而 div5 = result >> 1, 後面的實做和原本 bitwise.c 皆相同。

實驗結果

改進後效能比原先的 bitwise.c 來得好, 但指令的部份比較多

  • naive: 0.002719 sec
  • bitwise: 0.0011213 sec
  • bitwise_v2: 0.0008026 sec
include <stdio.h>
#include <stdbool.h>
#include <stdint.h>
#include <string.h>

uint64_t M15 = UINT64_C(0xFFFFFFFFFFFFFFFF) / 15 + 1;

uint8_t divisible(uint32_t n) {
    uint64_t lowbits = M15 * n;
    uint32_t remain = ((__uint128_t)lowbits * 15) >> 64; 
    uint8_t mask = 0;
    mask |= !!(remain % 3);
    mask |= !!(remain % 5 == 0) << 1;
    return mask;
}

int main(int argc, char **argv)
{
    for (size_t i = 1; i <= 100; i++) {
        uint32_t div = divisible(i);
        uint8_t div3 = div & 1;
        uint8_t div5 = div >> 1;
        unsigned int length = (2 << div3) << div5;
        char fmt[9];
        strncpy(fmt, &"FizzBuzz%u"[(8 >> div5) >> (div3 << 2)], length);
        fmt[length] = '\0';
        printf(fmt, i);
        printf("\n");
    }
    return 0;
}

TODO:
解釋 fastmod

reference: