代码改变世界

RapidJSON 代码剖析(四):优化 Grisu

2015-06-30 13:03 by Milo Yip, ... 阅读, ... 评论, 收藏, 编辑

我曾经在知乎的一个答案里谈及到 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),fpconvgrisu2milo 都是 Grisu2 的实现,doubleconv 是 V8 的 Grisu3 实现。milosprintf 的加速比约是 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.