使用四舍五入到最近模式的区间算术（n 的第 1 部分）

摘要

传统上，区间算术一直受到利用其优势编写的程序性能低下的困扰。我提出了一种基于 IEEE Std 754-2008 的方法，该方法具有显著提高速度的潜力。

引言

区间算术是试图回答计算精度问题的一种尝试——给定一组浮点运算，结果的精度如何？

区间算术用一组数字的范围替换单个浮点数，我们知道这个范围包含了该值。例如，使用公制卷尺测量一个 3m x 4m 的房间，结果可能精确到 0.5cm，这将被表示为 [2.995, 3.005] x [3.995, 4.005]。区间上的任何运算都以保持此属性的方式执行。例如，

add(a, b) == add([a_inf, a_sup], [b_inf, b_sup]) == [a_inf + b_inf, a_sup + b_sup].

基本实现

区间算术的基本思想是对每次操作执行两次——一次在“向下舍入”（RD）模式下，一次在“向上舍入”（RU）模式下。虽然这保证能正常工作，但存在以下缺点：

• 每次操作需要三次舍入模式的更改——切换到 RU、切换到 RD，然后再切换回默认模式（通常是“舍入到最近”或 RN）。

• 在大多数现代处理器上，这些切换中的每一次都非常耗时。例如，在较新的 Intel 处理器上，读取 SSE 舍入模式很慢，而设置它则是一个串行化操作！

RU 方法

对于某些操作来说，一个很好的方法是利用 RD( a ) == -RU( -a ) 的事实。为了做到这一点，我们存储区间为 [-a, b] 而不是 [a, b]。显然，加法/减法可以直接用 RU 执行，而乘法和除法可以通过一些额外的符号位翻转来执行。

这具有以下优点：

• 每次操作只需两次舍入模式切换。

• 能够为加法/减法使用 SIMD 指令。

然而，它不能普遍使用。例如，sqrt(-a) == NaN，而不是 -sqrt(a)。

“舍入到最近”方法

区间算术的问题归结为选择问题。给定一个区间边界的无限精度结果 A，当计算区间下界时，我们需要向下舍入；当计算区间上界时，我们需要向上舍入。

我们知道 RN( A ) 等于 RD( A ) 或 RU( A )，但哪一个呢？

事实证明，这些信息是可以获得的。此外，使用“舍入到最近”模式是使此方法奏效的必需条件。

精确变换

给定两个浮点数 a, b，可以如下计算精确变换：

• TwoSum() - (a, b) 转换为 (s, e)，其中 s = RN(a + b) 且 e = a + b - s

s = RN(a + b)
a' = RN(s - b)
b' = RN(s - a')
da = RN(a - a')
db = RN(b - b')
e = RN(da + db)

TwoSum() 算法需要 6 次浮点运算。

• TwoProd - (a, b) 转换为 (p, e)，其中 p = RN(a * b) 且 e = a * b - p

Split(x, n)

Require: C = 2^n + 1, where n = (bits in mantissa + 1)/2

t1 = RN(C * x)
t2 =  RN(x - t1)
xh = RN(t1 + t2)
xl = RN(x - xh)

TwoProd(x, y)

(xh, xl) = Split(x, s)
(yh, yl) = Split(y, s)
p = RN(x * y)
t1 = RN(-p + RN(xh * yh))
t2  = RN(t1 + RN(xh * yl))
t3 = RN(t2 + RN(xl * yh))
e = RN(t3 + RN(xl * yl))

TwoProd() 算法需要 17 次浮点运算。

如果 e 为负，则 RN( a + b ) == RU( a + b )，并且 RD( a + b ) == prev( s )。如果 e 为正，则 RN( a + b ) == RD( a + b )，并且 RU( a + b ) == succ( s )。类似的论证也适用于 TwoProd()。

基本函数

TwoSum() 和 TwoProd() 过程足以执行区间加法/减法和乘法。如果存在 fma() 操作（在 IEEE 754-2008 中是强制性的），则可以计算以下函数：

• TwoProdFMA(a, b)

p = RN(a * b)
e = fma(a, b, -p)

TwoProdFMA() 算法执行 TwoProd()，只需要 2 次浮点运算，而不是 17 次。

• Division() - q = RN(a / b); r = -fma(q, b, -a)
• Reciprocal() - q = RN(1.0 / b); r = -fma(q, b, -1.0)
• Sqrt() - q = RN( sqrt(a) ); r = -fma(q, q, -a)

如果 r 为负，则结果 = RN( A ) == RU( A ) 且 RD( A ) == prev( 结果 )。如果 r 为正，则结果 = RN( A ) == RD( A ) 且 RU( A ) == succ( 结果 )。

请注意

• prev(x) 和 succ(x) 是 IEEE Std 754 所必需的，在 C++ 2011 中实现为 std::nextafter()。
• 如果硬件中没有 fma() 操作，则可以使用 TwoProd() 和 TwoSum() 操作在软件中模拟它。这总共需要 35 次浮点运算——17 次用于 TwoProd()，18 次用于所需的 3 次 TwoSum() 操作。

测试

测试了以下区间算术的 C++ 实现：

• “基本”区间算术——根据需要将舍入模式设置为 RU、RD 和默认模式。

• “RU”区间算术——将 [a, b] 存储为 [-a, b]，在许多情况下消除了将舍入模式设置为 RD 的要求。

• “舍入到最近”区间算术——实现了本文定义的算法。

在所有情况下，如果 SIMD 指令可用且有用，则会使用它们。

测试涉及 100,000,000 次具有随机操作数的运算，如下所示：

1. 加法
2. 乘法
3. 除法（分母区间不包含 0.0）
4. 平方 sqr(x) = (x*x)
5. 平方根
6. hypot(x, y) = sqrt(sqr(x), sqr(y))。

选择此 hypot() 函数是因为它包含了加法（可能通过 RU 方法优化）、平方（也可能通过 RU 方法优化）以及平方根（只能通过“舍入到最近”方法优化）。它也可以用简单的 C++ 代码轻松计算。

使用 Visual Studio 2013 以 x64 Release 模式和 /fp:strict 编译器选项编译了代码。请注意，如果需要获得精确的结果，使用此选项至关重要。

在 Intel Core i7 2670QM 处理器上的结果（每操作纳秒）如下：

 

速度测试结果

	Basic	RU	舍入到最近
加法	104	79	25
乘法	230	194	61
除法	280	217	94
平方	238	185	42
平方根	139	83	44
hypot()	804	702	173

 

 

 

 

 

 

 

 

 

 

 

 

注释

• 这些时间是在减去循环开销和生成随机输入所需的时间后得出的。

• 与最新的 Intel 处理器不同，这款 2670QM 处理器没有硬件 fma()，必须在软件中模拟它。这显著降低了“舍入到最近”方法在乘法、除法、平方和平方根运算方面的优势。

摘要

使用舍入到最近模式的区间算术库是实用的，并且在当前流行的硬件上——比定向舍入方法快得多（3 倍以上）。

未来方向

一个包含 IEEE-1788 标准所需的大部分前向函数的区间数学包目前正在开发中。

参考文献

1. IEEE Std 754-2008

2. IEEE Std 1788-2015

3. Handbook of Floating-Point Arithmetic, Jean-Michel Muller et al., Birkhäuser

4. Interval Arithmetic: from Principles to Implementation, T. Hickey, Q. Ju, M.H. van Emden

5. Interval Operations in Rounding to Nearest, S.M. Rump, P. Zimmerman, S. Boldo, G. Melquiond