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

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

一旦数字滚动到 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

第一行是该值的十六进制表示形式。因此，它是指数为 52 的隐藏 1。

这是一个棘手的问题......

至少我们知道，直到 ，结果总是处于四舍五入到最接近的偶数模式。
例如，参见 3*x+x 总是准确的吗？n=5n*f

但后来我们逐渐失去了精确度......以什么速度？

我认为我们可以预测溢出的情况，所以我不会为此烦恼。因此，我们可以忽略指数部分，而专注于有效部分。

渐进下溢的情况也不是很有趣，因为在偏置指数达到 1 之前，加法将是精确的。

因此，我们可以关注 （1,2（ 区间中的有效数，1 表示无趣。

我尝试遵循 Smalltalk 片段以了解错误边界：

sample := Array new: 60000.
x := 1.0.
00001 to: 10000 do: [:i | sample at: i put: (x := x successor)].
x := 2.0.
10001 to: 20000 do: [:i | sample at: i put: (x := x predecessor)].
x := 1.5.
20001 to: 30000 do: [:i | sample at: i put: (x := x predecessor)].
x := 1.5.
30001 to: 40000 do: [:i | sample at: i put: (x := x successor)].
x := 1.125.
40001 to: 50000 do: [:i | sample at: i put: (x := x predecessor)].
x := 1.125.
50001 to: 60000 do: [:i | sample at: i put: (x := x successor)].
(3 to: 10) collect: [:i |
    | setOfErrors |
    setOfErrors := sample collect: [:f | ((1 to: 1<<i) inject: 0 into: [:acc :ignored | acc+f]) - (1<<i*f) / f ulp] as: Set.
    {i. setOfErrors size. setOfErrors max log2}].

1<<i是左移，所以等价于 .
因此，当 n 是 2 的幂时，对于上述有效样本，误差范围为：2^i

 #(#(3 5 4.0) #(4 9 6.0) #(5 17 8.0) #(6 33 10.0) #(7 65 12.0) #(8 129 14.0) #(9 257 16.0) #(10 513 18.0)) 

这是 时的最大误差，以及所有可能的误差组合。2^(2*(i-1))*ulp(f)n=2^i2^(i-1)+1k*2^i*ulp(f)k=-2^(i-2)...2^(i-2)

在我们达到 2 提高到有效位数之前，有许多不同的可能性......因此，我怀疑您是否可以找到位模式的任何简单功能，仅适用于默认舍入模式。

编辑

使用 Squeak Smalltalk 的 ArbitraryPrecisionFloat 包，我尝试用 p 有效位推广 Float 的公式。

例如，使用 11 位有效：

p := 11.
(3 to: p + 1) collect: [:power |
    | n setOfErrors |
    n := 1 << power.
    setOfErrors := (1 << p to: 2<<p - 1) collect: [:significand |
        | f acc |
        f := (significand asArbitraryPrecisionFloatNumBits: p) timesTwoPower: 0 - p.
        acc := 0. n timesRepeat: [acc := acc + f].
        (n*f-acc/(acc ulp)) asFraction] as: Set.
    {power. setOfErrors max log2 printShowingDecimalPlaces: 2. setOfErrors size}].

我们得到这个结果：

#(#(3 '1.00' 5) #(4 '2.00' 9) #(5 '3.00' 17) #(6 '3.91' 32) #(7 '4.91' 61) #(8 '5.88' 120) #(9 '6.77' 220) #(10 '7.41' 354) #(11 '7.41' 512) #(12 '10.00' 1024))  

请注意，过去的操作数，总和不会再改变。2^(p+1)

到目前为止，存在与上一个样本所示的不同结果。power=floor(p/2)2^(power-1)+1

然后，不同结果的数量增加得更慢，但最后，它达到 ，即每 1 个不同的有效数有 2 个不同的结果。2^(p-1)

这表明这将是相当复杂的。notQuiteFMA

询问 的模式，我们看到它似乎取决于 f 有效性的最后一位，直到总和达到一个新的 binade......(n*f-acc/(acc ulp))ceil(log2(n))

所以，也许你可以问这样的话：

• 的 Binade（指数）与 的 Binade（指数）是多少？n*ff
• 对于 F 的后继者，误差序列是否遵循可预测的模式？n*(f+k*ulp(f))

如果你能肯定地回答，那么你就有机会设计出一个加速功能......notQuiteFMA

结论

我们可以通过单步执行 binades 来实现 O（log n） 解决方案。这是因为一个 binade 中的所有加法，可能除了第一个加法之外，都会以相同的量改变累积总和，因此我们可以很容易地准确地计算出重复加法会产生什么总和，直到 binade 结束。

预赛

binade 是某个整数 e 的区间 [2^^e， 2^e+1） 或 （−2^^e+1， −2^e]。就此答案而言，整个最小指数范围（包括最小指数法线和所有次正态值）将被视为单个二元组。

对于浮点数 x（可以浮点格式表示的数字），将 E（x） 定义为表示 x 时使用的指数。鉴于使用此答案的 binade 的定义，E（x） 表示 binade x 在，除了它的符号。

代码格式的表达式表示浮点运算。x+y 是将 x 和 y 相加的实数算术结果。 是将 和 相加的浮点结果。x 并用于适当的相同值。x+yxyx

浮点加法对于有限值是弱单调的：< 表示 <= ，在任何舍入模式下，除非是无穷大和/或为相反符号的无穷大，从而产生 NaN。yzx+yx+zxyz

讨论

设 A ₀、A ₁ 和 A₂ 是 执行 之前、之后的值。 考虑 E（A ₀） ≥ E（A ₁） = E（A₂） 的情况。设 g 为 A₂−A₁。g 是可表示的，如果 A 1 是正态的，或者如果 A 1 是次正态的，则由 Sterbenz 引理表示，因为 A ₁、A 2 和都是次正态（或零），并且没有低于 A₂ 的 ULP 的位在加法中丢失。（由于 g 是可表示的，因此在计算时不会出现四舍五入误差。AA += ffA2-A1

我们将证明，所有后续的相加都会产生 E（） = E（A₁） 的值增加相同的量 g。fAAA

设 u 为 A₁ 的 ULP。设 r 是 ||模 U.如果 f 为正数，则 s 为 +1，如果 f 为负数，则为 −1。（f 为零的情况是微不足道的。使用有向舍入模式（朝向 −∞、朝向零或朝向 +∞），根据舍入方法和 的符号，任何不改变 E（） 的后续加法中的舍入误差始终为 −r 或始终为 u−r。使用四舍五入到最接近的领带 如果 r < 1/2u，则舍入误差始终为 −sr，否则始终为 s（u−r）。当四舍五入到最接近时，平局为偶数，如果 r ≠ 1/2u，则舍入误差始终为 −sr 或始终为 s（u−r），与平局为 away 一样。fAf

这样一来，四舍五入到最接近的，与 r = 1/2u 有偶数的联系。设 b 是 u 位置的 f 位。在这种情况下，A 1 的低位必须是偶数，因为它是将剩余模 u 为 1/2 u 的 f 与 A 0（必然是 0 模 u，因为 E（A ₀） ≥ E（A 1））相加以产生 A ₁ 的结果.此外，实数结果正好介于两个可表示值之间，因此它被四舍五入到一个偶数（2u 的倍数）。

假设 A₁ 的低位是偶数，我们有两种情况：

b 为 0，因此加法的实数算术结果是 0 u + 1/2 u 模 2 u ≡ 1/2 u，并且四舍五入到偶数，因此如果结果为正，则舍入误差为 −1/2 u，如果结果为负，则舍入误差为 +1/2u。然后对后续添加重复此操作，生成 e （） = E（A₁）。fAA

b 为 1，因此加法的实数算术结果是 1 u + 1/2 u 模 2 u ≡ 11/2 u，并且四舍五入到偶数，因此如果结果为正，则舍入误差为 +1/2 u，如果结果为负，则舍入误差为 −1/2 u。然后对后续添加重复此操作，生成 e （） = E（A₁）。fAA

（要了解为什么我们在断言所有后续加法都加相同量之前先进行二次加法，请观察产生 A₀ 的加法可能涉及具有较低 ULP 的先验位，在这种情况下，实数算术结果可能并不完全介于两者之间 可表示的值，因此 A 0 可以有一个奇数的低位数偶数，使用四舍五入到最接近，并列到偶数，并且 A _2-A 1 不等于 A _{1-A 0}。 A

实现

鉴于上述情况，我们可以计算出在达到下一个 binade 之前可以执行多少次 g，比如 t。然后我们可以计算出 A 在退出当前 binade 之前的确切值 A+tg。然后，一个加法为我们提供了该比合体中的第一个 A，另一个加法要么将我们带入一个新的加法，要么为我们提供执行当前比合计算所需的信息。然后，我们可以对每个 binade 重复此操作，直到达到所需的添加次数。

下面是演示该概念的代码。根据评论，此代码需要一些改进。我希望在时间允许的情况下进行研究，这可能比现在还要几周以上。代码中的函数从 中获取指数，而不是如上面讨论中所述将其限制为最小指数。这可能会减慢代码速度，但不应使其不正确。我希望代码处理次正常现象。无限和 NaN 可以通过普通测试添加。Efrexp

t 的计算，即在退出当前 binade 之前可以执行的加法次数，可以使用更多关于所需精度以及舍入如何影响它的讨论。

#include <math.h>
#include <stdio.h>
#include <stdlib.h>


static float Sum0(float A, float f, int n)
{
    if (n < 0) { f = -f, n = -n; }
    while (n--) A += f;
    return A;
}


//  Return the exponent of x, designating which binade it is in.
static int E(float x)
{
    int e;
    frexp(x, &e);
    return e;
}


static float Sum1(float A, float f, int n)
{
    if (n < 0) { f = -f, n = -n; }

    if (n == 0) return A;

    while (1)
    {
        //  Record initial exponent.
        int e0 = E(A);

        A += f;
        if (--n == 0) return A;

        //  If exponent changed, continue to do another step.
        if (E(A) != e0) continue;

        /*  Note it is possible that A crossed zero here and is in the "same"
            binade with a different sign.  With rounding modes toward zero or
            round to nearest with ties to away, the rounding caused by adding f
            may differ.  However, if we have crossed zero, the next addition
            will put A in a new binade, and the loop will exit this iteration
            and start a new step.
        */

        float A0 = A;
        A += f;
        if (--n == 0) return A;

        //  If exponent changed, continue to do another step.
        if (E(A) != e0) continue;

        //  Calculate exact change from one addition.
        float g = A - A0;

        //  If g is zero, no future additions will change the sum.
        if (g == 0) return A;

        /*  Calculate how many additions can be done until just before the
            next binade.

            If A and f have the same sign, A is increasing in magnitude, and
            the next binade is at the next higher exponent.  Otherwise, the
            next binade is at the bottom of the current binade.

            We use copysign to get the right sign for the target.  .5 is used
            with ldexpf because C's specifications for frexp and ldexpf
            normalize the significand to [.5, 1) instead of IEEE-754's [1, 2).
        */
        int e1 = e0 + (signbit(A) == signbit(f));

        double t = floor((ldexpf(copysignf(.5, A), e1) - A) / (double) g);

        /*  Temporary workaround to avoid incorrectly steppping to the edge of
            a lower binade:
        */
        if (t == 0) continue;
        --t;

        if (n <= t) return A + n*g;

        A += t*g;
        n -= t;
    }
}


static int Check(float A, float f, int n)
{
    float A0 = Sum0(A, f, n);
    float A1 = Sum1(A, f, n);
    return A0 != A1;
}


static void Test(float A, float f, int n)
{
    float A0 = Sum0(A, f, n);
    float A1 = Sum1(A, f, n);
    if (A0 != A1)
    {
        printf("Adding %a = %.99g to %a = %.99g %d times:\n", f, f, A, A, n);
        printf("\tSum0 -> %a = %.99g.\n", A0, A0);
        printf("\tSum1 -> %a = %.99g.\n", A1, A1);
        exit(EXIT_FAILURE);
    }
}


int main(void)
{
    srand(1);

    for (int i = 0; i < 1000000; ++i)
    {
        if (i % 10000 == 0) printf("Tested %d cases.\n", i);
        Test(rand() - RAND_MAX/2, rand() - RAND_MAX/2, rand() % (1<<20));
    }
}