计算谐波级数部分和时累积的计算错误

Question

我目前使用两种方法来计算调和级数的部分和:

1. 暴力破解(迭代)
2. 使用某些公式可以提供解决问题的捷径

困扰我的问题是:double 精确到小数点后 15 位，而 long double 精确到小数点后 19 位。我正在处理从 5 到大约 50 的部分和。例如，要达到 5 的总和，您需要调和级数的大约 83 个成员，而要达到 50 的总和，您需要大约 2911068851283175669760 个成员谐波系列。我目前使用的是 double 。

计算大总和(例如 >20)时，使用迭代和公式计算此类总和时是否会失去准确性？

如果是这样，使用一些高精度库会有好处吗？

Answer 1

使用任何为(无界)大整数和分数提供精确算术的语言，我们可以尝试使用非常简单的代码轻松评估舍入误差，至少对于小 n(精确算术可能会变得令人望而却步)

我在这里随意选择 Squeak Smalltalk 来说明这一点。

首先，让我们定义一个闭包，它将返回精确的部分和

hn := [:n | (1 to: n) detectSum: [:i | i reciprocal]].

然后另一个将通过迭代返回 Float( double )，从较小的项开始:

fn := [:n | (n to: 1 by: -1) detectSum: [:i | i asFloat reciprocal]].

现在另一个闭包将根据 ulp 估计误差:

un := [:n | | f | ((f := fn value: n) asFraction - (hn value: n)) abs / f ulp].

现在，让我们搜索前 1000 个部分和的最坏情况:

i := (1 to: 1000) detectMax: [:n | un value: n ].
{i. un value: i}.

结果为#(807 6.101840327204507)，H807的ulp错误约为6。这并不像您所看到的那样。

如果我们尝试先用更大的谐波求和:

fn := [:n | (1 to: n) detectSum: [:i | i asFloat reciprocal]].

然后结果是#(471 5.714808765246509)，5.7 ulp，惊喜一点。

我们可以尝试用后面的 fn 看得更远:

un value: 10000.
-> 1.8242087614993667

使用内存，我发现对于 100,000 个第一项 (H ~ 7)，最大误差约为 54 ulp，如果从较小的项开始更仔细地求和，则最大误差约为 26 ulp。

因此，最多有几千个术语，没有用从较小的术语开始求和(如果您要扫描增加的术语，则可以进行内存)，并且也没有必要使用卡汉求和，除非您真的关心最后几个术语ulp。

Kahan 对于数百万到数十亿的项会变得有趣，但对于更大的总和，如果你有 Euler Mascheroni 的良好近似，公式将会更快，也更准确，自然对数永远不会超过几个 ulp 。但你应该展示你的代码，否则这只是猜测。

例如，这个简单的二阶公式

fn := [:n | n ln + (0.5772156649015328606 + ((n*2.0) reciprocal - (n squared*12.0) reciprocal))].

对于较小的 n 值，将给出大约 5.5 ulp 的错误，但在使用像样的 MacOs libm 时，100,000 以上的 ulp 误差小于 1.5 ulp。