Skia’s New Approach to SIMD
===
延續之前針對[「浮點數的運算」](/s/H1SVhbETQ#)比較定點數與浮點數的效能差異，根據 jserv 老師的建議去了解 Skia 在 SIMD 上的應用，並且嘗試改用 SSE/AVX 這類 SIMD 指令進行實驗

### 什麼是 Skia ？ 跟 SIMD 又是什麼關係 ？
Skia 本身是 Google 家的開源專案， 主要是跨平台的 2D 圖像函式庫，因為相對其他函式庫有著極大的效能最佳化，所以被各大應用採用，像是 Google Chrome, Android, Mozilla Firefox 等等

關於 Skia 怎麼用請看這篇共筆 ([十分鐘讓你用 Skia 畫一張圖](/s/SyCFz7Kb4)) ，這邊著重討論以 SIMD 實做 Skia 的部份

不過如果你現在嘗試去找 Skia + SIMD 的[官方文件](https://github.com/google/skia/blob/ce5831e004e1feba48d4b4fbc70cb7c73ef0edec/site/dev/contrib/simd.md) (此文件是過去 github 留存的舊版本，2018/12/17 之後這個文件就被砍掉了)， 你會獲得一個精美的 Error 404


這主要是因為就在上個禮拜 (2018/12/17)，這份文件因為 " simd .md describes a whole system that doesn't exist" 所以被砍掉了，這讓我非常緊張，一來是我這樣到底要怎麼研究這個實做 （萬一 Skia 就只是「恩理論上我覺得應該可以用到 SIMD !」 但是實務上做不出來直接砍掉） ，二來就是原來 Google 這樣的大公司也會這樣畫大餅然後裝死


根據 Skia's New Approach to SIMD 的說法，要建立一個跨平台通用的 C++ 程式， x86 架構上的 SSE 家族指令/ ARM 架構下的 NEON 指令/ MIPS32 架構下的 DSP 指令可以達到最佳效能表現，可到達甚至 16 被的效能差異

不過這些 SIMD 底層指令通常都相當晦澀難懂（光是要讀懂就要花費一段時間了），因此在 Skia 上 SIMD 指令覆蓋率並不完全，此後講解的方向可能偏向以 Arm NEON 或 MIPS32 DSP 實現； 同時即使指令覆蓋率足夠，該實現的正確性也令人擔憂，就算知道程式本身可能存在 bug， 但是由於 SIMD 指令的難以閱讀因此修復 bug 也便得相當困難，使得某些溢位錯誤長時間已知但找不出原因； 最後是，SIMD 指令並不能保證 「總是」 可以獲得效能增益

:::info
那到底什麼是 SIMD 指令呢？
:::

SIMD (Single Instruction Multiple Data)，表示對同一種數據型態進行相同指令操作，屬於指令集層級的平行化運算，像是 Intel x86 架構下的 SSE

那使用了這些 SIMD 有什麼優點呢？ 要怎麼使用這些 SIMD 指令呢？

1. 簡單版 → 交給 GCC 的強大最佳化編譯選項 -O3

### 架構規格 （畫個餅因為還沒做出來）
整理來說， Skia 用到資料型態主要為 32位元/ 16位元/ 8位元整數以及浮點數運算， 8位元主要出現在 ARGB 像素資料上，可以說得上是最常用的 SIMD 對象，但是由於運算中常用到乘法運算可能會導致溢位，所以需要使用 16位元整數暫存，同時，對於一些特別需要高精度運算的指令，我們也需要提供 「8 位元整數轉換成浮點數」的指令；最後因為 32位元資料型態不常使用到因此並不在這次優化目標中

也就是說，這個架構主要目標在於以 「建立一個跨架構平台（EX. x86, ARM ......）的上層抽象 API ，透過這層 API 包裝成對應平台的 SIMD 指令，作為處理浮點數以及 8/16 位元整數運算指令界面」

值得注意的是，**x86 SSE 與 ARM NEON 都使用 128-bits vector 作為輸入**，正好可以塞進去 4 個 32-bits float 做 SIMD 運算；同時兩者對於**類似運算幾乎都可以找到互為對應的函式**，像是 NEON vmulq_f32() 對應到 SSE \_mm_mul_ps()

這樣一來就可以透過簡單的 #define 把不同架構上的指令集根據本地架構設定轉換成對應的指令集；不過需要注意到幾個問題

1. 部份指令不存在對應關係，像是 \_mm_movemask_ps()
2. 以 SSE 或 NEON 指令寫成的數學運算還是很難懂
3. 部份指令的參數數量可能會隨著使用場景需要調整（有時候我們希望是 2 個 float，有時後是 4 個等等）

因此 `SkNf <N>` 這個 wrapper class 需要可以透過後面的 N 參數建立一個「有 N 個浮點數的向量」，同時這個向量物件需要能夠提供對應參數量的運算函式，像是求最小/大，倒數，平方根，比較大小等

同時為了 Skia 在圖像上的支援， SkNf 也必須要提供「讀取 N byte 值作為 float 讀入」並且「無條件捨去該浮點數成 0 ~ 255 存入」的功能

因此為了支援不同架構支援不同位元數的浮點數運算
1. SkNf<1> 是最基本的運算單元，代表只有一個浮點數的向量單元；所有大於 1 個浮點數的向量都會以遞迴的方式對半切成兩個小向量，直到剩下的向量都是 SkNf <1>
2. 根據不同指令集有不同加速實做；以 ARM NEON 為例， SkNf <2> 可以透過 float32x2_t 與 64 位元的 SIMD 運算實做；而 SkNf <4> 則可以透過 float32x4_t 加上 128 位元 SIMD 實做
3. x86 SSE 則可以透過 \_\_m128 取全部位元實做 SkNf <4> 或取一半位元實做 SkNf <2>
4. 由於 「二維座標」「以四個浮點數表示的像素值」 這兩個使用情境是最常出現的，因此分別提供對應的物件 Sk2f， Sk4f；同時還有其對應到的向量模式 Sk2s，Sk4s

### SkNf

struct SkNx {
    typedef SkNx<N/2, T> Half;

    Half fLo, fHi;

    AI SkNx() = default;
    AI SkNx(const Half& lo, const Half& hi) : fLo(lo), fHi(hi) {}

    AI SkNx(T v) : fLo(v), fHi(v) {}

    AI SkNx(T a, T b)           : fLo(a)  , fHi(b)   { static_assert(N==2, ""); }
    AI SkNx(T a, T b, T c, T d) : fLo(a,b), fHi(c,d) { static_assert(N==4, ""); }
    AI SkNx(T a, T b, T c, T d,  T e, T f, T g, T h) : fLo(a,b,c,d), fHi(e,f,g,h) {
        static_assert(N==8, "");
    }
    AI SkNx(T a, T b, T c, T d,  T e, T f, T g, T h,
            T i, T j, T k, T l,  T m, T n, T o, T p)
        : fLo(a,b,c,d, e,f,g,h), fHi(i,j,k,l, m,n,o,p) {
        static_assert(N==16, "");
    }

    AI T operator[](int k) const {
        SkASSERT(0 <= k && k < N);
        return k < N/2 ? fLo[k] : fHi[k-N/2];
    }

    AI static SkNx Load(const void* vptr) {
        auto ptr = (const char*)vptr;
        return { Half::Load(ptr), Half::Load(ptr + N/2*sizeof(T)) };
    }
    AI void store(void* vptr) const {
        auto ptr = (char*)vptr;
        fLo.store(ptr);
        fHi.store(ptr + N/2*sizeof(T));
    }

    AI static void Load4(const void* vptr, SkNx* a, SkNx* b, SkNx* c, SkNx* d) {
        auto ptr = (const char*)vptr;
        Half al, bl, cl, dl,
             ah, bh, ch, dh;
        Half::Load4(ptr                  , &al, &bl, &cl, &dl);
        Half::Load4(ptr + 4*N/2*sizeof(T), &ah, &bh, &ch, &dh);
        *a = SkNx{al, ah};
        *b = SkNx{bl, bh};
        *c = SkNx{cl, ch};
        *d = SkNx{dl, dh};
    }
    AI static void Load3(const void* vptr, SkNx* a, SkNx* b, SkNx* c) {
        auto ptr = (const char*)vptr;
        Half al, bl, cl,
             ah, bh, ch;
        Half::Load3(ptr                  , &al, &bl, &cl);
        Half::Load3(ptr + 3*N/2*sizeof(T), &ah, &bh, &ch);
        *a = SkNx{al, ah};
        *b = SkNx{bl, bh};
        *c = SkNx{cl, ch};
    }
    AI static void Load2(const void* vptr, SkNx* a, SkNx* b) {
        auto ptr = (const char*)vptr;
        Half al, bl,
             ah, bh;
        Half::Load2(ptr                  , &al, &bl);
        Half::Load2(ptr + 2*N/2*sizeof(T), &ah, &bh);
        *a = SkNx{al, ah};
        *b = SkNx{bl, bh};
    }
    AI static void Store4(void* vptr, const SkNx& a, const SkNx& b, const SkNx& c, const SkNx& d) {
        auto ptr = (char*)vptr;
        Half::Store4(ptr,                   a.fLo, b.fLo, c.fLo, d.fLo);
        Half::Store4(ptr + 4*N/2*sizeof(T), a.fHi, b.fHi, c.fHi, d.fHi);
    }
    AI static void Store3(void* vptr, const SkNx& a, const SkNx& b, const SkNx& c) {
        auto ptr = (char*)vptr;
        Half::Store3(ptr,                   a.fLo, b.fLo, c.fLo);
        Half::Store3(ptr + 3*N/2*sizeof(T), a.fHi, b.fHi, c.fHi);
    }
    AI static void Store2(void* vptr, const SkNx& a, const SkNx& b) {
        auto ptr = (char*)vptr;
        Half::Store2(ptr,                   a.fLo, b.fLo);
        Half::Store2(ptr + 2*N/2*sizeof(T), a.fHi, b.fHi);
    }

    AI T min() const { return SkTMin(fLo.min(), fHi.min()); }
    AI T max() const { return SkTMax(fLo.max(), fHi.max()); }
    AI bool anyTrue() const { return fLo.anyTrue() || fHi.anyTrue(); }
    AI bool allTrue() const { return fLo.allTrue() && fHi.allTrue(); }

    AI SkNx    abs() const { return { fLo.   abs(), fHi.   abs() }; }
    AI SkNx   sqrt() const { return { fLo.  sqrt(), fHi.  sqrt() }; }
    AI SkNx  rsqrt() const { return { fLo. rsqrt(), fHi. rsqrt() }; }
    AI SkNx  floor() const { return { fLo. floor(), fHi. floor() }; }
    AI SkNx invert() const { return { fLo.invert(), fHi.invert() }; }

    AI SkNx operator!() const { return { !fLo, !fHi }; }
    AI SkNx operator-() const { return { -fLo, -fHi }; }
    AI SkNx operator~() const { return { ~fLo, ~fHi }; }

    AI SkNx operator<<(int bits) const { return { fLo << bits, fHi << bits }; }
    AI SkNx operator>>(int bits) const { return { fLo >> bits, fHi >> bits }; }

    AI SkNx operator+(const SkNx& y) const { return { fLo + y.fLo, fHi + y.fHi }; }
    AI SkNx operator-(const SkNx& y) const { return { fLo - y.fLo, fHi - y.fHi }; }
    AI SkNx operator*(const SkNx& y) const { return { fLo * y.fLo, fHi * y.fHi }; }
    AI SkNx operator/(const SkNx& y) const { return { fLo / y.fLo, fHi / y.fHi }; }

    AI SkNx operator&(const SkNx& y) const { return { fLo & y.fLo, fHi & y.fHi }; }
    AI SkNx operator|(const SkNx& y) const { return { fLo | y.fLo, fHi | y.fHi }; }
    AI SkNx operator^(const SkNx& y) const { return { fLo ^ y.fLo, fHi ^ y.fHi }; }

    AI SkNx operator==(const SkNx& y) const { return { fLo == y.fLo, fHi == y.fHi }; }
    AI SkNx operator!=(const SkNx& y) const { return { fLo != y.fLo, fHi != y.fHi }; }
    AI SkNx operator<=(const SkNx& y) const { return { fLo <= y.fLo, fHi <= y.fHi }; }
    AI SkNx operator>=(const SkNx& y) const { return { fLo >= y.fLo, fHi >= y.fHi }; }
    AI SkNx operator< (const SkNx& y) const { return { fLo <  y.fLo, fHi <  y.fHi }; }
    AI SkNx operator> (const SkNx& y) const { return { fLo >  y.fLo, fHi >  y.fHi }; }

    AI SkNx saturatedAdd(const SkNx& y) const {
        return { fLo.saturatedAdd(y.fLo), fHi.saturatedAdd(y.fHi) };
    }

    AI SkNx mulHi(const SkNx& m) const {
        return { fLo.mulHi(m.fLo), fHi.mulHi(m.fHi) };
    }
    AI SkNx thenElse(const SkNx& t, const SkNx& e) const {
        return { fLo.thenElse(t.fLo, e.fLo), fHi.thenElse(t.fHi, e.fHi) };
    }
    AI static SkNx Min(const SkNx& x, const SkNx& y) {
        return { Half::Min(x.fLo, y.fLo), Half::Min(x.fHi, y.fHi) };
    }
    AI static SkNx Max(const SkNx& x, const SkNx& y) {
        return { Half::Max(x.fLo, y.fLo), Half::Max(x.fHi, y.fHi) };
    }
};