Skip to content
Snippets Groups Projects

Changes

Page history
  • New page
  • Templates
    • Edit Libraries/Proposals/ExpandFloating authored by Edward Kmett's avatar Edward Kmett
    Hide whitespace changes
    Inline Side-by-side
    libraries/proposals/expand-floating.md
    View page @ 58c219fa
    # 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](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
    ```wiki
    >>> exp 0.0000003
    1.000000300000045
    ......@@ -35,6 +46,7 @@ Consider
    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.
    ```wiki
    >>> exp 0.0000003 - 1
    3.0000004502817035e-7
    ......@@ -43,6 +55,7 @@ If we subtract out the leading term after it has destroyed all of our precision
    has lost a lot of precision relative to:
    ```wiki
    >>> expm1 0.0000003
    3.000000450000045e-7
    ......@@ -52,15 +65,21 @@ has lost a lot of precision relative to:
    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.
    ```
    classFractional a =>Floating a where...
    class Fractional a => Floating a where
    ...
    log1p :: a -> a
    expm1 :: a -> a
    log1pexp :: a -> a
    ......@@ -70,70 +89,105 @@ classFractional a =>Floating a where...
    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](https://groups.google.com/forum/#!msg/haskell-core-libraries/qXYMfV8JZ6k/tTuFrBMdDgAJ) around breaking builds. Consequently, at least at start we'll need to supply default definitions.
    ```
    log1p a = log (1+ a)expm1 a = exp a -1log1pexp a = log1p (exp a)log1mexp a = log1p (negate (exp a))
    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 #-}
    | 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
    | 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 #-}
    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.