log1p.hs

log1p :: Double -> Double
log1p x
    | x == 0               = 0
    | x == -1              = -1 / 0
    | x < -1               = 0 / 0
    | x' < m_epsilon * 0.5 = x
    | (x > 0 && x < 1e-8) || (x > -1e-9 && x < 0)
                           = x * (1 - x * 0.5)
    | x' < 0.375           = x * (1 - x * chebyshevBroucke (x / 0.375) coeffs)
    | otherwise            = log (1 + x)
  where
    m_epsilon = encodeFloat (signif + 1) expo - 1.0
        where (signif, expo) = decodeFloat (1.0::Double)
    x' = abs x
    coeffs = [0.10378693562743769800686267719098e+1,
              -0.13364301504908918098766041553133e+0,
               0.19408249135520563357926199374750e-1,
              -0.30107551127535777690376537776592e-2,
               0.48694614797154850090456366509137e-3,
              -0.81054881893175356066809943008622e-4,
               0.13778847799559524782938251496059e-4,
              -0.23802210894358970251369992914935e-5,
               0.41640416213865183476391859901989e-6,
              -0.73595828378075994984266837031998e-7,
               0.13117611876241674949152294345011e-7,
              -0.23546709317742425136696092330175e-8,
               0.42522773276034997775638052962567e-9,
              -0.77190894134840796826108107493300e-10,
               0.14075746481359069909215356472191e-10,
              -0.25769072058024680627537078627584e-11,
               0.47342406666294421849154395005938e-12,
              -0.87249012674742641745301263292675e-13,
               0.16124614902740551465739833119115e-13,
              -0.29875652015665773006710792416815e-14,
               0.55480701209082887983041321697279e-15,
              -0.10324619158271569595141333961932e-15]
    chebyshevBroucke i = fini . foldr step [0, 0, 0]
        where
            step k [b0, b1, _] = [(k + i * 2 * b0 - b1), b0, b1]
            fini   [b0, _, b2] = (b0 - b2) * 0.5