浮点运算符正确舍入实现

Floating Point operators correctly rounded implementation

我正在用 C 实现一个浮点库,基本上它是基于 IEEE 标准的。我已经开始添加了。

我在理解四舍五入方面遇到了一些问题。

也许一个例子可以帮助

x = -3.652248e-11 
y = 1.263346e-10 
Cz = 8.981213e-11
Mz = 8.981214e-11
Cz = 8.98121e-11
Mz = 8.98121e-11
x = ae20a0a7
y = 2f0ae808
Cf = 2ec57fbc
Mf = 2ec57fbd
Cz = 457fbc
Mz = 457fbd
Ce = 5d
Me = 5d

上面的xy是随机输入值,Cz是C运算的结果+最后是Mz我实施的结果。

CfMf是最终结果的位表示。如您所见,最后一点不同,我不明白为什么。我从一本浮点运算手册中获得了实现灵感。

我想我不明白的是舍入的实际执行方式。我的加法算法是基于恒等式

x + y = (-1)^{sx}2^Ex(|x| + (-1)^(sx xor sy) |y| 2^{Ey-Ex})

如果我命名数量

|z| = (|x| + (-1)^(sx xor sy) |y| 2^{Ey-Ex})

基本上,当我需要 post 使用左移对结果进行归一化时,问题就出现了,在我的情况下要小心 |z| 总是正数。在这种情况下应该应用哪种舍入技术?

我的穆勒等人的副本。是借给朋友的,所以我无法仔细检查您具体使用的算法,但会逐步添加您列出的值:

x = 0xae20a0a7 = -b1.01000001010000010100111 * 2^-35
y = 0x2f0ae808 = +b1.00010101110100000001000 * 2^-33

如果我们将 xy 归一化为一个公共指数并相加,我们将得到非归一化的无限精确结果:

  b100.01010111010000000100000 * 2^-35
-   b1.01000001010000010100111 * 2^-35
-------------------------------------
   b11.00010101111111101111001 * 2^-35

现在标准化而不四舍五入:

b1.10001010111111110111100 1 * 2^-34
                          ^
                          rounding point

无限精确的结果恰好是最接近的两个浮点数的一半,所以我们选择偶数一个,向下舍入到

b1.10001010111111110111100 * 2^-34 = 0x2ec57fbc

鉴于这是一个恰好一半的情况,最可能解释为什么你没有得到正确的答案是你没有处理这个 与 even 部分的关系正确的舍入规则。如果您尝试通过简单地添加半个 ulp 并截断来四舍五入,您将得到您正在观察的结果。