65.9K
CodeProject 正在变化。 阅读更多。
Home

“双精度”类型

Daniel Pfeffer
starIconstarIconstarIconstarIcon
emptyStarIcon
starIcon

4.96/5 (24投票s)

2015年3月14日

CPOL

11分钟阅读

viewsIcon

67053

downloadIcon

1379

一种廉价实现高浮点精度的方法。

引言

许多科学应用需要比“双精度”更高的精度。然而，目前还没有硬件实现“四精度”，并且软件实现成本通常是 prohibitive 的。本文介绍了一种使用“双精度”类型实现的有效折衷方案。

免责声明

当前的实现尚未经过彻底测试。尽管已尽一切努力确保准确性，但仍可能存在代码行为不符合预期的情况。在未进行彻底测试之前，请勿将您的职业生涯或声誉托付给此代码！

增强功能

版本 1.1

  1. 大多数函数（acos()、atan()、asinh() 函数仍需改进）的精度有所提高。
  2. 已向解决方案添加了一个名为 TestQD 的测试程序，以提供函数的“基本健全性检查”。
  3. 错误修复 - fmod()、round()、流输出。
  4. fma() 函数已重命名为 Fma()，以避免与 VS2013 运行时库不兼容。

版本 1.1.1

  1. 错误修复 - pow(负数, 整数) 现在返回正确的值
  2. 错误修复 - 0.0 / 0.0 == NaN
  3. 提供的 Visual Studio 解决方案使用 Visual Studio 2013。对于仍在使用 VS 2012/2010 的用户，版本 1.1 的解决方案和项目文件仍应有效（未更改/添加任何文件名）。

版本 1.1.3

  1. 错误修复 - (unsigned) int、(unsigned) long 和 (unsigned) long long 构造函数
  2. 错误修复 - toInt()、toLong() 和 toLongLong() 转换器

背景

在模拟布朗运动时，我偶然发现了这个扩展精度的解决方案。该程序要求通过极其微小的量更新模拟中原子的位置，这意味着增量被四舍五入（有时甚至完全移位）。确保计算准确的唯一方法是使用更高的精度。

使用代码

代码以 Visual Studio 2010 项目的形式提供。文件 qd\config.h 和 include\qd\qd_config.h 包含宏，应进行修改以支持不同的编译器。特别是，其中一些 API 在标准的 C++11 库中存在（但在早期版本中不存在）。

除一个例外外，该代码被设计为“双精度”类型的直接替换。例如，如果您的原始代码是

#include <iostream>
#include <math>

int main(void)
{
    double two = 2.0;

    std::cout << std::sqrt( two );
    return 0;
}

那么代码可以重写为使用“双精度”类型，如下所示

#include <iostream>
#include "dd_real.h"

int main(void)
{
    dd_real two = "2.0";

    std::cout << std::sqrt( two );
    return 0;
}

例外情况是初始化为常量。为了将变量初始化为完整的“双精度”，必须将该值作为字符串传递。这是因为未作为字符串传递的值将在编译时四舍五入到“双精度”，这就失去了意义。

IEEE 754 浮点数算术标准 vs. “双精度”

IEEE 754-2008 浮点数算术标准[1] 定义了三种二进制算术类型——binary32、binary64 和 binary128。它们的属性总结在表 1 中。

  Sign 指数 尾数 硬件？
表 1 - IEEE 754 二进制类型
binary32 1 位 8 位 24 位 是 - 'float'
binary64 1 位 11 位 53 位 是 - 'double'
binary128 1 位 15 位 113 位 N
  1. 在二进制表示中，尾数的 MSB 是隐式的。这解释了上述表示中的“额外”一位。
  2. 虽然 binary32 和 binary64 类型通常在硬件中实现，但 binary128 目前没有硬件实现。鉴于 binary64 的软件实现通常比硬件实现慢 10 到 100 倍，这意味着 binary128 的软件实现将 prohibitive 缓慢。

双精度类型是 binary128 的额外精度和 binary64 的速度之间的折衷。它依赖于 IEEE 754 浮点数的属性，通过成对的 binary64 值来表示数字，以便该对的总和代表该数字。为了获得最佳结果，对中的位不应有重叠，即较小数的绝对值不应大于较大数绝对值的 2-53。

通过这种方式实现的精度略低于 binary128（106 位而不是 113 位），但指数范围略小于 binary64。原因是许多操作仅在对中的任一数字都不发生下溢时才有效。这意味着较大的数字必须至少是 253*epsilon（最小的 binary64 数字）。

在本文中，我将不带证明地介绍基本算法。有兴趣了解证明的读者可以参考相关文献。

前置条件

双精度类型的实现依赖于以下条件

  1. 所有舍入都采用“最近或偶数”

  2. 计算中不涉及“隐藏位”——操作始终舍入到 binary64，不经过中间步骤

[2] 暗示使用 80x87 的标准模式（64 位精度）将不起作用；必须将其设置为 53 位精度。另一方面，使用 SSE2 指令的程序将正常工作。请注意，Windows 默认会将 80x87 设置为 53 位模式，以尝试与 SSE/SSE2 指令兼容。

基本运算[2]

“双精度”算术的实现依赖于能够精确执行某些浮点运算的能力。在以下所有内容中，RN( ) 表示应使用“最近或偶数”舍入执行操作。**必须强调的是，只有在操作过程中不发生溢出或下溢的情况下，这些基本运算才有效。** 例如，两个同号无穷大的加法将产生 (inf, NaN) 的结果，而不是预期的 (inf, 0)。

有可能使这些基本运算的行为符合预期，但这需要额外的测试和（可能的）分支。

基本运算如下

加/减法[3, 4, 5]

对两个 binary64 数字进行加/减运算，返回精确结果——幅度较大的数是加/减运算的舍入结果，较小的数是精确结果的余数。

存在两种变体

(s, t) = Fast2Sum(a, b)[3] – 要求 exponent(a) >= exponent(b)

s = RN( a + b )
t1 = RN( s - a )
t = RN( b - t1 )

(s, t) = 2Sum(a, b)[4,5] – 对 a, b 无先决条件

s = RN( a + b )
t1 = RN( s - b )
t2 = RN( s - t1 )
t3 = RN( a - t1 )
t4 = RN( b - t2 )
t = RN( t3 + t4 )

Veltkamp 分裂[3]

取一个精度为 *p* 的浮点数 x，将其分裂成两个浮点数 (xh, xl)，使得

  • xh 具有 *p-s* 位有效数字
  • xl 具有 *s* 位有效数字

这假设计算过程中没有发生溢出。二进制浮点数的一个怪癖是 xl 实际上可以放入 *s-1* 位（xl 的符号位用作附加位）。对于 binary64 (*p*=53)，使用 *s*=27 意味着数字的两个部分都可以放入 26 位。

(xh, xl) = Split(x, s)

Define: C = 2^s + 1
t1 = RN( C * x )
t2 = RN( x – t1 )
xh = RN( t1 + t2 )
xl = RN( x – xh )

乘法[2, 3]

将两个 binary64 数字相乘，返回精确结果——幅度较大的数是乘法的舍入结果，而较小的数是精确结果的余数。

(r1, r2) = 2Mult(x, y) – 也称为 Dekker Product

Define: s = 27
(xh, xl) = Split( x, s )
(yh, yl ) = Split( y, s )
r1 = RN( x * y )
t1 = RN( -r1 + RN( xh * yh ) )
t2 = RN( t1 + RN( xh * yl ) )
t3 = RN( t2 + RN( xl * yh ) )
r2 = RN( t3 + RN( xl * yl ) )

如果存在硬件融合乘加指令，则可以实现以下速度更快的方法

(r1, r2) = 2MultFMA(x, y)

r1 = RN( x * y )
r2 = fma( x, y, -r1 )

“双精度”类型

在以下内容中，我将基于 Bailey 等人在其 QD 库[6,7] 中发布的代码。虽然基本实现是他们的，但我对代码进行了以下增强

函数式

  • 三角函数和双曲函数的实现已重写，以提高精度和速度。特别是，现在使用 fdlibm[8] 中的优秀代码来执行三角函数参数的归约
  • 添加了原始 QD 库中未定义的其他函数（例如 expm1、logp1、cbrt 等）

代码组织

  • 原始代码可以配置为使用更快但不那么精确的基本函数（加法、乘法、除法）。已删除这些函数的不精确版本
  • 除类 *dd_real* 外，所有“双精度”类型的代码都已移入 *qd* 命名空间
  • 所有辅助函数都已移入 *dd_real* 类，以免污染全局命名空间
  • 如果 C++ 标准库为“双精度”类型定义了数学函数，那么它也会在 *std* 命名空间中为 *dd_real* 类型定义。

左键单击 gadget 并拖动以移动它。左键单击 gadget 的右下角并拖动以调整其大小。右键单击 gadget 以访问其属性。

必须强调的是，“双精度”类型是对 106 位浮点类型的近似。它在许多方面与 IEEE 754 binary128（113 位）类型不同

  • 不保证正确舍入

尽管已尽一切努力提供“忠实”舍入（即返回最接近无限精度值的两个值之一），但不能保证返回的值确实是两个值中最接近的。

  • “抖动”精度

像 1 + 2-200 这样的数字很容易在“双精度”格式中表示，但不是有效的 IEEE 754 binary128 数字。这种“抖动”精度会使错误分析和其他问题复杂化。

  • 无穷大和 NaN

无穷大和 NaN 算术的完整实现需要对所有运算阶段的正常操作数和正常结果进行测试。可以编写一个廉价的测试，用于处理非正常（无穷大/NaN）和大多数正常操作数，但幅度最大的 (1 - 2-53) 内的结果可能会触发错误的溢出条件。

  • 与浮点数相关的定理和引理

IEEE 754 浮点类型遵循大量的定理，这些定理使得错误分析等更容易处理。大多数这些定理不适用于“双精度”

基本算术

基本思想是使用上面定义的精确运算来计算精确运算的近似值。请注意，只有结果的两个分量之和才有意义。特别是，结果的高位不一定等于对高位进行运算的结果，例如 zh <> xh + yh！

请注意，存在运算符的优化版本，例如将“双精度”添加到“双精度”等等。

加/减法

(zh, zl) = ddAdd( (xh, xl), (yh, yl) )

(sh, sl) = 2Sum(xh, yh)
(th, tl) = 2Sum(xl, yl)
c = RN(sl + th)
(vh, vl) = Fast2Sum(sh, c)
w = RN( tl + vl )
(zh, zl) = Fast2Sum(vh, w)

乘法

(ph, pl) = ddMul( (xh, xl), (yh, yl) )

请注意，xl 和 yl 的乘积未计算。这是因为它永远不会对结果的低位部分做出贡献。

(ph, pl) = 2Prod(xh, yh)
pl = RN(pl + RN(xh * yl))
pl = RN(pl + RN(xl * yh))


除法

(ph, pl) = ddDiv( (xh, xl), (yh, yl) )

不存在执行除法的精确方法。取而代之的是，通过 xh / yh 计算初步近似值，然后通过牛顿-拉夫逊迭代进行细化。

多项式函数

sqrt() 和 cbrt() 使用牛顿-拉夫逊迭代的变体来计算。

hypot() 使用一种确保所有值都精确的算法来计算。

pow() 计算如下

  1. 如果指数是整数且小于 LLONG_MAX，则通过连续平方和乘法计算该值。
  2. 如果指数是非整数但小于 LLONG_MAX，则将指数分解为整数和非整数部分。整数指数按上述方法计算，分数指数计算为 exp( x* log( frac(y) )。
  3. 如果指数大于 LLONG_MAX，则计算为 exp( x * log( y ) )

三角函数

对于三角函数，首先使用 fdlibm[8] 中提供的 pi/2 余数算法将参数归约到 [-pi/4, pi/4] 范围。然后将其除以 8，得到 [-pi/32, pi/32] 范围内的值。正弦和余弦计算为多项式，正切计算为正弦与余弦之比。

对于反函数，执行牛顿-拉夫逊迭代以计算例如 tan( y ) – x = 0 的 y 值。

对数

对数函数（log、log10、log2）通过将 x 乘以适当的常数，然后计算 log2(C * x) 来计算。我们使用 log2() 作为基本函数，因为它允许我们轻松地提取结果的整数部分（使用 frexp() API），只剩下分数对数的计算。

指数

函数 exp、exp2、expm1 都通过将参数归约到小于 1 并计算泰勒级数来计算。对于 expm1，不添加泰勒级数的零项。

请注意，如果 expm1 的参数足够大，结果将与 exp(x) - 1.0 相同。

双曲函数

除 sinh(x) 在 x 接近零时外，这些函数都通过使用适当的 exp() 函数来计算。对于 sinh(x) 在 x 接近零时，使用函数的泰勒级数来尝试减少灾难性抵消。

实用函数

ldexp、frexp、fmod、modf、fma、copysign 等都存在，并且其定义与 C++ 标准中的类似函数相同。

std::numeric_limits 已完全定义。

常量

函数 dd_pi()、dd_e()、dd_ln10() 等返回预期的值，四舍五入到最近或偶数。

已知bug

  1. 需要提高某些反三角函数（acos()、atan()）和反双曲函数（asin()）的精度。

参考文献

  1. IEEE 浮点数算术标准。IEEE 标准 754-2008，2008 年 8 月

  2. 浮点数算术手册。J.M. Muller 等人，Birkhauser 201x

  3. 扩展可用精度的浮点技术。T.J. Dekker，《Numerische Mathematik》18(3):224-242, 1971

  4. 《计算机程序设计艺术》第 2 卷。D. Knuth，Addison-Wesley Reading MA，第 3 版，1998 年

  5. 准双精度浮点加法。O. Moller. BIT 5:37-50, 1965

  6. 四精度浮点算术算法。Y. Hida、X.S. Li、D.H. Bailey，《第 15 届 IEEE 计算机算术会议论文集》第 155-162 页，2001 年

  7. 原始 QD 库可在 http://crd.lbl.gov/~dhbailey/mpdist 获取

  8. fdlibm 库可在 http://www.netlib.org/fdlibm 获取

  9. 基本函数软件手册。W. J. Cody Jr. & W. Waite，Prentice-Hall，1980 年

历史

2015-03-14 Original version
2015-03-31 Bug fixes, basic test program
© . All rights reserved.