在Haskell中使用系列的自然对数给出了不精确的结果

Question

所以我编写了两个函数来计算变量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的日志函数)正确计算自然对数？

Answer 1

浮点数学很棘手.导致精度损失的一个原因是增加了数量级别不同的数字.例如,在你的算法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)表现都是一样的数字.

但总的来说,我建议阅读一本数值分析书,以了解更多有关此类问题的信息.