Natural Logarithm using series in Haskell Gives Imprecise Results

Question

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?

(1) Floating-point arithmetic is imprecise by nature. (2) Part of the imprecision in your functions arises merely because you are using `Float` rather than `Double`. — duplode, Apr 09 '17 at 04:16
@duplode thanks for the response. It's indeed more precise, but after I changed the functions type signature from float to Double, then made a comparison to Prelude log, the last few digits still differ vastly. — Leo Mingo, Apr 09 '17 at 04:23

hugomg · Accepted Answer · 2017-04-09T17:38:27.710

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.

Natural Logarithm using series in Haskell Gives Imprecise Results

1 Answers1