content: {-# LANGUAGE CPP #-}
  {-# LANGUAGE ForeignFunctionInterface #-}
  --------------------------------------------------------------------
  -- |
  -- Copyright :  (c) Edward Kmett 2013
  -- License   :  BSD3
  -- Maintainer:  Edward Kmett <ekmett@gmail.com>
  -- Stability :  experimental
  -- Portability: non-portable
  --
  -- These functions provide wildly inaccurate but very fast
  -- approximations to common transcendental functions.
  --
  -- The algorithms here are based on Martin Ankerl's optimized 'pow',
  -- <http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/>
  -- which is in turn based on
  -- <http://nic.schraudolph.org/pubs/Schraudolph99.pdf>
  --------------------------------------------------------------------
  module Data.Approximate.Numerics
    ( Fast(..)
    , blog
    ) where
  
  class Floating a => Fast a where
    -- | Calculate an approximate log.
    flog    :: a -> a
    flog_lb :: a -> a
    flog_ub :: a -> a
    -- | Calculate an approximate exp.
    fexp    :: a -> a
    fexp_lb :: a -> a
    fexp_ub :: a -> a
    -- | Calculate an approximate pow.
    fpow    :: a -> a -> a
    fpow_lb :: a -> a -> a
    fpow_ub :: a -> a -> a
  
  instance Fast Float where
    flog     = logf_fast
    flog_lb  = logf_fast_lb
    flog_ub  = logf_fast_ub
    fexp     = better_expf_fast
    fexp_lb  = expf_fast_lb
    fexp_ub  = expf_fast_ub
    fpow     = better_powf_fast_precise
    fpow_lb  = powf_fast_lb
    fpow_ub  = powf_fast_ub
  
  instance Fast Double where
    flog     = log_fast
    flog_lb  = log_fast_lb
    flog_ub  = log_fast_ub
    fexp     = better_exp_fast
    fexp_lb  = exp_fast_lb
    fexp_ub  = exp_fast_ub
    fpow     = better_pow_fast_precise
    fpow_lb  = pow_fast_lb
    fpow_ub  = pow_fast_ub
  
  -- | Borchardt’s Algorithm from “Dead Reckoning: Calculating without instruments”.
  --
  -- This is a remarkably bad approximate logarithm.
  --
  -- 'flog' had better outperform it! It is provided merely for comparison.
  blog :: Floating a => a -> a
  blog x = 6 * (x - 1) / (x + 1 + 4 * sqrt x)
        -- 6 / (x + 1 + 4 * sqrt x) / (x-1) -- herbie's suggested replacement.
  #ifdef HERBIE
  {-# ANN blog "NoHerbie" #-}
  #endif
  
  foreign import ccall unsafe pow_fast_lb       :: Double -> Double -> Double
  foreign import ccall unsafe pow_fast_ub       :: Double -> Double -> Double
  foreign import ccall unsafe better_pow_fast_precise  :: Double -> Double -> Double
  foreign import ccall unsafe powf_fast_lb      :: Float -> Float -> Float
  foreign import ccall unsafe powf_fast_ub      :: Float -> Float -> Float
  foreign import ccall unsafe better_powf_fast_precise :: Float -> Float -> Float
  
  foreign import ccall unsafe better_exp_fast   :: Double -> Double
  foreign import ccall unsafe exp_fast_lb       :: Double -> Double
  foreign import ccall unsafe exp_fast_ub       :: Double -> Double
  foreign import ccall unsafe better_expf_fast  :: Float -> Float
  foreign import ccall unsafe expf_fast_lb      :: Float -> Float
  foreign import ccall unsafe expf_fast_ub      :: Float -> Float
  
  foreign import ccall unsafe log_fast        :: Double -> Double
  foreign import ccall unsafe log_fast_lb     :: Double -> Double
  foreign import ccall unsafe log_fast_ub     :: Double -> Double
  foreign import ccall unsafe logf_fast       :: Float -> Float
  foreign import ccall unsafe logf_fast_lb    :: Float -> Float
  foreign import ccall unsafe logf_fast_ub    :: Float -> Float
  