1

无理数集N的大小的严格下限是多少，在 64 位机器上的 Matlab 中表示为双精度数，我在对乘积的k个十进制数字有信心的情况下将其相乘？例如，在将编码不同随机 pi 块的 ~10^12 双打相乘后，我可以期望什么精度？

4

3 回答 3

3

如果您要求严格限制，如果所有操作都以 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 足够小，分布大约是 m 在 [-e,e] 中均匀分布的总和，并且这个总和非常接近方差的正态分布（高斯）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

请参阅http://en.wikipedia.org/wiki/Normal_distribution以计算 P(abs(error)>k*sigma)

  • 在 68% 的情况下（1 sigma），您将拥有 10 位或更多位的精度。
  • erfc(10/sqrt(2))给你小于 9 位精度的概率，大约 6*10^22 中的 1 个案例，所以我让你计算只有 4 位精度的概率（你不能用双精度来评估它，它下溢）！！！

我的实验标准偏差比理论偏差小一点（2e-15 9e-16 4e-16 for 20 100 & 500 double）但这一定是由于我的输入错误 i/3137 i=1..10000 的偏差分布...

这是提醒结果将由输入中的错误分布主导的好方法，如果它们是由诸如 M_PI*num/den 之类的浮点运算产生的，则可能会超过 e

此外，正如 Eric 所说，仅使用 * 是一个非常理想的情况，如果混合 +，事情可能会退化得更快。

最后一点：我们可以制作一个达到最大误差界限的输入列表，将所有元素设置为 (1+e) 将四舍五入为 1.0，我们得到最大理论误差界限，但我们的输入分布有很大偏差！HEM 错误，因为所有乘法都是精确的，我们只得到 (1+e)^n，而不是 (1+e)^(2n-1)，所以大约只有一半的错误......

更新2：逆问题

既然你想要倒数，那么序列的长度 n 是多少，这样我就可以得到 k 位精度并具有一定的置信度 10^-c

我将只回答k>=8，因为(m*e) << 1在上述近似值中是必需的。

假设c=7，您得到 k 位，置信度为 10^ -75.3*sigma < 10^-k.
sigma = 2*e*sqrt((2*n-1)/12)意味着n=0.5+1.5*(sigma/e)^2。 因此 n ~ 3*2^105*sigma^2，因为 sigma^2 < 10^-2k/5.3^2，我们可以写e=2^-53
n < 3*2^105*10^-(2*k)/5.3^2

AN 对于长度 n=4.3e12，少于 k=9 位的概率小于 10^-7，对于 10 位，n=4.3e10 左右。

对于 15 位，我们将达到 n=4 的数字，但这里我们的正态分布假设非常粗略且不成立，尤其是 5 sigma 处的分布尾部，因此请谨慎使用（Berry-Esseen theorem bounds how far from normal is such distribution http ://en.wikipedia.org/wiki/Berry-Esseen_theorem）

于 2012-08-21T09:54:11.537 回答
2

假设所有输入值、中间值和最终值都没有下溢或上溢，所描述的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•10 12 个单位的 2 -53，并且小于 0.0002220693。所以计算结果对四位小数以下的东西是好的。（但请记住，您需要避免上溢和下溢。）

（注意上述计算的严格性：我使用 Maple 精确计算了二项式 (1+2 -53 ) 2•10 12 -1的 1000 项（已删除最初的 1 项）并添加一个可证明的值大于所有剩余项的总和。然后我让 Maple 将该精确结果评估为 1000 个十进制数字，它小于我上面报告的界限。）

于 2012-08-21T03:28:43.083 回答
1

对于 64 位浮点数，假设标准 IEEE 754，具有 52+1 位尾数。

这意味着相对精度在 1.0000...0 和 1.0000...1 之间，其中小数点后的二进制位数为 52。（您可以将 1.000...0 视为以二进制形式存储在尾数 AKA 有效位）。

误差是 1/2 的 52 次方除以 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 太大，则运行时间太长。1-(1 - 1.11E-16)^N可能的误差（中间舍入）为 2.2204e-012，与没有中间舍入= 1.1102e-012相比是两倍。

大约，我们可以说中间舍入使误差加倍。

如果将 pi 乘以 10^12 次，并且没有扩展精度 FPU。这可能是因为您在继续之前将中间步骤写入内存（并且可能执行其他操作）（只需确保编译器没有重新排序您的指令，以便没有 FPU 结果累积），然后对您的亲戚有一个严格的上限错误是：

2.22e-004

请注意，对小数位的信心并不意味着它有时会完全是小数位。

例如，如果答案是：

1.999999999999，错误是1E-5，实际答案可能是2.000001234。

在这种情况下，即使是第一个十进制数字也是错误的。但这真的取决于你有多幸运（答案是否落在这样的边界上）。

该解决方案假定双打（包括答案​​）都已标准化。对于非规范化的结果，显然，对其进行非规范化的二进制数字会使精度降低那么多位数。

于 2012-08-21T02:07:00.620 回答