Expand Floating Proposal

THIS DESCRIPTION IS STILL WORK IN PROGRESS

This was originally proposed in https://mail.haskell.org/pipermail/libraries/2014-April/022667.html for inclusion into GHC 7.8.

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

The proposal is to add them to the Floating class where it is defined in GHC.Float.

We do not have to export these from Prelude. 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.

Final API

The final API (suitable for Haskell Report inclusion) we want to end up with 4 additional members to the Floating class.


class Fractional a => Floating a where
    ...
    log1p :: a -> a
    expm1 :: a -> a
    log1pexp :: a -> a
    log1mexp :: a -> a

The original proposal included defaults, but in the original libraries mailing list discussion it was decided that having no defaults would be a better end state. However, it isn't possible to proceed without defaults and comply with a 3 release policy around breaking builds. Consequently, at least at start we'll need to supply default definitions.

log1p a = log (1 + a)
expm1 a = exp a - 1  
log1pexp a = log1p (exp a)
log1mexp a = log1p (negate (exp a))

Migration plan

Phase 1 (GHC 8.0)

Add log1p, expm1, log1pexp, log1mexp to Floating in GHC.Float.
Do not export these from the Prelude. However, add a re-export of Floating with all of its members from the existing Numeric module.
Include the naive default definitions in the Floating class:
Add smart overrides for instance Floating Float and instance Floating Double:

We can import efficient, accurate implementations from libm and the like for log1p and expm1. They may be available for log1pexp and log1mexp as well. If not:

  log1mexp a
    | a <= log 2 = log (negate (expm1 a))
    | otherwise  = log1p (negate (exp a))
  {-# INLINE log1mexp #-}
  log1pexp a
    | a <= 18   = log1p (exp a)
    | a <= 100  = a + exp (negate a)
    | otherwise = a

will serve for both Float and Double, combining multiple locally accurate approximations.

Add a smart override for instance Floating a => Floating (Complex a):

  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 #-}

We've decided to move ahead with Phase 1 and defer decisions on further phases. Further actions are open for discussion. These include working out a roadmap for removing the default definitions and/or moving the definitions into the Prelude.

Removing Defaults:

Phase 2

Implement a warning that we're going to remove the default definitions of log1p and expm1 (and possibly log1pexp and log1mexp.) Add this warning to -Wcompat.

Phase 3 (GHC 8.4 or later)

Move the warning from -Wcompat to -Wall one release prior to making the change. (simply expanding the MINIMAL pragma of Floating to include log1p and expm1 would do.)

Phase 4 (GHC 8.6 or later)

Remove the default definitions of log1p, expm1 (and possibly log1pexp and log1mexp).

Moving into the Prelude

Phase 2 (GHC 8.4 or later)

Warn about additional members coming into the Prelude

Phase 3 (GHC 8.6 or later)

Expose the additional members of Floating from the Prelude.

Dependencies

If we decide to proceed with both of these refinements of the original plan then Moving into the Prelude Phase 3 should likely coincide with Removing Defaults Phase 3 or 4.

expand floating