大数与小数的求和算法
原文首发于我的微信公众号:GeekArtT .
在计算机求和的过程中,一个大数和小数的相加会因为浮点数的有限精度,而导致截断误差的出现。所以在构建计算网格的时候,都要极力避免这样情形的发生,将计算统一在相对较近的数量级上。所以,当需要对一系列的数值做加法时,一个好的技巧是将这些数由大到小做排列,再逐个相加。
而如果一定要做出这样的大数与小数的求和,一个直观想法就是:大数部分和小数部分的高位相加,将剩余的小数部分作为单独的“补全”部分相加。这种直观想法的官方名称叫做Kahan求和法。
假设当前的浮点数变量可以保存6位的数值。那么,数值12345与1.234相加的理论值应该是12346.234。但由于当前只能保存6位数值,这个正确的理论值会被截断为12346.2,这就出现了0.034的误差。当有很多这样的大数与小数相加时,截断误差就会逐步累积,导致最后的计算结果出现大的偏差。
Kahan算法:
def KahanSum(input): var sum = 0.0 var c = 0.0 for i = 1 to input.length do var y = input[i] - c // Initially, c is zero; then it compensates previous accuracy. var t = sum + y // low-order digits of y are lost c = (t - sum) - y // recover the low-order digits of y, with negative symbol sum = t next i return sum
在上述伪代码中,变量c表示的即是小数的补全部分compensation,更严格地说,应该是负的补全部分。随着这个补全部分的不断积累,当这些截断误差积累到一定量级,它们在求和的时候也就不会被截断了,从而能够相对好地控制整个求和过程的精度。
以下,先用一个具体的理论例子来说明。比如,用10000.0 + pi + e来说明,我们依旧假设浮点型变量只能保存6位数值。此时,具体写出求和算式应该是:10000.0 + 3.14159 + 2.71828,它们的理论结果应该是10005.85987,约等于10005.9。
但由于截断误差,第一次求和10000.0 + 3.14159只能得到结果10003.1;这个结果再与2.71828相加,得到10005.81828,被截断为10005.8。此时结果就相差了0.1。
运用Kahan求和法,我们的运行过程是(记住,我们的浮点型变量保存6位数值),
第一次求和:
y = 3.14159 - 0.00000 t = 10000.0 + 3.14159 = 10003.14159 = 10003.1 // low-order digits have lost c = (10003.1 - 10000.0) - 3.14159 = 3.10000 - 3.14159 = - (.0415900) // recover the negative parts of compensation errors sum = 1003.1
第二次求和:
y = 2.71828 - (-.0415900) = 2.75985 // Add previous compensated parts to current small number t = 10003.1 + 2.75987 = 10005.85987 = 10005.9 // As the low-order digits have been accumulated large enought, it won't be canceled by big number c = (10005.9 - 10003.1) - 2.75987 = 2.80000 - 2.75987 = .040130 sum = 10005.9
以上是理论分析。再举一个可以运行的Python代码示例,方便感兴趣的朋友做研究。这个例子曾经出现于Google的首席科学家Vincent Vanhoucke在Udacity上开设的deep learning课程。这个求和算式是:在10^9的基础上,加上10^(-6),重复10^6次,再减去10^9,即10^9 + 10^6*10^(-6) - 10^9,理论值应该为1。
Python Code:
summ = 1000000000 for indx in xrange(1000000): summ += 0.000001 summ -= 1000000000 print summ
运行后,可以得到结果是0.953674316406。可以看到,在10^6次求和后,截断误差的累积量已经非常可观了。
如果我们用Kahan求和法来做改进,可以得到:
Python Code with Kahan method:
summ = 1000000000 c = 0.0 for indx in xrange(1000000): y = 0.000001 - c t = summ + y c = (t - summ) - y summ = t summ -= 1000000000 print summ
运行后,我们可以欣喜地看到正确结果:1.0。
如果你喜欢我的文章或分享,请长按下面的二维码关注我的公众号,谢谢!