64 位机器上的 Matlab 中表示为双精度的无理数集合 N 的大小的严格下限是多少，我在乘积的 k 个十进制数字时相乘？ 例如，在编码不同的随机 pi 块的 ~10^12 双精度相乘后，我可以期望什么精度？
如果要求严格限制，则如果所有操作都以 IEEE 754 双精度执行，则 @EricPostpischil 的响应是绝对误差边界。
如果你要求信心，我把它理解为错误的统计分布。假设 [-e/2，e/2] 中的误差分布均匀，您可以尝试在数学堆栈交换的 M 操作后询问误差的理论分布......我想紧绷的界限在某种程度上是非常保守的。
让我们用一些 Smalltalk 代码来说明对这些统计数据的实验估计(任何具有大整数/分数算术的语言都可以做到(：
nOp := 500.
relativeErrorBound := ((1 + (Float epsilon asFraction / 2)) raisedTo: nOp * 2 - 1) - 1.0.
nToss := 1000.
stats := (1 to: nToss)
collect: [:void |
| fractions exactProduct floatProduct relativeError |
fractions := (1 to: nOp) collect: [:e | 10000 atRandom / 3137].
exactProduct := fractions inject: 1 into: [:prod :element | prod * element].
floatProduct := fractions inject: 1.0 into: [:prod :element | prod * element].
relativeError := (floatProduct asFraction - exactProduct) / exactProduct.
relativeError].
s1 := stats detectSum: [:each | each].
s2 := stats detectSum: [:each | each squared].
maxEncounteredError := (stats detectMax: [:each | each abs]) abs asFloat.
estimatedMean := (s1 /nToss) asFloat.
estimatedStd := (s2 / (nToss-1) - (s1/nToss) squared) sqrt.
我得到 nOp=20 双倍乘法的这些结果：
relativeErrorBound -> 4.440892098500626e-15
maxEncounteredError -> 1.250926201710214e-15
estimatedMean -> -1.0984634797115124e-18
estimatedStd -> 2.9607828266493842e-16
对于 nOp=100：
relativeErrorBound -> 2.220446049250313e-14
maxEncounteredError -> 2.1454964094158273e-15
estimatedMean -> -1.8768492273800676e-17
estimatedStd -> 6.529482793500846e-16
对于 nOp=500：
relativeErrorBound -> 1.1102230246251565e-13
maxEncounteredError -> 4.550696454362764e-15
estimatedMean -> 9.51007740905571e-17
estimatedStd -> 1.4766176010100097e-15
您可以观察到标准差的增长比误差边界的增长慢得多。
更新：在第一次近似(1+e)^m = 1+m*e+O((m*e)^2)，所以只要 m*e 足够小，分布在 [-e，e] 中大约是 m 个均匀的和，并且这个和非常接近方差m*(2e)^2/12的正态分布(高斯(。您可以在 Matlab 中检查std(sum(rand(100,5000)))是否靠近sqrt(100/12)。
我们可以认为它对于m=2*10^12-1仍然是正确的，即大约m=2^41，m*e=2^-12。在这种情况下，全局误差是准正态分布，全局误差的标准差是sigma=(2^-52*sqrt(2^41/12))或近似sigma=10^-10
请参阅计算 P(abs(error(>k*sigma 的 http://en.wikipedia.org/wiki/Normal_distribution(
- 在 68% 的情况下(1 西格玛(，您将拥有 10 位或更高的精度。
-
erfc(10/sqrt(2))给出了精度小于 9 位的概率，大约 6*10^22 中的 1 种情况，所以我让你计算只有 4 位精度的概率(你不能用双精度评估它，它会下溢(!!
比理论标准差小一点(2e-15 9e-16 4e-16 表示 20 100 和 500 双倍(，但这一定是由于我的输入误差 i/3137 i=1..10000
...这是提醒结果的好方法，如果错误是由浮点运算(如 M_PI*num/den (引起的，则结果将取决于输入中的错误分布，则误差可能会超过 e。
此外，正如 Eric 所说，仅使用 * 是一个非常理想的情况，如果您混合 +，事情可能会退化得更快。
最后注意：我们可以制作一个达到最大误差界限的输入列表，将所有元素设置为 (1+e(，四舍五入为 1.0，我们得到最大理论误差界限，但是我们的输入分布非常偏颇！下摆错了，因为所有的乘法都是精确的，我们只得到 (1+e(^n，而不是 (1+e(^(2n-1(，所以大约只有一半的误差......
更新2：逆问题
既然你想要相反的，序列的长度 n 是多少，这样我就可以得到 k 位精度，具有一定的置信度 10^-c
我只会回答k>=8，因为在上面的近似中需要(m*e) << 1。
让我们取c=7，你得到 k 位数字的置信度为 10^-7 意味着5.3*sigma < 10^-k.
sigma = 2*e*sqrt((2*n-1)/12)这与e=2^-53 n=0.5+1.5*(sigma/e)^2.
因此 n ~ 3*2^105*sigma^2，作为 sigma^2 <10^-2k/5.3^2，我们可以写n < 3*2^105*10^-(2*k)/5.3^2
A.N. 对于长度 n=4.3e12，少于 k=9 位的概率小于 10^-7，对于 10 位数字，大约 n=4.3e10。
对于 15 位数字，我们会达到 n=4 个数字，但这里的正态分布假设非常粗糙并且不成立，尤其是 5 西格玛处的分布尾部，因此请谨慎使用(贝里-埃森定理限制了这种分布离正态有多远 http://en.wikipedia.org/wiki/Berry-Esseen_theorem(
假设所有输入、中间值和最终值不下溢或溢出，所描述的 M 操作中的相对误差最多为 (1+2-53(M-1。
考虑将实数 a0 转换为双精度。结果是某个数字 a0•(1+e(，其中 -2-53 ≤ e ≤ 2-53(因为转换为双精度应始终产生最接近的可表示值，而双精度值的量程是最高位的 2-53，最接近的值始终在半个量子内(。为了进一步分析，我们将考虑 e 的最坏情况值，2-53。
当我们将一个(先前转换的(值乘以另一个值时，数学上精确的结果是 a0•(1+e( • a1•(1+e(。计算结果有另一个舍入误差，因此计算结果为 a0•(1+e( • a1•(1+e( • (1+e( = a0 • a1 • (1 +e(3。显然，这是 (1+e(3 的相对误差。对于这些操作，我们可以看到误差累积为 (1+e(M：每个运算将所有先前的误差项乘以 1+e。
给定 N 个输入，将有 N 次转换和 N-1 次乘法，因此最坏的误差将是 (1+e(2 N - 1。
此错误的相等性仅在 N≤1 上实现。否则，误差必须小于此界限。
请注意，如此简单的错误边界仅在简单问题中才有可能，例如具有齐次运算的问题。在典型的浮点运算中，混合了加法、减法、乘法和其他运算，通常不可能如此简单地计算界限。
对于 N=10 12 (M=2•10 12-1(，上述界限小于 2.000222062•1012 单位的 2-53，并且小于 .0002220693。因此，计算结果对四个十进制数字以下的东西很好。(但请记住，您需要避免溢出和下溢。
(注意上述计算的严格性：我使用 Maple 精确计算二项式 (1+2-53(2•1012-1 的 1000 项(删除了最初的 1 项(，并添加一个可证明大于所有剩余项之和的值。然后我让 Maple 将确切的结果评估为 1000 个十进制数字，它小于我上面报告的界限。
对于 64 位浮点数，假设标准 IEEE 754 具有 52+1 位
尾数。这意味着相对精度介于 1.0000...0 和 1.0000...1 之间，其中小数点后的二进制位数为 52。(您可以将 1.000...0 视为存储在尾数 AKA 有效数中的二进制内容(。
误差是 52 次方的 1/2 除以 2(分辨率的一半(。注意我选择尽可能接近 1.0 的相对精度，因为它是最坏的情况(否则在 1.111..11 和 1.111..01 之间，它更精确(。
在十进制中，双精度的最坏情况相对精度为 1.11E-16。
如果将 N 乘以此精度加倍，则新的相对精度(假设中间舍入没有额外的误差(为：
1 - (1 - 1.11E-16)^N
因此，如果将 pi(或任何双精度 10^12(乘以，则误差上限为：
1.1102e-004
这意味着您可以对大约 4-5 位数字充满信心。
如果您的 CPU 支持中间结果的扩展精度浮点数，则可以忽略中间舍入误差。
如果未使用扩展精度 FPU(浮点单位(，则中间步骤中的舍入会引入额外的误差(与乘法相同(。这意味着严格的下限计算公式为：
1 -
((1 - 1.11E-16) * (1 - 1.11E-16) * (1 - 1.11E-16)
* (1 - 1.11E-16) * (1 - 1.11E-16) % for multiplication, then rounding
... (another N-4 lines here) ...
* (1 - 1.11E-16) * (1 - 1.11E-16))
= 1-(1-1.11E-16)^(N*2-1)
如果 N 太大，则运行时间太长。可能的误差(使用中间舍入(为 2.2204e-012，与没有中间舍入1-(1 - 1.11E-16)^N =1.1102e-012 相比是两倍。
大约，我们可以说中间舍入会使误差加倍。
如果将 pi 乘以 10^12 倍，并且没有扩展精度 FPU。这可能是因为您在继续之前将中间步骤写入内存(也许还会做其他事情((只需确保编译器没有对您的指令重新排序，以便没有 FPU 结果累积(，那么相对错误的严格上限是：
2.22e-004
的置信度并不意味着它有时完全是小数位。
例如，如果答案是：
1.999999999999，错误为 1E-5，实际答案可能是 2.000001234。
在这种情况下，即使是第一个十进制数字也是错误的。但这实际上取决于你有多幸运(答案是否落在这样的边界上(。
此解决方案假定双精度(包括答案(都已规范化。对于非规范化结果，显然，非规范化所用的二进制位数会将准确性降低这么多位。