TLW*_*TLW 3 floating-point precision ieee-754
是否有一种亚 O(n) 方法1可以在 IEEE 754 算术中复制浮点值的重复求和？
也就是说，是否有一种更快的方法可以精确复制以下伪代码的结果：
f64 notQuiteFMA(f64 acc, f64 f, bigint n) {
for (bigint i = 0; i < n; i++) {
acc = acc + f;
}
return acc;
}
？3
在实际算术中，答案很简单——return acc + f*n就足够了。然而，在有限精度下，答案要复杂得多。作为一个例子，请注意notQuiteFMA(0, 1, 1 << 256)大约是4 2^53，而不是人们期望的无限精度的 2^256。（这个例子还说明了为什么加快计算速度可能是可取的 - 计数到 2^256 是相当不切实际的。）
acc（鸽巢原理 - 在进入循环之前，重复加法不能进行超过 2^64 步 - 因此“只需”为和的每个初始值存储整个前缀和最终循环中的步数f。）x + 1 == x当加法不足以增加尾数并且结果向下舍入时，第一次发生在 2^53 左右。我相信确切的值取决于舍入模式。我们可以通过逐步执行 binades 来实现 O(log n ) 解决方案。这是因为二进制序列中的所有加法（可能除了第一个）都会将累积和更改相同的量，因此我们可以轻松地准确计算出重复加法直到二进制序列结束时将产生的总和。
\nbinade是某个整数e的区间 [2^ e , 2^ e +1 ) 或 (\xe2\x88\x922^ e +1 , \xe2\x88\x922 e ] ] 。就这个答案而言，整个最小指数范围（包括最小指数法线和所有次法线）将被视为单个二进制。
\n对于浮点数x（以浮点格式表示的数字），将 E( x ) 定义为x表示中使用的指数。鉴于此答案使用的二元数据的定义，E( x ) 表示x所在的二元数据，除了其符号。
\n代码格式的表达式表示浮点运算。x + y是x和y相加的实数算术结果。是和x+y相加的浮点结果。x和视情况使用相同的值。xyx
对于有限值，浮点加法是弱单调的：y<z意味着x+y<= x+z，在任何舍入模式下，除非x是无穷大并且yorz是相反符号的无穷大，导致 NaN。
令A 0、A 1和A 2A分别为 的两次执行之前、之间和之后的值A += f。考虑 E( A 0 ) \xe2\x89\xa5 E( A 1 ) = E( A 2 )的情况。令g为A 2 \xe2\x88\x92 A 1。g可以用 Sterbenz\xe2\x80\x99 引理表示，如果A 1是正规的，或者如果A 1是次正规的，因为A 1、A 2和f都是次正规的（或零）并且没有低于 ULP 的位A 2在加法中丢失。（由于g是可表示的，A2-A1因此计算它时没有舍入误差。）
我们将证明，所有后续添加都会产生E( ) = E( A 1 )f的值，增加相同的量g。AAA
令u为A 1的 ULP 。令r为 | 的余数 f| 对u取模。如果f为正，则令s为 +1；如果f为负，则令 s 为 +1。（ f为零时的情况是微不足道的。）使用定向舍入模式（朝向 \xe2\x88\x92\xe2\x88\x9e、朝向零或朝向 +\xe2\x88\x9e），任何舍入误差根据舍入方法和 的符号，后续不改变 E( ) 的加法要么始终为 \xe2\x88\x92 r ，要么始终为 u \xe2\x88\x92 r。使用舍入到最接近的关系，如果r < \xc2\xbd u ，则舍入误差始终为 \xe2\x88\x92 sr，否则始终为 s ( u \xe2\x88\x92 r )。舍入到最接近且与偶数相关时，如果r \xe2\x89\xa0 \xc2\xbd u，舍入误差始终为 \xe2\x88\x92 sr或始终s ( u \xe2\x88\x92 r )，如下与远方有联系。Af
这就留下了最接近的圆与r = \xc2\xbd u的联系。令b为f在u位置的位。在这种情况下， A 1的低位必须是偶数，因为它是f与A 0相加的结果， f 的余数为 \xc2\xbd u的u模（必然为 0 模u因为E( A 0 ) \xe2\x89\xa5 E( A 1 )) 产生A 1。除此之外，实数结果恰好位于两个可表示值的中间，因此它被四舍五入到偶数低位（2 u的倍数）。
\n鉴于A 1的低位是偶数，我们有两种情况：
\nb为 0，因此 的f加法结果为 0 u + \xc2\xbd u modulo 2 u \xe2\x89\xa1 \xc2\xbd u，并且四舍五入为偶数低位数，因此如果结果为正，则舍入误差为 \xe2\x88\x92\xc2\xbd u；如果结果为负，则舍入误差为+\xc2\xbd u 。A然后重复此操作以产生E( A) = E( A 1 )的后续加法。
b为 1，因此 的f加法结果为 1 u + \xc2\xbd u modulo 2 u \xe2\x89\xa1 1\xc2\xbd u，并且四舍五入为偶数低位如果结果为正，则舍入误差为 +\xc2\xbd u；如果结果为负，则舍入误差为 \xe2\x88\x92\xc2\xbd u。A然后重复此操作以产生E( A) = E( A 1 )的后续加法。
（要了解为什么我们在断言所有后续加法添加相同数量之前先进行二进制加法，请观察产生A 0的加法可能涉及来自先前A具有较低 ULP 的位，在这种情况下，实数-算术结果可能并不正好位于可表示值之间的中间位置，因此即使使用四舍五入到最接近的偶数，A 0也可能有奇数低位数字，并且A 2 \xe2\x88\x92 A 1不等于A 1 \xe2 \x88\x92 A 0。）
考虑到上述情况，我们可以计算出g的加法次数，比如说t可以执行到下一个二进制文件到达之前。然后我们可以在退出当前二进制文件A + tg之前计算A的准确值。然后，一次加法为我们提供了该二进制文件中的第一个A，另一个加法要么将我们带入一个新的二进制文件，要么为我们提供了对当前二进制文件执行计算所需的信息。然后我们可以对每个二进制重复此操作，直到达到所需的添加数量。
\n下面是演示这个概念的代码。根据评论，此代码需要一些改进。我希望在时间允许的情况下进行这项工作，这可能需要几周后的时间。代码中的函数E从 中获取指数frexp，而不是如上面讨论中所述将其限制为最小指数。这可能会减慢代码速度，但不会导致代码错误。我希望代码能够处理次正常情况。可以通过普通测试添加无穷大和 NaN。
t的计算（退出当前二进制之前可以执行的加法次数）可以使用更多关于所需精度以及舍入可能如何影响它的讨论。
\n#include <math.h>\n#include <stdio.h>\n#include <stdlib.h>\n\n\nstatic float Sum0(float A, float f, int n)\n{\n if (n < 0) { f = -f, n = -n; }\n while (n--) A += f;\n return A;\n}\n\n\n// Return the exponent of x, designating which binade it is in.\nstatic int E(float x)\n{\n int e;\n frexp(x, &e);\n return e;\n}\n\n\nstatic float Sum1(float A, float f, int n)\n{\n if (n < 0) { f = -f, n = -n; }\n\n if (n == 0) return A;\n\n while (1)\n {\n // Record initial exponent.\n int e0 = E(A);\n\n A += f;\n if (--n == 0) return A;\n\n // If exponent changed, continue to do another step.\n if (E(A) != e0) continue;\n\n /* Note it is possible that A crossed zero here and is in the "same"\n binade with a different sign. With rounding modes toward zero or\n round to nearest with ties to away, the rounding caused by adding f\n may differ. However, if we have crossed zero, the next addition\n will put A in a new binade, and the loop will exit this iteration\n and start a new step.\n */\n\n float A0 = A;\n A += f;\n if (--n == 0) return A;\n\n // If exponent changed, continue to do another step.\n if (E(A) != e0) continue;\n\n // Calculate exact change from one addition.\n float g = A - A0;\n\n // If g is zero, no future additions will change the sum.\n if (g == 0) return A;\n\n /* Calculate how many additions can be done until just before the\n next binade.\n\n If A and f have the same sign, A is increasing in magnitude, and\n the next binade is at the next higher exponent. Otherwise, the\n next binade is at the bottom of the current binade.\n\n We use copysign to get the right sign for the target. .5 is used\n with ldexpf because C\'s specifications for frexp and ldexpf\n normalize the significand to [.5, 1) instead of IEEE-754\'s [1, 2).\n */\n int e1 = e0 + (signbit(A) == signbit(f));\n\n double t = floor((ldexpf(copysignf(.5, A), e1) - A) / (double) g);\n\n /* Temporary workaround to avoid incorrectly steppping to the edge of\n a lower binade:\n */\n if (t == 0) continue;\n --t;\n\n if (n <= t) return A + n*g;\n\n A += t*g;\n n -= t;\n }\n}\n\n\nstatic int Check(float A, float f, int n)\n{\n float A0 = Sum0(A, f, n);\n float A1 = Sum1(A, f, n);\n return A0 != A1;\n}\n\n\nstatic void Test(float A, float f, int n)\n{\n float A0 = Sum0(A, f, n);\n float A1 = Sum1(A, f, n);\n if (A0 != A1)\n {\n printf("Adding %a = %.99g to %a = %.99g %d times:\\n", f, f, A, A, n);\n printf("\\tSum0 -> %a = %.99g.\\n", A0, A0);\n printf("\\tSum1 -> %a = %.99g.\\n", A1, A1);\n exit(EXIT_FAILURE);\n }\n}\n\n\nint main(void)\n{\n srand(1);\n\n for (int i = 0; i < 1000000; ++i)\n {\n if (i % 10000 == 0) printf("Tested %d cases.\\n", i);\n Test(rand() - RAND_MAX/2, rand() - RAND_MAX/2, rand() % (1<<20));\n }\n}\n