Search code examples
haskellmathlibm

Natural Logarithm using series in Haskell Gives Imprecise Results


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?


Solution

  • 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.