Proposal: Add log1p and expm1 to GHC.Float.Floating

[Edward Kmett]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.haskell.org/pipermail/libraries/attachments/20140417/aa3c7d54/attachment.html>