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

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

Grisu是什麼

Grisu 是把浮點數轉換為字串的演算法。在 Chrome 裡執行這段 JavaScript 實際上就呼叫了 Grisu：

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

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

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

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

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

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

// 除 +/-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" 作結尾。

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）與除法成對出現時，而運算元相同，那麼編譯器會把模除運算生成一個乘法及減法：

這一節沒有使用 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.