log1p and expm1 are C standard library functions that are important for work with exponentials and logarithms.

I propose adding them to the Floating class where it is defined in GHC.Float.

We do not have to export these from Prelude. My knee-jerk reaction is that we probably shouldn't. The names are kind of awful, but are standard across the rest of the industry. We already have a precedent of not exporting clutter in the classes in the existing Data.Functor containing (<$), but it not currently coming into scope from the Prelude.

They are critical for any reasonably precise work with logarithms of values near 1, and of exponentials for small values of x, and it is somewhat embarrassing explaining to someone coming from the outside with a numerical background whom expects to find them to exist why we don't have them.

These arise all over the place in any work on probabilities in log-scale.

Backgrounder:

Consider

>>> exp 0.0000003
1.000000300000045

As the argument x get small, this gets very close to 1 + x. However, that leading 1 means you've consumed most of the precision of the floating point number you are using. 6 decimal places is ~18 bits of your significand that are just gone because of bad math.

If we subtract out the leading term after it has destroyed all of our precision it is too late.

>>> exp 0.0000003 - 1
3.0000004502817035e-7

has lost a lot of precision relative to:

>>> expm1 0.0000003
3.000000450000045e-7

Now every decimal place we get closer to 0 doesn't destroy a decimal place of precision!

Similar issues arise with logs of probabilities near 1. If you are forced to use log, as your probability gets closer to 1 from below you throw away most of your accuracy just by encoding the argument to the function with the same kind of error rate.

Here is straw man documentation ripped from my log-domain package that is probably way too technical, but serves as a starting point for discussion.

class ... => Floating a where

-- | The Taylor series for @'exp' x@ is given by
--
-- @
-- 'exp' x = 1 + x + x^2/2! + x^3/3! ...
-- @
--
-- When @x@ is small, the leading @1@ consumes virtually all of the available precision,

-- because subsequent terms are very small.
--
-- This computes:
--
-- @
-- exp x - 1 = x + x^2/2! + ..
-- @
--
-- For many types can afford you a great deal of additional precision if you move
-- things around algebraically to provide the 1 by other means.
expm1 :: Floating a => a -> a
expm1 x = exp x - 1

-- | Computes @log(1 + x)@
--
-- This is away from 0 so the Taylor series is defined, but it also provides an inverse to 'expm1'.
--
-- This can provide much more accurate answers for logarithms of numbers close to 1 (x near 0).
--
-- @
-- log1p (expm1 x) = log (1 + exp x - 1) = log (exp x)
-- @
log1p :: Floating a => a -> a
log1p x = log (1 + x)

They can be given definitions in terms of the standard C library functions for the CFloat, CDouble, Float and Double, either by foreign import or adding a pair of foreign prims.

Finally, here is a robust implementation for Data.Complex from the same package that properly deals with the subtleties involved in not losing precision.

expm1 x@(a :+ b)
| a*a + b*b < 1, u <- expm1 a, v <- sin (b/2), w <- -2*v*v = (u*w+u+w) :+ (u+1)*sin b
| otherwise = exp x - 1
{-# INLINE expm1 #-}

log1p x@(a :+ b)
| abs a < 0.5 && abs b < 0.5, u <- 2*a+a*a+b*b = log1p (u/(1+sqrt (u+1))) :+ atan2 (1 + a) b
| otherwise = log (1 + x)
{-# INLINE log1p #-}

Discussion Period: 2 weeks

-Edward Kmett