RapidJSON 代碼剖析（四）：優化 Grisu

時間 2019-11-06

原文原文鏈接

我曾經在知乎的一個答案裏談及到 V8 引擎裏實現了 Grisu 算法，我先引用該文的內容簡單介紹 Grisu。而後，再談及 RapidJSON 對它作了的幾個底層優化。git

（配圖中的《Grisù》是一套1970年代的意大利卡通短片，主角 Grisù 是一隻想成爲消防員的小龍。估計 Grisu 算法以龍命名，是與上一代的 Dragon4 算法相關。）github

Grisu是什麼

Grisu 是把浮點數轉換爲字符串的算法。在 Chrome 裏執行這段 JavaScript 實際上就調用了 Grisu：算法

document.write(1/3); // 0.3333333333333333

這個問題看似簡單，其實是很複雜的事情。json

在1980年以前，許多 C 語言標準庫中的 printf() 都會產生「不正確」的結果。Coonen 在那時候作了相關的博士研究[1]，但並無受到普遍的知悉和應用。1990年 Steele 等人發表了名爲 Dragon4 的算法[2]，經過使用大數運算來確保精確的轉換。而這個算法應用在大部分 C 程序語言庫的 printf()，以及其餘用 printf() 來實現這功能的語言之中。c#

這樣就20年過去，雖然中間有一些小改進[3][4]，但基本的思路仍然是這樣。到了2010年，Loitsch 發表了 Grisu 算法[5]，而 V8 也實現了這個算法。api

該篇文章闡述了3個算法，分別是原始的 Grisu，及其改進版 Grisu2 和 Grisu3。Grisu 算法比 Dragon4 等算法優越的地方就是一個字──快。如下簡單介紹一下它們的特色。架構

首先，什麼是「正確」的轉換？其定義是，一個浮點數轉換成的十進位字符串以後，該字符串能夠完美的轉換回去原來的浮點數，若是用C語言來描述的話：app

// 除 +/-inf 及 NaN 外的浮點數都應該傳回 true。
bool Verify(double d) {
    assert(!isnan(d) && !isinf(d));
    char buffer[32]; // 足夠大麼？
    dtoa(buffer, d); // 雙精度浮點數轉換至字符串，是double-to-ascii
    char* end;
    double d2 = strtod(buffer, &end); // 字符串轉換至雙精度浮點數
    return *end == '\0' && d == d2;
}

如前所述，Dragon4 使用大數運算，還涉及動態內存分配（你知道printf()裏可能會作malloc()麼？）。而 Grisu 則不須要，只須要用到64位整數運算即可以實現轉換！因此那篇文章題目就以 "... with integers" 做結尾。ide

Grisu 的算法很是簡單，但它有一個缺點，就是其結果並不像是給人看的。如文中的例子，Grisu 把 0.3 的打印成 29999999999999998e-17。這是「正確的」轉換結果，它能夠經過 Verify() 驗證。函數

雖然該算法很是快，但通常人大概不會接受這樣的結果。做者也所以研發出改進版本 Grisu2，在使用64位整數實現 double 的轉換時，能夠利用額外的 64 - 53 = 11 位去縮減轉換的結果（53 爲 double 的尾數位數）。Grisu2 能夠把 >99.9% 的浮點數轉換成最短的「正確」字符串，其餘<0.1%的浮點數仍然是「正確」的，但不是最短的答案。

也許通常人就見好就收了，畢竟已證實算法的正確性，只是有那麼 <0.1% 狀況未達至最完美的結果。不過該做者仍是繼續研究出 Grisu3。Grisu3 並不能解決那一小撮麻煩製造者，但它能在計算期間偵查到哪些 double 在這算法中未能得出最優的結果。既然辦事快捷的小部門搞不定，就能夠把它交給 Dragon4 或其餘較慢的算法。

V8 裏實現了 Grisu3 以及大整數的算法（我不願定是Dragon4仍是其餘），後來Google也把它分離成爲一個獨立的C++程序庫double-conversion。

爲了優化 RapidJSON 的浮點數轉換，也因爲 RapidJSON 是僅需頭文件的 JSON 庫，我就按 Loitsch 的實現編寫了一個 Grisu2 的頭文件庫，能夠在 dtoa-benchmark，那裏還比較了多個 dtoa() 實現的性能。由於 Grisu3 須要另外一個更龐大的大數實現，而暫時 RapidJSON 不須要最優結果，因此就僅實現了一個性能更好及更簡短的 Grisu2。

測試不一樣算法／實現下的性能（VC2013 64-bit），fpconv、grisu2、milo 都是 Grisu2 的實現，doubleconv 是 V8 的 Grisu3 實現。milo 對 sprintf 的加速比約是 9x。

加入了經優化的 Grisu2 以後，RapidJSON 的 JSON 字符串化（stringify）性能遠超其餘 JSON 庫。

整數除法優化

在原始論文[5] Grisu2 的digit_gen()函數（對應於 double-conversion 實現中的DigitGen()）中，有一段代碼是用於生成足夠的十進位數位：

uint32_t p1 = /*...*/;
int kappa = 10;
uint32_t div = 1000000000;
while (kappa > 0) {
    d = p1 / div; // 第一個除法
    if (d || *len)
        buffer[(*len)++] = static_cast<char>('0' + static_cast<char>(d));
    p1 %= div;
    kappa--;
    div /= 10;    // 第二個除法
    // ...
}

在性能剖析時，我發現這段代碼的32位無號整數除法成爲一個瓶頸。

翻查一下資料[6]，Intel Haswell 架構的 DIV 32位除法指令的延遲（latency）是 28 個週期，吞吐率是 10 個週期。做爲比較，同一架構下 MUL 32位乘法指令的延遲只是 4 個週期，吞吐率只是半個週期。

在許多書籍（如[7]）也會談及，當除數爲常數時，能夠把除法變成乘以除數的倒數。如今的編譯器都會自動作這個優化。事實上，在上面的代碼裏，第二個除法（div /= 10）中的除數（10）就是常數，編譯器會自動把它優化成64位乘法及右移指令，例如 clang 在 x86-64 目標下：

mov    esi, 0xcccccccd  ; 3435973837
imul   rsi, rax
shr    rsi, 0x23        ; 35

注意到 3435973837 \cdot 2^{-35} = 0.10000000000582076609134674072265625，而這個精度足以應付任意32位無號整數除以10。

咱們再分析原來的代碼，會發現，其實除數 div 等於 10^{\kappa - 1}，咱們可使用常數除數令編譯器進行上述的優化：

while (kappa > 0) {
    uint32_t d = 0;
    switch (kappa) {
        case  9: d = p1 /  100000000; p1 %=  100000000; break;
        case  8: d = p1 /   10000000; p1 %=   10000000; break;
        case  7: d = p1 /    1000000; p1 %=    1000000; break;
        case  6: d = p1 /     100000; p1 %=     100000; break;
        case  5: d = p1 /      10000; p1 %=      10000; break;
        case  4: d = p1 /       1000; p1 %=       1000; break;
        case  3: d = p1 /        100; p1 %=        100; break;
        case  2: d = p1 /         10; p1 %=         10; break;
        case  1: d = p1;              p1 =           0; break;
        default:;
    }
    if (d || *len)
        buffer[(*len)++] = static_cast<char>('0' + static_cast<char>(d));
    kappa--;
    // ...
}

這樣的話，編譯器就會把除法都變成乘法及右移。因爲這個switch的 cases 是密集的，編譯器也可以使用 branch table 很好地優化它。不過缺點就是生成的機器碼較原來多。

順便一提，當整數模除運算（modulo operation）與除法成對出現時，而操做數相同，那麼編譯器會把模除運算生成一個乘法及減法：

\begin{align*} c &= \left\lfloor a \mathbin{/} b \right\rfloor \\ r &= a - b \cdot c \end{align*}

這一節沒有使用 intrinsic 或其餘底層優化，只是手動把除法用另外一個方式表達，就能達到有效的性能提高。

底層優化

然而，爲了進一步提高性能，RapidJSON 也會盡可能針對編譯器／目標平臺作一些優化。

例如，在 Grisu 算法中，須要實現一個自定義的浮點數類型（DiyFp），這個類型的乘法須要使用到一個 64 位整數乘法，並得到 128 位結果，而後進行數值修約（rounding）。然而，標準 C++ 中並無得到 128 位乘法結果的方法，所以[5]提供了一個通用實現方法：

struct DiyFp {
    DiyFp operator*(const DiyFp& rhs) const {
        const uint64_t M32 = 0xFFFFFFFF;
        const uint64_t a = f >> 32;
        const uint64_t b = f & M32;
        const uint64_t c = rhs.f >> 32;
        const uint64_t d = rhs.f & M32;
        const uint64_t ac = a * c;
        const uint64_t bc = b * c;
        const uint64_t ad = a * d;
        const uint64_t bd = b * d;
        uint64_t tmp = (bd >> 32) + (ad & M32) + (bc & M32);
        tmp += 1U << 31;  /// mult_round
        return DiyFp(ac + (ad >> 32) + (bc >> 32) + (tmp >> 32), e + rhs.e + 64);
    }
    // ...
    uint64_t f;
    int e;
};

上面的通用實現可用於支持32位乘法並得到64位結果的編譯器。然而，在許多 64 位 CPU 架構下，64位乘數是能夠得到128位結果的。咱們能夠針對各編譯器，使用 intrinsic 或擴展類型來實現這個函數：

DiyFp operator*(const DiyFp& rhs) const {
#if defined(_MSC_VER) && defined(_M_AMD64)
        uint64_t h;
        uint64_t l = _umul128(f, rhs.f, &h);
        if (l & (uint64_t(1) << 63)) // rounding
            h++;
        return DiyFp(h, e + rhs.e + 64);
#elif (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6)) && defined(__x86_64__)
        __extension__ typedef unsigned __int128 uint128;
        uint128 p = static_cast<uint128>(f) * static_cast<uint128>(rhs.f);
        uint64_t h = static_cast<uint64_t>(p >> 64);
        uint64_t l = static_cast<uint64_t>(p);
        if (l & (uint64_t(1) << 63)) // rounding
            h++;
        return DiyFp(h, e + rhs.e + 64);
#else
        // 通用實現
#endif
    }

除此之外，這個 DiyFp 類型還要支持浮點數正規化（normalization）的操做，即把尾數（mantissa）的最高位變成 1（這個浮點數沒有隱藏最高位）：

DiyFp Normalize() const {
        DiyFp res = *this;
        while (!(res.f & (static_cast<uint64_t>(1) << 63))) {
            res.f <<= 1;
            res.e--;
        }
        return res;
    }

許多 CPU 也支持 Find first set 指令，執行一個指令便能掃瞄到最高爲1的位：

DiyFp Normalize() const {
#if defined(_MSC_VER) && defined(_M_AMD64)
        unsigned long index;
        _BitScanReverse64(&index, f);
        return DiyFp(f << (63 - index), e - (63 - index));
#elif defined(__GNUC__) && __GNUC__ >= 4
        int s = __builtin_clzll(f);
        return DiyFp(f << s, e - s);
#else
        // 通用實現
#endif
    }

因爲 gcc/clang 的內置函數能對不一樣目標平臺生成最優的代碼，使用起來更爲方便。

結語

因爲篇幅的關係，本文並無仔細地解釋 Grisu 的算法，而我也不認爲能比原文[5]更淺白地介紹當中的原理。本文只是談到兩個優化方式，一個是利用常數除數令編譯器能進行優化，而另外一種優化則是因爲 C++ 標準沒法使用一些 CPU 提供的功能，而要採用編譯器或平臺相關的優化方法。

參考文獻

[1] Coonen, Jerome T. "an Implementation Guide to a Proposed Standard for Floating-Point Arithmetic." Computer 13.1 (1980): 68-79.

[2] Steele Jr, Guy L., and Jon L. White. "How to print floating-point numbers accurately." ACM SIGPLAN Notices. Vol. 25. No. 6. ACM, 1990. http://kurtstephens.com/files/p372-steele.pdf

[3] Gay, David M. "Correctly rounded binary-decimal and decimal-binary conversions." Numerical Analysis Manuscript 90-10 (1990). http://ampl.com/REFS/rounding.pdf

[4] Burger, Robert G., and R. Kent Dybvig. "Printing floating-point numbers quickly and accurately." ACM SIGPLAN Notices. Vol. 31. No. 5. ACM, 1996. http://www.cs.indiana.edu/~dyb/pubs/FP-Printing-PLDI96.pdf

[5] Loitsch, Florian. "Printing floating-point numbers quickly and accurately with integers." ACM Sigplan Notices 45.6 (2010): 233-243. http://www.cs.tufts.edu/~nr/cs257/archive/florian-loitsch/printf.pdf

[6] Granlund, Torbjörn. "Instruction latencies and throughput for AMD and Intel x86 processors, 2014." https://gmplib.org/~tege/x86-timing.pdf

[7] Warren, Henry S. Hacker's delight. Pearson Education, 2012.