# 2020q3 Homework10 (quiz 10)
contributed by < `JimmyLiu0530` >

###### tags: `進階電腦系統理論與實作`

## 測驗1: RGB to Gray
給定每個像素 (pixel) 為 32-bit 的 RGBA 的位元圖 (bitmap)，RGBA 代表紅綠藍三原色的首字母，Alpha 值則表示顏色的透明度/不透明度，其轉換為灰階影像 (gray scale) 的函式為：

void rgba_to_bw(uint32_t *bitmap, int width, int height, long stride) {
    int row, col;
    uint32_t pixel, r, g, b, a, bw;
    
    for (row = 0; row < height; row++) {
        for (col = 0; col < width; col++) {
            pixel = bitmap[col + row * stride / 4];
            
            a = (pixel >> 24) & 0xff;
            r = (pixel >> 16) & 0xff;
            g = (pixel >> 8) & 0xff;
            b = pixel & 0xff;
            bw = (uint32_t) (r * 0.299 + g * 0.587 + b * 0.114);
            
            bitmap[col + row * stride / 4] = (a << 24) + (bw << 16) + (bw << 8) + (bw);
        }
    }
}

人類視網膜包含三類錐形對光敏感的細胞，分別對不同波長範圍的光有反應:
> 上圖中的縱軸代表人眼對某對應波長的光的敏感度。其中 1 Å (埃) = $10^{-10}$米 = $0.1nm$ [[出處](http://www.phy.ntnu.edu.tw/demolab/html.php?html=everydayPhysics/color)]

人眼吸收綠色比其他顏色敏感，也可說人眼最容易捕捉到綠色，所以當影像變成灰階時，僅僅將紅色、綠色、藍色加總取平均，不足以反映出人眼所見。常見的方法是將 $Red×77,Green×151,Blue×28$，這三個除數的總和為 `256` (即 $2^8$)，可使除法變簡單 (等同於右移 8 位元)。

如果我們將所有的運算都建立表格，於是

bwPixel = table[rgbPixel & 0x00ffffff] + rgbPixel & 0xff000000;

$→$ 16 MB; 表格太大

如果先計算針對「乘上 0.299」一類的運算，計算後再建立表格呢？

bw = (uint32_t) mul_299[r] + (uint32_t) mul_587[g] + (uint32_t) mul_144[b];
bitmap[col + row * strike / 4] = (a << 24) + (bw << 16) + (bw << 8) + bw;

$→$ 降到 32 KB 以內; cache friendly

接著嘗試使用 NEON 指令集加速。將 RGB 三色的權重存入 r3, r4, r5 等暫存器。

`vdup.8` (Vector Duplicate)，分別複製到大小為 8 bit 的 NEON register d0 - d2

    mov     r3, #77
    mov     r4, #151
    mov     r5, #28
    vdup.8  d0, r3
    vdup.8  d1, r4
    vdup.8  d2, r5

`vld4.8` (Vector Load)，載入 pixel 的資料到 4 個 8-bit 的 NEON register d4-d7，其中那個 4 為 interleave，因為我們有 ARGB，所以 gap = 4。

再來就是計算 weighted average。可用 Vector Multiply 和 Vector Multiply Accumulate 指令:

   @ (alpha,R,G,B) = (d7,d6,d5,d4)
   vld4.8      {d4-d7}, [r0]!
   vmull.u8    q10, d6, d0
   vmlal.u8    q10, d5, d1
   vmlal.u8    q10, d4, d2

將值除以 256 就是我們要的灰階值。

`vrshrn` (Vector Shift Right by immediate value)
> vrshrn.u16 d4, q10, #8

最後儲存結果。

`vst` (Vector Store)
> vst4.8 {d4-d7}, [r3]!

效能評比如下: (對照 `v0` 和 `v4`)

從這張表即可理解，使用 NEON 指令加速只用原本 4.6％ 的時間即可完成 RGBA 轉灰階的處理。

[Neon Intrinsics](https://developer.arm.com/architectures/instruction-sets/simd-isas/neon/intrinsics) 允許我們避免直接撰寫組合語言，而是有如 C 函式呼叫，例如 `vmull.u8` 和 `vmlal.u8` 分別對應以下 intrinsics: (提示: 在 Neon Intrinsics 網頁搜尋 NEON intrinsics 可獲得更詳盡的描述)
- [vmull_u8](https://developer.arm.com/architectures/instruction-sets/simd-isas/neon/intrinsics?search=vmull_u8)
- [vmlal_u8](https://developer.arm.com/architectures/instruction-sets/simd-isas/neon/intrinsics?search=vmlal_u8)

考慮以下使用 `float` 來保存像素的程式碼:

#include <stdint.h>
#define MIN(a, b) (((a) < (b)) ? (a) : (b))
#define MAX(a, b) (((a) > (b)) ? (a) : (b))

uint8_t *to_rgb(const float *m, int w, int h, uint8_t *rgb)
{
    int size = w * h;
    const float *ptr0 = m + (0 * size);
    const float *ptr1 = m + (1 * size);
    const float *ptr2 = m + (2 * size);
    uint8_t *ret = rgb;

#define SATURATE_CAST_UCHAR(X) (uint8_t) MIN(MAX((int) (X), 0), 255);
    for (int remain = size; remain > 0; remain--) {
        rgb[0] = SATURATE_CAST_UCHAR(*ptr0);
        rgb[1] = SATURATE_CAST_UCHAR(*ptr1);
        rgb[2] = SATURATE_CAST_UCHAR(*ptr2);
        rgb += 3;
        ptr0++, ptr1++, ptr2++;
    }
#undef SATURATE_CAST_UCHAR

    return ret;
}

函式 `to_rgb` 可將原本透過連續的 `float` 空間保存的資料轉為 RGB 像素，每個像素佔用 3 個位元組，即 24 位元。針對 Aarch64 (Arm 64-bit) 架構 (假設為 little-endian)，下方的程式碼則可將 RGB 轉換為灰階:

#include <arm_neon.h>
float *rgb_to_gray(const uint8_t *rgb, int w, int h, float *m)
{
    const uint8_t Y_shift = 8;
    const uint8_t R2Y = 77, G2Y = 151, B2Y = 28;
    float *ptr = m;
    int size = w * h;

    int nn = size >> 3;
    int remain = size - (nn << 3);

    uint8x8_t _R2Y = vdup_n_u8(R2Y);
    uint8x8_t _G2Y = vdup_n_u8(G2Y);
    uint8x8_t _B2Y = vdup_n_u8(B2Y);

    for (; nn > 0; nn--) {
        uint8x8x3_t _rgb = vld3_u8(rgb);

        uint16x8_t y16 = vmull_u8(_rgb.val[0], _R2Y);
        y16 = vmlal_u8(y16, _rgb.val[1], _G2Y);
        y16 = vmlal_u8(y16, _rgb.val[2], _B2Y);
        y16 = vshrq_n_u16(y16, Y_shift);

        float32x4_t _ylow = vcvtq_f32_u32(vmovl_u16(vget_low_u16(y16)));
        float32x4_t _yhigh = vcvtq_f32_u32(vmovl_u16(vget_high_u16(y16)));

        vst1q_f32(ptr + OFFSET1, _ylow);
        vst1q_f32(ptr + OFFSET2, _yhigh);

        rgb += 3 * 8, ptr += 8;
    }

    for (; remain > 0; remain--) {
        *ptr = (rgb[0] * R2Y + rgb[1] * G2Y + rgb[2] * B2Y) >> Y_shift;
        rgb += 3, ptr++;
    }

    return m;
}

請補完程式碼，使上述程式碼得以符合預期地運作。

### 程式說明
line 4-7 初始化等等會用到的變數，像是 RGB 的權重、位移量...等。

因為本題是在 64-bit 的環境下執行，因此暫存器也是 64-bit，代表一輪可以處理 8 個像素。因此， `nn` 紀錄需要做的輪數，剩餘不足一輪的會在 line 33-35 單獨處理。

line 12-14 將 RGB 各顏色的權重分別載入到 64-bit 暫存器，其中載入的方法為複製 8 份 1-byte 權重(e.g. `R2Y`、`G2Y`、`B2Y`)到暫存器中。

每輪執行以下步驟:
1. 先將 8 個像素一次載入到三個暫存器。

    uint8x8x3_t _rgb = vld3_u8(rgb);

載入的方式是採 `interleaving`，也就是說從 `rgb` 取到的第一個 byte 會載入到第一個暫存器，取到的第二個 byte 會載入到第二個暫存器，取到的第三個 byte 會載入到第三個暫存器，取到的第四個 byte 會載入到第一個暫存器。

2. 計算出灰階值

    uint16x8_t y16 = vmull_u8(_rgb.val[0], _R2Y);
    y16 = vmlal_u8(y16, _rgb.val[1], _G2Y);
    y16 = vmlal_u8(y16, _rgb.val[2], _B2Y);
    y16 = vshrq_n_u16(y16, Y_shift);

`vmull_u8` 會將第一個 argument 與第二個 argument 的各個 byte 相乘，因此乘出來的結果就會變成 16 bits。全部的結果就會變成 16 x 8 = 128 bits。`vmulalu8` 就只是 `vmull_u8` 乘完再與另一個暫存器相加而已。
最後，將計算出來的加權值除以 `256`，也就是向右位移 `8` 後，即為灰階值。

3. 存回結果

    float32x4_t _ylow = vcvtq_f32_u32(vmovl_u16(vget_low_u16(y16)));
    float32x4_t _yhigh = vcvtq_f32_u32(vmovl_u16(vget_high_u16(y16)));

    vst1q_f32(ptr + OFFSET1, _ylow);
    vst1q_f32(ptr + OFFSET2, _yhigh);

    rgb += 3 * 8, ptr += 8;

現在要將每個像素的 RGB 從原本 16-bit unsigned 轉成 32-bit float，因為我們原本就是使用 `float` 來保存資料。
轉換的方法如下:
- 將原本存於 `y16` 中 8 個像素的資料分成兩群各自去做轉換，因此 `_ylow`、`_yhigh` 裏頭各自有 4 個像素的資料。
- 最後將 `_ylow`、`_yhigh` 存回 `ptr` 中。顯然，`_ylow` 要存回的初始位置就是 `ptr`，所以`OFFSET1 = 0`。因為 `ylow` 佔 32*4 bits，也就是 4 個 float，所以 `_yhigh` 要從 `ptr` + 4 開始存，因此 `OFFSET1 = 4`。

### 延伸問題
#### 1. 整合 [libattopng](https://github.com/misc0110/libattopng)，輸出 RGBA 和灰階的 PNG 圖片
> - 可能會用到 [ImageMagick](https://imagemagick.org/) 事先處理輸入檔案
> - 參照 [Bomb Lab](https://hackmd.io/@posutsai/SyNzl72f4) 實作紀錄 得知如何安裝 GNU Toolchain 及 QEMU

----------------------------------------------

## 測驗2: LeetCode [1617. Count Subtrees With Max Distance Between Cities](https://leetcode.com/problems/count-subtrees-with-max-distance-between-cities/)

給定 $n$ 個城市，編號從 $1$ 到 $n$，再給予一個大小為 $n-1$ 的陣列 `edges`，其中 $edges[i]=[u_i,v_i]$ 表示城市 $u_i$ 和 $v_i$ 之間存在雙向邊。題目保證任兩個城市間有唯一的路徑，也就是說所有城市構成一棵樹 (tree)。


> 上圖的聯通關係為:
> - 1, 4
> - 1, 3
> - 3, 2
> - 1, 3, 2
> - 4, 1, 3
> - 4, 1, 3, 2

一棵子樹 (subtree) 是上述所有城市的子集。在子集中，任意城市間可透過子集的其他城市和邊到達。兩個子樹被認為不一樣的條件是，至少有個城市在其中一棵子樹中存在，卻又不存在於另一棵子樹中。

給定 $d$ 範圍介於 $1$ 到 $n-1$，請找出城市間最大距離恰好為 $d$ 的所有子樹數目。請你回傳一個大小為 $n-1$ 的陣列，其中第 $d$ 個元素（下標從 1 開始）是城市間最大距離恰好等於 $d$ 子樹數目。

> 兩個城市間距離定義是，它們之間需要經過的邊的數目。


第一組輸入:
- $n=4$
- $edges=[[1,2],[2,3],[2,4]]$

預期輸出：
- $[3,4,0]$

> 既然 $n=4$，於是輸出的陣列會有 3 個元素
> - 子樹 `{1,2}`, `{2,3}` 和 `{2,4}` 最大距離都為 `1`
> - 子樹 `{1,2,3}`, `{1,2,4}`, `{2,3,4}` 和 `{1,2,3,4}` 最大距離都為 `2`
> - 不存在城市間最大距離為 `3` 的子樹。

第二組輸入:
- $n=2$
- $edges=[[1,2]]$

預期輸出：
- $[1]$

> 既然 $n=2$，於是輸出的陣列只會有 1 個元素，且題目已保證任意城市之間只有唯一的路徑

第三組輸入:
- $n=3$
- $edges=[[1,2],[2,3]$

預期輸出：
- $[2,1]$

> 既然 $n=3$，於是輸出的陣列會有 2 個元素
> - 子樹 `{1,2}` 和 `{2,3}` 最大距離都為 `1`
> - 子樹 `{1,2,3}` 最大距離為 `2`

限制:
- $2 \leq n \leq 15$
- $edges.length==n-1$
- $edges[i].length==2$
- $1 \leq u_i,v_i \leq n$
- 題目保證 $(u_i,v_i)$ 所表示的邊互不相同

思路:
- 因為 $n \leq 15$，所以就往狀態壓縮 + dynamic programming 的技巧去思考，但本題狀態間的轉移無法進行，只能窮舉所有情況。
- 題目的樹可看成是無向無環連通圖（二者等價）
> 在圖論中，環是一條只有第一個和最後一個頂點重複的非空路徑
> - 不具備環的圖被稱作無環圖
> - 不具備有向環的有向圖被稱做有向無環圖 (DAG,directed acyclic graph)
> - 無環的連通圖被稱作樹
> 

- 窮舉狀態，首先要判斷是否無效
1. 該狀態僅有 $1$ 個節點，因題目要求的 $d$ 最小是 $1$，即距離為 $1$，所以可縮減
2. 狀態下不是一個子樹，即非一個連通圖，那怎麼判斷是否是連通圖呢？可利用 BFS 或 DFS 判斷所能到達的節點的個數，是否和集合內元素個數 (二進位 $1$ 的個數)相同即可

- 判斷該狀態是有效後，如何求最大距離？所謂的最大距離是無向無環連通圖的最大直徑

以下程式碼是可能的解法:

#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#define MIN(a, b) (((a) < (b)) ? (a) : (b))

#define ITERATE(i, begin, end) \
    for (int i = (begin), i##_end_ = (end); i < i##_end_; i++)

int *countSubgraphsForEachDiameter(int n,
                                   int **edges,
                                   int edgesSz,
                                   int *edges0sz,
                                   int *rsz)
{
    uint8_t d[n][n]; /* distance */
    memset(d, 0xFF, sizeof(d));
    uint16_t conn[n];
    ITERATE (i, 0, n)
        conn[i] = 1 << i;

    for (int z = 0; z < edgesSz; z++) {
        int *e = edges[z];
        int u = e[0] - 1, v = e[1] - 1;
        d[u][v] = d[v][u] = 1;
        conn[u] |= 1 << v, conn[v] |= 1 << u;
    }

    ITERATE (k, 0, n)
        ITERATE (i, 0, n)
            ITERATE (j, 0, n)
                d[i][j] = MIN(d[i][j], d[i][k] + d[k][j]);
    int *rv = calloc(n - 1, sizeof(int));

    ITERATE (S, 1, (1 << n)) {
        int ctmp = 1 << __builtin_ctz(S);
        while (1) {
            int conn_nxt = ctmp;
            ITERATE (d, 0, n)
                if ((ctmp >> d) & 1)
                    conn_nxt |= conn[d];
            conn_nxt &= S;
            if (COND1)
                break;
            ctmp = conn_nxt;
        }
        if (COND2)
            continue;
        uint8_t ret = 0;
        ITERATE (i, 0, n)
            if ((S >> i) & 1)
                ITERATE (j, 0, i)
                    if (((S >> j) & 1) && (ret < d[i][j]))
                        ret = d[i][j];
        if (!ret)
            continue;
        rv[ITER]++;
    }

    *rsz = n - 1;
    return rv;
}

參考資訊:
- [花花醬解說: Count Subtrees With Max Distance Between Cities](https://www.youtube.com/watch?v=qag7m9VRf-A&feature=youtu.be)
- [Bit manipulation and diameter of a tree](https://jimmy-shen.medium.com/bit-manipulation-and-diameter-of-a-tree-3856d549c5f8)

### 程式說明
首先，介紹一些變數的意義:
- `d[i][j]`:節點 `i` 與節點 `j` 之間的最短距離
- `conn[i]`: 與節點 `i` 相鄰的點，用二進位表示

根據 input，將 `conn`、`d` 做初始化，然後計算兩點間的最短距離。

上面提到，我們採用窮舉所有子樹來求解，不過我們要先去檢查該子樹是否為連通? 作法說明如下:

`S` 代表一顆子樹，從這子樹中 index 最小的點開始，不斷新增與此點相鄰的點進 `conn_nxt`，因為我們只在乎這顆子樹的點而已，所以其他不在這子樹的點我們都忽略，對應到程式碼

    conn_nxt &= S;

透過不斷的從現有的點中往外擴展，直到收斂，因此，`COND1 為 conn_nxt == ctmp`。
此時，若該子樹為連通，則 `ctmp` 一定包含該子樹所有的點。反之，如果
`ctmp` 沒有包含該子樹所有的點，也就是說 `~ctmp & S != 0`，則該子樹不連通，因此 `COND2 == ~ctmp & S`。

接著確認是連通圖後，針對該子樹去找出這棵樹的 diameter，即 $max_{u,v∈V} \ d[u][v]$，其中 $V$ 為該子樹的點集合。找法說明如下:

看過所有該子圖中的點與子圖中其他點的距離，就可以找出最大的距離，即 diameter。
最後，根據該子樹的 diameter 去遞增相對應 `rv[]` 的位置，因此，`ITER == --ret`。

- [ ] 指出改進方式
- 更有效率找 all pair shortest path 的方法-[BFS](https://mathoverflow.net/questions/59680/all-pairs-shortest-paths-in-trees): $O(n^3)$ to $O(n^2)$

### 延伸問題
#### 1. 用 C99/C11 (可搭配 GNU extension) 重新實作符合 LeetCode 要求但效率優於上述表現的實作，一併分析時間複雜度。