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

[Michael Snoyman]

+1 to the proposal as it stands.


On Thu, Apr 17, 2014 at 8:15 PM, Edward Kmett <ekmett at gmail.com> wrote:

> 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
>
> _______________________________________________
> Libraries mailing list
> Libraries at haskell.org
> http://www.haskell.org/mailman/listinfo/libraries
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.haskell.org/pipermail/libraries/attachments/20140418/0e6ef266/attachment.html>