Compensated.hs

{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE PatternGuards #-}
{-# LANGUAGE Trustworthy #-}

#ifndef MIN_VERSION_vector
#define MIN_VERSION_vector(x,y,z) 1
#endif

--------------------------------------------------------------------
-- |
-- Copyright :  (c) Edward Kmett 2013-2015
-- License   :  BSD3
-- Maintainer:  Edward Kmett <ekmett@gmail.com>
-- Stability :  experimental
-- Portability: non-portable
--
-- This module provides a fairly extensive API for compensated
-- floating point arithmetic based on Knuth's error free
-- transformation, various algorithms by Ogita, Rump and Oishi,
-- Hida, Li and Bailey, Kahan summation, etc. with custom compensated
-- arithmetic circuits to do multiplication, division, etc. of compensated
-- numbers.
--
-- In general if @a@ has x bits of significand, @Compensated a@ gives
-- you twice that. You can iterate this construction for arbitrary
-- precision.
--
-- References:
--
-- * <http://web.mit.edu/tabbott/Public/quaddouble-debian/qd-2.3.4-old/docs/qd.pdf>
--
-- * <http://www.ti3.tuhh.de/paper/rump/OgRuOi05.pdf>
--
-- * Donald Knuth's \"The Art of Computer Programming, Volume 2: Seminumerical Algorithms\"
--
-- * <http://en.wikipedia.org/wiki/Kahan_summation_algorithm>
--------------------------------------------------------------------
module Numeric.Compensated
  ( Compensable(..)
  , _Compensated
  , Overcompensated
  , primal
  , residual
  , uncompensated
  , fadd
  -- * lifting scalars
  , add, times, squared, divide, split
  , kahan, (+^), (*^)
  -- * compensated operators
  , square
  ) where

#if __GLASGOW_HASKELL__ < 710
import Control.Applicative
#endif
import Control.Lens as L
import Control.DeepSeq
import Data.Binary as Binary
import Data.Data
import Data.Foldable as Foldable
import Data.Function (on)
import Data.Hashable
import Data.Ratio
import Data.SafeCopy
import Data.Serialize as Serialize
import Data.Bytes.Serial as Bytes
#if !(MIN_VERSION_base(4,11,0))
import Data.Semigroup
#endif
import Data.Vector.Unboxed as U
import Data.Vector.Generic as G
import Data.Vector.Generic.Mutable as M
import Foreign.Ptr
import Foreign.Storable
import Numeric (Floating(..))
import Text.Read as T
import Text.Show as T

-- $setup
-- >>> import Numeric.Compensated

{-# ANN module "hlint: ignore Use -" #-}
{-# ANN module "hlint: ignore Use curry" #-}
{-# ANN module "hlint: ignore Eta reduce" #-}
{-# ANN module "hlint: ignore Unused LANGUAGE pragma" #-}

-- | @'add' a b k@ computes @k x y@ such that
--
-- > x + y = a + b
-- > x = fl(a + b)
--
-- Which is to say that @x@ is the floating point image of @(a + b)@ and
-- @y@ stores the residual error term.
add :: Num a => a -> a -> (a -> a -> r) -> r
add a b k = k x y where
  x = a + b
  z = x - a
  y = (a - (x - z)) + (b - z)
{-# INLINE add #-}

-- | @'fadd' a b k@ computes @k x y@ such that
--
-- > x + y = a + b
-- > x = fl(a + b)
--
-- but only under the assumption that @'abs' a '>=' 'abs' b@. If you
-- aren't sure, use 'add'.
--
-- Which is to say that @x@ is the floating point image of @(a + b)@ and
-- @y@ stores the residual error term.
fadd :: Num a => a -> a -> (a -> a -> r) -> r
fadd a b k = k x (b - (x - a)) where
  x = a + b
{-# INLINE fadd #-}

-- | @'times' a b k@ computes @k x y@ such that
--
-- > x + y = a * b
-- > x = fl(a * b)
--
-- Which is to say that @x@ is the floating point image of @(a * b)@ and
-- @y@ stores the residual error term.
--
-- This could be nicer if we had access to a hardware fused multiply-add.
times :: Compensable a => a -> a -> (a -> a -> r) -> r
times a b k =
  split a $ \a1 a2 ->
  split b $ \b1 b2 ->
  let x = a * b in k x (a2*b2 - (((x - a1*b1) - a2*b1) - a1*b2))
  -- let x = a * b in k x (((a1*b1-x)+a1*b2+b2*b1)+a2*b2)
{-# INLINEABLE times #-}

-- this is a variant on a division algorithm by Liddicoat and Flynn
divide :: Compensable a => a -> a -> (a -> a -> r) -> r
divide a b = with (aX * ms) where
  x0   = recip b
  aX   = times a x0 compensated -- calculate aX
  m    = 1 +^ negate (times b x0 compensated)
  mm   = m*m
  ms   = 1+((m+mm)+m*mm)
{-# INLINEABLE divide #-}

-- | Priest's renormalization algorithm
--
-- @renorm a b c@ generates a 'Compensated' number assuming @a '>=' b '>=' c@.
renorm :: Compensable a => a -> a -> a -> Compensated a
renorm a b c =
  fadd b c $ \x1 y1 ->
  fadd a x1 $ \x2 y2 ->
  fadd x2 (y1 + y2) compensated
{-# INLINE renorm #-}

-- | @'squared' a k@ computes @k x y@ such that
--
-- > x + y = a * a
-- > x = fl(a * a)
--
-- Which is to say that @x@ is the floating point image of @(a * a)@ and
-- @y@ stores the residual error term.
squared :: Compensable a => a -> (a -> a -> r) -> r
squared a k =
  split a $ \a1 a2 ->
  let x = a * a in k x (a2*a2 - ((x - a1*a1) - 2*(a2*a1)))
{-# INLINE squared #-}

-- | Calculate a fast square of a compensated number.
square :: Compensable a => Compensated a -> Compensated a
square m =
  with m $ \a b ->
  squared a $ \x1 y1 ->
  times a b $ \x2 y2 ->
  add y1 (x2*2) $ \x3 y3 ->
  renorm x1 x3 (b*b + 2*y2 + y3)
{-# INLINE square #-}

-- | error-free split of a floating point number into two parts.
--
-- Note: these parts do not satisfy the `compensated` contract
split :: Compensable a => a -> (a -> a -> r) -> r
split a k = k x y where
  c = magic*a
  x = c - (c - a)
  y = a - x
{-# INLINEABLE split #-}

-- | Calculate a scalar + compensated sum with Kahan summation.
(+^) :: Compensable a => a -> Compensated a -> Compensated a
a +^ m = with m $ \b c -> let y = a - c; t = b + y in compensated t ((t - b) - y)
{-# INLINE (+^) #-}

-- | Compute @a * 'Compensated' a@
(*^) :: Compensable a => a -> Compensated a -> Compensated a
c *^ m =
  with m $ \ a b ->
  times c a $ \x1 y1 ->
  times c b $ \x2 y2 ->
  fadd x1 x2 $ \x3 y3 ->
  add y1 y3 $ \x4 y4 ->
  renorm x3 x4 (y4 + y2)
{-# INLINE (*^) #-}

class (RealFrac a, Floating a) => Compensable a where
  -- | This provides a numeric data type with effectively doubled precision by
  -- using Knuth's error free transform and a number of custom compensated
  -- arithmetic circuits.
  --
  -- This construction can be iterated, doubling precision each time.
  --
  -- >>> round (Prelude.product [2..100] :: Compensated (Compensated (Compensated Double)))
  -- 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
  --
  -- >>> Prelude.product [2..100]
  -- 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
  data Compensated a

  -- | This extracts both the 'primal' and 'residual' components of a 'Compensated'
  -- number.
  with :: Compensable a => Compensated a -> (a -> a -> r) -> r

  -- | Used internally to construct 'compensated' values that satisfy our residual contract.
  --
  -- When in doubt, use @'add' a b 'compensated'@ instead of @'compensated' a b@
  compensated :: Compensable a => a -> a -> Compensated a

  -- | This 'magic' number is used to 'split' the significand in half, so we can multiply
  -- them separately without losing precision in 'times'.
  magic :: a

instance Compensable Double where
  data Compensated Double = CD {-# UNPACK #-} !Double {-# UNPACK #-} !Double
  with (CD a b) k = k a b
  {-# INLINE with #-}
  compensated = CD
  {-# INLINE compensated #-}
  magic = 134217729
  {-# INLINE magic #-}

instance Compensable Float where
  data Compensated Float = CF {-# UNPACK #-} !Float {-# UNPACK #-} !Float
  with (CF a b) k = k a b
  {-# INLINE with #-}
  compensated = CF
  {-# INLINE compensated #-}
  magic = 4097
  {-# INLINE magic #-}

instance Compensable a => Compensable (Compensated a) where
  data Compensated (Compensated a) = CC !(Compensated a) !(Compensated a)
  with (CC a b) k = k a b
  {-# INLINE with #-}
  compensated = CC
  {-# INLINE compensated #-}
  magic = times (magic - 1) (magic - 1) $ \ x y -> compensated x (y + 1)
  {-# INLINE magic #-}

#if __GLASGOW_HASKELL__ < 707
instance Typeable1 Compensated where
  typeOf1 _ = mkTyConApp (mkTyCon3 "analytics" "Data.Analytics.Numeric.Compensated" "Compensated") []
#else
deriving instance Typeable Compensated
#endif

instance (Compensable a, Hashable a) => Hashable (Compensated a) where
  hashWithSalt n m = with m $ \a b -> hashWithSalt n (a,b)

instance (Compensable a, Data a) => Data (Compensated a) where
  gfoldl f z m = with m $ \a b -> z compensated `f` a `f` b
  toConstr _ = compensatedConstr
  gunfold k z c = case constrIndex c of
    1 -> k (k (z compensated))
    _ -> error "gunfold"
  dataTypeOf _ = compensatedDataType
  dataCast1 f = gcast1 f

compensatedConstr :: Constr
compensatedConstr = mkConstr compensatedDataType "compensated" [] Prefix
{-# NOINLINE compensatedConstr #-}

compensatedDataType :: DataType
compensatedDataType = mkDataType "Data.Analytics.Numeric.Compensated" [compensatedConstr]
{-# NOINLINE compensatedDataType #-}

instance (Compensable a, NFData a) => NFData (Compensated a) where
  rnf m = with m $ \x y -> rnf x `seq` rnf y
  {-# INLINE rnf #-}

instance (Compensable a, Show a) => Show (Compensated a) where
  showsPrec d m = with m $ \a b -> showParen (d > 10) $
    showString "compensated " . T.showsPrec 11 a . showChar ' ' . T.showsPrec 11 b

instance (Compensable a, Read a) => Read (Compensated a) where
  readPrec = parens $ prec 10 $ do
    Ident "compensated" <- lexP
    a <- step T.readPrec
    b <- step T.readPrec
    return $ compensated a b

type Overcompensated a = Compensated (Compensated a)

-- | This provides the isomorphism between the compact representation we store these in internally
-- and the naive pair of the 'primal' and 'residual' components.
_Compensated :: Compensable a => Iso' (Compensated a) (a, a)
_Compensated = iso (`with` (,)) (uncurry compensated)
{-# INLINE _Compensated #-}

-- | This 'Lens' lets us edit the 'primal' directly, leaving the 'residual' untouched.
primal :: Compensable a => Lens' (Compensated a) a
primal f c = with c $ \a b -> f a <&> \a' -> compensated a' b
{-# INLINE primal #-}

-- | This 'Lens' lets us edit the 'residual' directly, leaving the 'primal' untouched.
residual :: Compensable a => Lens' (Compensated a) a
residual f c = with c $ \a b -> compensated a <$> f b
{-# INLINE residual #-}

-- | Extract the 'primal' component of a 'compensated' value, when and if compensation
-- is no longer required.
uncompensated :: Compensable a => Compensated a -> a
uncompensated c = with c const
{-# INLINE uncompensated #-}

type instance Index (Compensated a) = Int
instance (Compensable a, Compensable b) => Each (Compensated a) (Compensated b) a b where
  each f m = with m $ \a b -> compensated <$> f a <*> f b
  {-# INLINE each #-}

instance Compensable a => Eq (Compensated a) where
  m == n = with m $ \a b -> with n $ \c d -> a == c && b == d
  m /= n = with m $ \a b -> with n $ \c d -> a /= c || b /= d
  {-# INLINE (==) #-}

instance Compensable a => Ord (Compensated a) where
  compare m n = with m $ \a b -> with n $ \c d -> compare a c <> compare b d
  {-# INLINE compare #-}
  m <= n = with m $ \a b -> with n $ \c d -> case compare a c of
    LT -> True
    EQ -> b <= d
    GT -> False
  {-# INLINE (<=) #-}
  m >= n = with m $ \a b -> with n $ \c d -> case compare a c of
    LT -> False
    EQ -> b >= d
    GT -> a >= c -- @compare x NaN@ and @compare naN x@ always return 'GT', but @m >= n@ should be 'False'
  {-# INLINE (>=) #-}
  m > n = with m $ \a b -> with n $ \c d -> case compare a c of
    LT -> False
    EQ -> b > d
    GT -> a > c -- @compare x NaN@ and @compare naN x@ always return 'GT', but @m >= n@ should be 'False'
  {-# INLINE (>) #-}
  m < n = with m $ \a b -> with n $ \c d -> case compare a c of
    LT -> True
    EQ -> b < d
    GT -> False
  {-# INLINE (<) #-}

instance Compensable a => Semigroup (Compensated a) where
  (<>) = (+)
  {-# INLINE (<>) #-}

instance Compensable a => Monoid (Compensated a) where
  mempty = compensated 0 0
  {-# INLINE mempty #-}
  mappend = (+)
  {-# INLINE mappend #-}

-- | Perform Kahan summation over a list.
kahan :: (Foldable f, Compensable a) => f a -> Compensated a
kahan = Foldable.foldr (+^) mempty
{-# INLINE kahan #-}

instance Compensable a => Num (Compensated a) where
  m + n =
    with m $ \a  b  ->
    with n $ \c d ->
    add  a  c $ \x1 y1 ->
    add y1  d $ \x2 y2 ->
    add  b x2 $ \x3 y3 ->
    add x1 x3 $ \x4 y4 ->
    add x4 (y2 + y3 + y4) compensated
  {-# INLINE (+) #-}

{-
  m + n =
    with m $ \a b ->
    with n $ \c d ->
    add a c $ \x1 y1 ->
    add b d $ \x2 y2 ->
    renorm x1 x2 (y1 + y2)
  {-# INLINE (+) #-}
-}


  m * n =
    with m $ \a b ->
    with n $ \c d ->
    times a c $ \x1 y1 ->
    times b c $ \x2 y2 ->
    times a d $ \x3 y3 ->
    add x1 x2 $ \x4 y4 ->
    add x3 x4 $ \x5 y5 ->
    add y1 y4 $ \x6 y6 ->
    add y5 x6 $ \x7 y7 ->
    add x5 x7 $ \x8 y8 ->
    add x8 (b*d + y2 + y3 + y6 + y7 + y8) compensated
  {-# INLINE (*) #-}

{-
  m * n =
    with m $ \a b ->
    with n $ \c d ->
    times a c $ \x1 y1 ->
    times b c $ \x2 y2 ->
    times a d $ \x3 y3 ->
    add y1 x2 $ \x4 y4 ->
    add x3 x4 $ \x5 y5 ->
    renorm x1 x5 (b * d + y2 + y4 + y3 + y5)
  {-# INLINE (*) #-}
-}

  negate m = with m (on compensated negate)
  -- {-# INLINE negate #-}

  x - y = x + negate y
  {-# INLINE (-) #-}

  signum m = with m $ \a _ -> compensated (signum a) 0
  {-# INLINE signum #-}

  abs m = with m $ \a b ->
    if a < 0
    then compensated (negate a) (negate b)
    else compensated a b
  {-# INLINE abs #-}

  fromInteger i = add x (fromInteger (i - round x)) compensated where
    x = fromInteger i
  {-# INLINE fromInteger #-}

instance Compensable a => Enum (Compensated a) where
  succ a = a + 1
  {-# INLINE succ #-}
  pred a = a - 1
  {-# INLINE pred #-}
  toEnum i = add x (fromIntegral (i - round x)) compensated where
    x = fromIntegral i
  {-# INLINE toEnum #-}
  fromEnum = round
  {-# INLINE fromEnum #-}
  enumFrom a = a : Prelude.enumFrom (a + 1)
  {-# INLINE enumFrom #-}
  enumFromThen a b = a : Prelude.enumFromThen b (b - a + b)
  {-# INLINE enumFromThen #-}
  enumFromTo a b
    | a <= b = a : Prelude.enumFromTo (a + 1) b
    | otherwise = []
  {-# INLINE enumFromTo #-}
  enumFromThenTo a b c
    | a <= b    = up a
    | otherwise = down a
    where
     delta = b - a
     up x | x <= c    = x : up (x + delta)
          | otherwise = []
     down x | c <= x    = x : down (x + delta)
            | otherwise = []
  {-# INLINE enumFromThenTo #-}

instance Compensable a => Fractional (Compensated a) where
  recip m = with m $ \a b -> add (recip a) (-b / (a * a)) compensated
  {-# INLINE recip #-}

  -- | A variant on a hardware division algorithm by Liddicoat and Flynn
  a / b = (a*x0) * (1+((m+mm)+m*mm)) where
    x0  = recip b
    m   = 1 - b*x0
    mm  = m*m
  -- {-# INLINE (/) #-}

  fromRational r = fromInteger (numerator r) / fromInteger (denominator r)
  -- {-# INLINE fromRational #-}

instance Compensable a => Real (Compensated a) where
  toRational m = with m (on (+) toRational)
  -- {-# INLINE toRational #-}

instance Compensable a => RealFrac (Compensated a) where
  properFraction m = with m $ \a b -> case properFraction a of
    (w, p) -> add p b $ \ x y -> case properFraction x of
      (w',q) -> (w + w', add q y compensated)
  -- {-# INLINE properFraction #-}

instance (Compensable a, Binary a) => Binary (Compensated a) where
  get = compensated <$> Binary.get <*> Binary.get
  put m = with m $ \a b -> do
    Binary.put a
    Binary.put b

instance (Compensable a, Serialize a) => Serialize (Compensated a) where
  get = compensated <$> Serialize.get <*> Serialize.get
  put m = with m $ \a b -> do
    Serialize.put a
    Serialize.put b

instance (Compensable a, Serial a) => Serial (Compensated a) where
  deserialize = compensated <$> Bytes.deserialize <*> Bytes.deserialize
  serialize m = with m $ \a b -> do
    Bytes.serialize a
    Bytes.serialize b

-- ಠ_ಠ this unnecessarily expects that the format won't change, because I can't derive a better instance.
instance (Compensable a, Serialize a) => SafeCopy (Compensated a) where
  -- safecopy-0.10.0 changed its default implementations for these methods.
  -- The implementations below are copied from the pre-0.10.0 defaults.
  errorTypeName _ = "<unknown type>"
  getCopy = contain Serialize.get
  putCopy = contain . Serialize.put

instance (Compensable a, Storable a) => Storable (Compensated a) where
  sizeOf _ = sizeOf (undefined :: a) * 2
  -- {-# INLINE sizeOf #-}
  alignment _ = alignment (undefined :: a)
  -- {-# INLINE alignment #-}
  peekElemOff p o | q <- castPtr p, o2 <- o * 2 =
    compensated <$> peekElemOff q o2 <*> peekElemOff q (o2+1)
  -- {-# INLINE peekElemOff #-}
  pokeElemOff p o m | q <- castPtr p, o2 <- o * 2 = with m $ \a b -> do
    pokeElemOff q o2 a
    pokeElemOff q (o2+1) b
  -- {-# INLINE pokeElemOff #-}
  peekByteOff p o | q <- castPtr p =
    compensated <$> peekByteOff q o <*> peekByteOff q (o + sizeOf (undefined :: a))
  -- {-# INLINE peekByteOff #-}
  pokeByteOff p o m | q <- castPtr p = with m $ \a b -> do
    pokeByteOff q o a
    pokeByteOff q (o+sizeOf (undefined :: a)) b
  -- {-# INLINE pokeByteOff #-}
  peek p | q <- castPtr p = compensated <$> peek q <*> peekElemOff q 1
  -- {-# INLINE peek #-}
  poke p m | q <- castPtr p = with m $ \a b -> do
    poke q a
    pokeElemOff q 1 b
  -- {-# INLINE poke #-}

newtype instance U.MVector s (Compensated a) = MV_Compensated (U.MVector s (a,a))
newtype instance U.Vector (Compensated a) = V_Compensated (U.Vector (a, a))

instance (Compensable a, Unbox a) => M.MVector U.MVector (Compensated a) where
  basicLength (MV_Compensated v) = M.basicLength v
  {-# INLINE basicLength #-}
  basicUnsafeSlice i n (MV_Compensated v) = MV_Compensated $ M.basicUnsafeSlice i n v
  {-# INLINE basicUnsafeSlice #-}
  basicOverlaps (MV_Compensated v1) (MV_Compensated v2) = M.basicOverlaps v1 v2
  {-# INLINE basicOverlaps #-}
  basicUnsafeNew n = MV_Compensated <$> M.basicUnsafeNew n
  {-# INLINE basicUnsafeNew #-}
  basicUnsafeReplicate n m = with m $ \x y -> MV_Compensated <$> M.basicUnsafeReplicate n (x,y)
  {-# INLINE basicUnsafeReplicate #-}
  basicUnsafeRead (MV_Compensated v) i = uncurry compensated <$> M.basicUnsafeRead v i
  {-# INLINE basicUnsafeRead #-}
  basicUnsafeWrite (MV_Compensated v) i m = with m $ \ x y -> M.basicUnsafeWrite v i (x,y)
  {-# INLINE basicUnsafeWrite #-}
  basicClear (MV_Compensated v) = M.basicClear v
  {-# INLINE basicClear #-}
  basicSet (MV_Compensated v) m = with m $ \ x y -> M.basicSet v (x,y)
  {-# INLINE basicSet #-}
  basicUnsafeCopy (MV_Compensated v1) (MV_Compensated v2) = M.basicUnsafeCopy v1 v2
  {-# INLINE basicUnsafeCopy #-}
  basicUnsafeMove (MV_Compensated v1) (MV_Compensated v2) = M.basicUnsafeMove v1 v2
  {-# INLINE basicUnsafeMove #-}
  basicUnsafeGrow (MV_Compensated v) n = MV_Compensated <$> M.basicUnsafeGrow v n
  {-# INLINE basicUnsafeGrow #-}
#if MIN_VERSION_vector(0,11,0)
  basicInitialize (MV_Compensated v) = M.basicInitialize v
  {-# INLINE basicInitialize #-}
#endif

instance (Compensable a, Unbox a) => G.Vector U.Vector (Compensated a) where
  basicUnsafeFreeze (MV_Compensated v) = V_Compensated <$> G.basicUnsafeFreeze v
  {-# INLINE basicUnsafeFreeze #-}
  basicUnsafeThaw (V_Compensated v) = MV_Compensated <$> G.basicUnsafeThaw v
  {-# INLINE basicUnsafeThaw #-}
  basicLength (V_Compensated v) = G.basicLength v
  {-# INLINE basicLength #-}
  basicUnsafeSlice i n (V_Compensated v) = V_Compensated $ G.basicUnsafeSlice i n v
  {-# INLINE basicUnsafeSlice #-}
  basicUnsafeIndexM (V_Compensated v) i
                = uncurry compensated <$> G.basicUnsafeIndexM v i
  {-# INLINE basicUnsafeIndexM #-}
  basicUnsafeCopy (MV_Compensated mv) (V_Compensated v)
                = G.basicUnsafeCopy mv v
  {-# INLINE basicUnsafeCopy #-}
  elemseq _ m z = with m $ \x y -> G.elemseq (undefined :: U.Vector a) x
                                 $ G.elemseq (undefined :: U.Vector a) y z
  {-# INLINE elemseq #-}

-- | /NB:/ Experimental and partially implemented.
--
-- Other than sqrt, the accuracy of these is basically uncalculated! In fact many of these are known to be wrong! Patches and improvements are welcome.
instance Compensable a => Floating (Compensated a) where
#ifdef SPECIALIZE_INSTANCES
  {-# SPECIALIZE instance Floating (Compensated Double) #-}
  {-# SPECIALIZE instance Floating (Compensated Float) #-}
  {-# SPECIALIZE instance Compensable a => Floating (Compensated (Compensated a)) #-}
#endif

  exp m =
    with m $ \a b ->
    times (exp a) (exp b) compensated

  sin m =
    with m $ \a b ->
    times (sin a) (cos b) $ \x1 y1 ->
    times (sin b) (cos a) $ \x2 y2 ->
    add x1 x2 $ \x3 y3 ->
    add y1 y2 $ \x4 y4 ->
    add x4 y3 $ \x5 y5 ->
    add x5 x3 $ \x6 y6 ->
    add (y4 + y5 + y6) x6 compensated

  cos m =
    with m $ \a b ->
    times (cos a) (cos b) $ \x1 y1 ->
    times (-sin b) (sin a) $ \x2 y2 ->
    add x1 x2 $ \x3 y3 ->
    add y1 y2 $ \x4 y4 ->
    add x4 y3 $ \x5 y5 ->
    add x5 x3 $ \x6 y6 ->
    add (y4 + y5 + y6) x6 compensated

  tan m =
    with m $ \a b ->
    add (tan a) (tan b) compensated /
    (1 +^ times (tan a) (tan b) compensated)

  sinh m =
    with m $ \a b ->
    times (sinh a) (cosh b) $ \x1 y1 ->
    times (cosh a) (sinh b) $ \x2 y2 ->
    add x1 x2 $ \x3 y3 ->
    add y1 y2 $ \x4 y4 ->
    add x4 y3 $ \x5 y5 ->
    add x5 x3 $ \x6 y6 ->
    add (y4 + y5 + y6) x6 compensated

  cosh m =
    with m $ \a b ->
    times (cosh a) (cosh b) $ \x1 y1 ->
    times (sinh b) (sinh a) $ \x2 y2 ->
    add x1 x2 $ \x3 y3 ->
    add y1 y2 $ \x4 y4 ->
    add x4 y3 $ \x5 y5 ->
    add x5 x3 $ \x6 y6 ->
    add (y4 + y5 + y6) x6 compensated

  tanh m =
    with m $ \a b ->
    fadd (tanh a) (tanh b) compensated /
    (1 +^ times (tanh a) (tanh b) compensated)

  -- This requires an accurate 'exp', which we currently lack.
  log m =
    with m $ \ a b -> let
      xy1 = add (log a) (b/a) compensated
      xy2 = xy1 + m * exp (-xy1) - 1 -- Newton Raphson step 1
    in xy2 + m * exp (-xy2) - 1      -- Newton Raphson step 2

  -- | Hardware sqrt improved by the Babylonian algorithm (Newton Raphson)
  sqrt m = with (z4 + m/z4) $ on compensated (/2) where
    z0 = sqrt (m^.primal)
    z1 = with (z0 +^ (m / compensated z0 0)) $ on compensated (/2)
    z2 = with (z1 + m/z1) $ on compensated (/2)
    z3 = with (z2 + m/z2) $ on compensated (/2)
    z4 = with (z3 + m/z3) $ on compensated (/2)

  -- TODO: do log1p, expm1, log1mexp, log1pexp right!
  log1p a = log (1 + a)
  {-# INLINE log1p #-}
  expm1 a = exp a - 1
  {-# INLINE expm1 #-}
  log1mexp a | a <= log 2 = log (negate (expm1 (negate a)))
             | otherwise  = log1p (negate (exp (negate a)))
  {-# INLINE log1mexp #-}
  log1pexp a
    | a <= 18   = log1p (exp a)
    | a <= 100  = a + exp (negate a)
    | otherwise = a
  {-# INLINE log1pexp #-}

  -- (**)    = error "TODO"
  pi      = error "TODO"
  asin    = error "TODO"
  atan    = error "TODO"
  acos    = error "TODO"
  asinh   = error "TODO"
  atanh   = error "TODO"
  acosh   = error "TODO"