精确复制浮点值的重复加法

Question

是否有一种亚 O(n) 的方法¹ 可以复制 IEEE 754 算术中浮点值的重复求和？

也就是说，是否有一种更快的方法来完全²复制以下伪代码的结果：

f64 notQuiteFMA(f64 acc, f64 f, bigint n) {
    for (bigint i = 0; i < n; i++) {
        acc = acc + f;
    }
    return acc;
}

？³

在实际算术中，答案很简单 -

return acc + f*n

就足够了。然而，在有限精度下，答案要复杂得多。作为一个例子，请注意

notQuiteFMA(0, 1, 1 << 256)

约为 ⁴ 2^53，而不是人们期望的无限精度的 2^256。 （这个例子还说明了为什么加快计算速度可能是可取的 - 计数到 2^256 是相当不切实际的。）

1. 我知道这整个事情在技术上可以在 O(1) 时间内完成⁵，使用 2^128 条目查找表，每个条目都有一个 2^64 元素的数组；请忽略此类银河算法。 （鸽巢原理 - 在进入循环之前，重复加法不能超过 2^64 步 - 因此“只需”为

   acc

   和

   f

   的每个初始值存储整个前缀和最终循环中的步数。 )
2. 给定当前舍入模式的位精确（至少忽略非信号 NaN 的蠕虫）。
3. 将此视为伪代码 - 即所描述的双精度浮点加法序列。我知道在某些编程语言中，浮点计算可以以更高的精度完成/重新排序以使事情更加SIMD友好/等等；我对这个问题的这些案例不感兴趣。
4. 64 位浮点数有 52 位尾数；当加法不足以增加尾数且结果向下舍入时，

   x + 1 == x

   首先出现在 2^53 左右。我相信确切的值取决于舍入模式。
5. 假设一个计算模型，其中 bigint 与常量有界值的模至少为 O(1)。

Answer 1

我可以回答这个问题，尽管我认为你不会喜欢它。

一旦数字达到 53 个有效位，该值将停止递增。我们可以计算一下。下面是一个 Python 脚本，显示了向浮点数加 1 不再改变浮点数的点：

import struct

j = 1
for i in range(1,58):
    j = 2*j
    r = float(j)
    s = r+1
    if r == s:
        pat = struct.pack('>d',r)
        print(pat.hex())
        print(i,j,r,s)
        break

输出：

4340000000000000
53 9007199254740992 9007199254740992.0 9007199254740992.0

第一行是该值的十六进制表示。