长浮点数(long double)的陷阱

前言

此文转载自 Prelert 的博客,作者为 David。这篇文章目前只能在 WebArchive 找到了,因此我把此文转载并翻译到我的博客上。

我是在搜索为什么 Rust 没有对应 C/C++ 中 long double 的数据类型的时候看到了这篇博文,Rust 不提供对应的数据类型造成了一些互操作性的问题(参见这里这里)。与此相对的是,Zig 和新发布的 Carbon 语言都支持 f16f128 数据类型(其中 Zig 还支持 f80,Carbon 还支持 bfloat16)。不过这倒是不令人意外,因为 Zig 和 Carbon 都以与 C/C++ 的极致兼容性为卖点。这篇博客也许能解释一部分 Rust 不支持更高精度浮点数的原因。

正文

C++ 提供三种浮点数据类型:float, doublelong double。关于这些类型,C++11 标准只提到了:

double 类型需要提供至少与 float 同等的精度,而 long double 需要提供至少与 double 同等的精度。

float 类型所有支持的数值是 duoble 的子集;double 类型所有支持的数值是 long duoble 的子集。各浮点类型的数据表示方法取决于具体的实现。

但是,几乎所有 C++ 编译器都带一个 C 编译器,而 C99 标准的 F 附件则有更详细的规定:

  • float 类型对应 IEC 60559 标准的单精度浮点数
  • double 类型对应 IEC 60559 标准的双精度浮点数
  • long double 类型对应 IEC 60559 标准的扩展精度浮点数,或者非 IEC 60559 的扩展精度浮点数,或者 IEC 60559 的双精度浮点数。

只有一个彻头彻尾的 M 才会编写一个在 C++ 和 C 部分使用不同浮点类型的 C++ 编译器,因此实际上 C++ 也要服从同样的规定。我所使用果的所有近 20 年内的 C++ 编译器都使用 IEC 60559(也就是 IEEE 754)所规定的单精度和双精度浮点数,但是在实现最后一个类型 ——long double—— 上这些编译器会有不一致,这也导致了一些问题。

在我的整个软件开发生涯中,我遇到过的几次与 long double 类型相关的问题,可以被归为两类:

  1. 缺乏测试
  2. 可移植性

缺乏测试

在去年底我记录了一个可以归为第一类的问题。一个 glibc 在 x86_64 平台上的 powl() 函数实现中的 bug 有五年多未得到修复。我感觉如果这个 bug 是在用的更广泛的 pow()函数中,那么会有更多人感到惊奇然后有人会更快修复它。因为这个函数的 long double 版本用的人较少,因此这个 bug 就烂在那了。

另一个 long double 缺乏测试的例子是我在加入 Prelert 之前碰到的,与 IBM 的 xlC/C++ 编译器相关在 AIX 系统上的问题。调用这个编译器的时候使用的名字(硬链接)决定了它的行为方式,当使用名称 xlC128_r 调用编译器时,它使用 128 位的 long double 表示。然后有一段时间,即使是最简单的程序编译都会崩(core dump)。尽管 bug 报告里面提到了一个调用 fork() 的例子,但其实如果打开了 -brtl 开关,最简单的 "hello world" 程序也能崩!显然所有的测试都是在使用其他更常用的名字调用编译器时完成的(这些情况下 long double 不是 128 位)。

可移植性

而关于可移植性,一些需要注意的坑有:

  1. 微软的 Visual C++ 使用 IEEE 754 双精度浮点数表示 long double—— 跟 double 一样(C99 标准所允许的第三种情况)。因此在你的代码中区分 long doubledouble 毫无意义如果你只用微软的 VC++ 进行编译。但如果你的代码需要支持别的平台并且你仍然使用 long double,那么你就给你的代码里面掺入了一个很关键的,与平台相关的行为区别,它会让你吃亏的。大部分其他 x86 编译器都把 long double 看作 x87 所使用的 80 位扩展精度浮点数。
  2. SPARC 芯片上(我知道它已经快挂了),long double 类型使用 128 位的表示,但是默认情况下,编译器会生成软件实现,而非硬件实现的(浮点数)操作。这个情况可以回溯到大部分 SPARC 芯片都不支持这样的操作,然后使用中断来实现的时候。在软件层面实现浮点操作比去响应这些中断要快。但是软实现的浮点操作比硬件加速的浮点操作慢好几个数量级 —— 我们发现一些单元测试会慢 20 倍,而且这些并不是单纯的做 long double 运算的测试。这是一个牺牲性能来换取代码(在编译和正确性方面)可移植的例子。

参考其他可移植的语言很有指导意义。Java 有与 IEEE754 的单 / 双精度浮点数对应的 floatdouble 类型(并且不像 C++,Java 标准对如何实现浮点数运算非常明确)。Java 没有给程序员提供 long double 类型,大概是因为我上面提到的可移植性问题(尽管这个标准允许在中间计算过程中使用 x87 扩展精度的格式)。Python 只有一个 float 类型,并且 “通常用 C 语言的 double 实现”。因此,如果你的整体系统包含使用其他语言编写的组件,那你弃用 long double 可以避免数据交换中的问题。同样的道理也适用于在数据库中存储浮点数 —— 例如 PostgreSQL 提供 realdouble,对应 IEEE 754 的单 / 双精度浮点数。

最后,一个在 x86 CPU 上只用 floatdouble 的好处是,编译器可以利用上 CPU 的 SSE 单元,然后两次或四次(浮点)运算有机会并行完成。这时使用 64 位调用约定的把函数参数传到寄存器中,那之后就可以直接使用 SSE 寄存器了。反之,long double 变量只能在 x87 浮点运算单元中使用,并且不会使用寄存器传参,从而让程序变慢了。

结论

有些人可能会说,使用 long double 可以提高结果的精度。这可能是对的,但是无论一个固定精度的浮点数有多少有效位数,它都会有有效位数损失的情况(如果使用的算法不好的话)。使用扩展精度而不是双精度可能会在某些情况下避免这些问题,但是长期来看唯一的解决方法是使用更适合计算机的算法,或者设法检测有效位数损失,并且把结果替换成合适的值。

在我看来,如果你想编写可移植的 C++ 代码,使得它不仅可以在多个平台运行,并且在某些平台上没有很夸张的性能问题,你最好避开 long duoble。这也是我们在 Prelert 做法 —— 我们的 C++ 代码不使用 long duoble,并且在使用 Boost 的时候我们定义 BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS 这个宏,让 Boost.Math 也不使用。