在Haskell中使用系列的自然对数给出不精确的结果(Natural Logarithm using series in Haskell Gives Imprecise Results)

在Haskell中使用系列的自然对数给出不精确的结果(Natural Logarithm using series in Haskell Gives Imprecise Results)

所以我写了两个函数来计算变量x的自然对数，将递增和的上限增加到33000，函数仍然返回在ghci中测试的不精确结果，与从Prelude导入的缺省日志函数相比，这里是代码定义：

lnOfx :: Float -> Float lnOfx x = netSum f 1 33000 where f i = 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm where oddTerm = 2*i-1 netSum f minB maxB = sum [f i | i <- [minB .. maxB]] lnOfx2 :: Float -> Float lnOfx2 x = netSum f 1 33000 where f i = (1/i)*((x-1)/x)**i netSum f minB maxB = sum [f i | i <- [minB .. maxB]]

测试结果如下：

log 3 1.0986122886681098 lnOfx 3 1.0986125 lnOfx2 3 1.0986122 log 2 0.6931471805599453 lnOfx 2 0.6931472 lnOfx2 2 0.6931473

那么，为什么结果会不一样，以及正确的方法（如Prelude的日志函数）能够正确计算自然对数？

So I wrote two functions to calculate natural logarithm of variable x, after increasing the upperbound of incremental sum to 33000 the functions still return imprecise results being tested in ghci, compared with the default log function imported from Prelude, here is the code definition:

lnOfx :: Float -> Float lnOfx x = netSum f 1 33000 where f i = 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm where oddTerm = 2*i-1 netSum f minB maxB = sum [f i | i <- [minB .. maxB]] lnOfx2 :: Float -> Float lnOfx2 x = netSum f 1 33000 where f i = (1/i)*((x-1)/x)**i netSum f minB maxB = sum [f i | i <- [minB .. maxB]]

And the test results:

log 3 1.0986122886681098 lnOfx 3 1.0986125 lnOfx2 3 1.0986122 log 2 0.6931471805599453 lnOfx 2 0.6931472 lnOfx2 2 0.6931473

So why do the results differ and what is the right way(like the log function from Prelude does) to calculate natural logarithm correctly?

最满意答案

浮点数学是棘手的。 可能导致精确度损失的一个因素是添加数量大小不同的数字。 例如，在你的算法中，从i=25开始，你的总和中的条款足够小，以至于他们不再有所作为：

-- 25t term: let f x i = let oddTerm = 2*i-1 in 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm let y = f 3.0 25 -- summation up to 24 item let s = 1.098612288668109 -- this will return True, surprisingly! s + y == s

你可以做的一件事是减少这种情况是以相反的顺序添加数字，所以在添加到大数字之前，小数字会加在一起。

lnOfx :: Float -> Float lnOfx x = netSum f 1 33000 where f i = 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm where oddTerm = 2*i-1 netSum f minB maxB = sum $ reverse [f i | i <- [minB .. maxB]]

在我的测试中，这足以让print (lnOfx 3.0)和print (log 3.0)显示所有相同的数字。

但总的来说，我会推荐阅读一本数值分析书来了解更多关于这类问题的内容。

Floating point math is tricky. One of the thing that can cause loss of precision is adding numbers with very different magnitudes. For example, in your algorithm starting around i=25 the terms in your sum are small enough that they stop making a difference:

-- 25t term: let f x i = let oddTerm = 2*i-1 in 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm let y = f 3.0 25 -- summation up to 24 item let s = 1.098612288668109 -- this will return True, surprisingly! s + y == s

One thing you can do to mitigate this is add the numbers in the reverse order, so the small numbers get added together before they are added to the big numbers.

lnOfx :: Float -> Float lnOfx x = netSum f 1 33000 where f i = 2*(1/oddTerm)*((x-1)/(x+1))**oddTerm where oddTerm = 2*i-1 netSum f minB maxB = sum $ reverse [f i | i <- [minB .. maxB]]

In my tests this was enough so that print (lnOfx 3.0) and print (log 3.0) showed all the same digits.

But in general I would recommend reading a numerical analysis book to learn more about this kind of problem.