我需要评估行的总和:1/1+1/2+1/3+...+1/n。考虑到在 C++ 中的评估并不完全准确,求和的顺序起着重要作用。1/n+1/(n-1)+...+1/2+1/1 表达式给出了更准确的结果。所以我需要找出求和的顺序,这提供了最大的准确性。我什至不知道从哪里开始。首选的实现语言是 C++。对不起,我的英语,如果有任何错误。
8 回答
对于较大的 n,您最好使用渐近公式,例如http://en.wikipedia.org/wiki/Harmonic_number上的公式;
另一种方法是使用 exp-log 转换。基本上:
H_n = 1 + 1/2 + 1/3 + ... + 1/n = log(exp(1 + 1/2 + 1/3 + ... + 1/n)) = log(exp(1) * exp(1/2) * exp(1/3) * ... * exp(1/n))。
您的标准库可以快速准确地计算指数和对数。使用乘法你应该得到更准确的结果。
如果这是您的作业并且您需要使用简单的加法,那么您最好按照其他人的建议从最小的加到最大的加法。
精度不足的原因是 float、double 和 long double 类型的精度。他们只存储这么多“小数”位。所以把一个很小的值加到一个很大的值上是没有效果的,小的项在较大的项中“丢失”了。
您要求和的系列有一条“长尾”,从某种意义上说,小项应该加起来有很大的贡献。但是如果你按降序求和,那么一段时间后每个新的小项都将不起作用(即使在此之前,它的大部分小数位都会被丢弃)。一旦达到这一点,您可以再添加十亿个术语,如果您一次执行一个,它仍然没有任何效果。
我认为按升序求和应该为此类系列提供最佳精度,尽管可能存在一些奇怪的极端情况,其中由于四舍五入(1/2)的幂而导致的错误可能恰好为某些人提供更接近的答案添加订单比其他订单。不过,您可能无法真正预测到这一点。
我什至不知道从哪里开始。
实际上,如果您要对大 N 进行求和,则从小到大的顺序相加并不是最好的方法——您仍然会遇到这样一种情况,即您相加的数字相对于总和而言太小而无法产生一个准确的结果。
以这种方式看待问题:无论排序如何,您都有 N 个求和,并且您希望总误差最小。因此,您应该能够通过最小化每个求和的误差来获得最小的总误差 - 并且通过添加尽可能接近彼此的值来最小化求和中的误差。我相信遵循这条逻辑链会给你一个部分和的二叉树:
Sum[0,i] = value[i]
Sum[1,i/2] = Sum[0,i] + Sum[0,i+1]
Sum[j+1,i/2] = Sum[j,i] + Sum[j,i+1]
依此类推,直到你得到一个答案。
当然,当 N 不是 2 的幂时,您最终会在每个阶段得到剩余部分,您需要将其结转到下一阶段的求和中。
(StackOverflow 的边距当然太小,无法证明这是最优的。部分原因是我没有花时间证明它。但它确实适用于任何 N,无论多么大,因为所有添加都是添加几乎相同大小的值。嗯,在最坏的非 2 次幂的情况下,除了 log(N) 之外的所有值,与 N 相比,这几乎是微乎其微的。)
http://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic 您可以找到具有现成 C/C++ 实现的库。
除非您使用一些准确的封闭形式表示,否则从小到大的有序求和可能是最准确的简单解决方案(我不清楚为什么 log-exp 会有所帮助 - 这是一个巧妙的技巧,但你不是据我所知,在这里赢得任何东西)。
通过意识到一段时间后,总和将变得“量化”,您可以进一步获得精度:实际上,当您有 2 位精度时,将 1.3 与 41 相加得到 42,而不是 42.3 - 但您通过保持几乎可以实现精度翻倍“错误”术语。这称为卡汉求和。您将计算误差项 (42-41-1.3 == -0.3) 并在下一次添加中通过在下一项添加 0.3 来纠正它,然后再将其添加。
除了从小到大的排序之外,Kahan Summation 可能会尽可能准确。我严重怀疑您是否需要更好的谐波级数 - 毕竟,即使经过 2^45 次迭代(疯狂的很多),您仍然只能处理至少 1/2^45 大的数字,并且一个数量级为 45 (<2^6) 的总和,数量级差异为 51 的二次幂 - 即,如果您添加“错误”的顺序,甚至仍可以用双精度变量表示。
如果你从小到大,并使用 Kahan Summation,那么太阳可能会在今天的处理器达到一定百分比的误差之前熄灭 - 而且你会遇到其他棘手的准确性问题,只是由于该规模的单个术语误差首先无论如何(无论如何,2^53 或更大的数量不能准确地表示为双精度数。)
我不确定求和的顺序是否起重要作用,我以前没听说过。我想您想在浮点运算中执行此操作,因此首先要考虑更多内联 (1.0/1.0 + 1.0/2.0+1.0/3.0) - 否则编译器将进行整数除法
确定评估顺序,可能是 for 循环或括号?
例如
float f = 0.0;
for (int i=n; i>0; --i)
{
f += 1.0/static_cast<float>(i);
}
哦忘了说,编译器通常会有开关来确定浮点评估模式。这可能与您按求和顺序所说的有关-在视觉 C+ 中,这些可以在代码生成编译设置中找到,在 g++ 中有选项 -float 可以处理此问题
实际上,另一个人是对的-您应该先按最小分量的顺序求和;所以 1/n + 1/(n-1) .. 1/1
这是因为浮点数的精度与比例相关,如果从 1 开始,相对于 1.0,您将拥有 23 位的精度。如果您从较小的数字开始,则精度相对于较小的数字,因此您将获得相对于 1xe-200 或其他任何值的 23 位精度。然后随着数字变大,会出现舍入误差,但总体误差会小于另一个方向
由于你所有的数字都是有理数,最简单的(也可能是最快的,因为它必须做更少的浮点运算)是用有理数(2个整数p,q的元组)进行计算,然后只做一个最后进行浮点除法。
更新以有效地使用此技术,您将需要对 p & q 使用 bigints,因为它们增长得非常快......
Lisp 中内置有理数的快速原型显示:
(defun sum_harmonic (n acc)
(if (= n 0) acc (sum_harmonic (- n 1) (+ acc (/ 1 n)))))
(sum_harmonic 10 0)
7381/2520
[2.9289682]
(sum_harmonic 100 0)
14466636279520351160221518043104131447711/278881500918849908658135235741249214272
[5.1873775]
(sum_harmonic 1000 0)
53362913282294785045591045624042980409652472280384260097101349248456268889497101
75750609790198503569140908873155046809837844217211788500946430234432656602250210
02784256328520814055449412104425101426727702947747127089179639677796104532246924
26866468888281582071984897105110796873249319155529397017508931564519976085734473
01418328401172441228064907430770373668317005580029365923508858936023528585280816
0759574737836655413175508131522517/712886527466509305316638415571427292066835886
18858930404520019911543240875811114994764441519138715869117178170195752565129802
64067621009251465871004305131072686268143200196609974862745937188343705015434452
52373974529896314567498212823695623282379401106880926231770886197954079124775455
80493264757378299233527517967352480424636380511370343312147817468508784534856780
21888075373249921995672056932029099390891687487672697950931603520000
[7.485471]
因此,下一个更好的选择可能是保留浮点列表并减少它,将每个步骤中的两个最小数字相加......