Data/Floating/Helpers.hs

   1 {-
   2  - Copyright (C) 2010 Nick Bowler.
   3  -
   4  - License BSD2:  2-clause BSD license.  See LICENSE for full terms.
   5  - This is free software: you are free to change and redistribute it.
   6  - There is NO WARRANTY, to the extent permitted by law.
   7  -}
   8
   9 -- | Internal helper functions needed by at least two modules.
  10 {-# LANGUAGE ForeignFunctionInterface #-}
  11 module Data.Floating.Helpers (
  12     binarySearch, scaleRational, formatDouble
  13 ) where
  14
  15 import Prelude hiding (Double, RealFloat(..), RealFrac(..))
  16 import Data.Floating.Types.Core
  17 import Data.Roundable
  18 import Data.Ratio
  19 import Data.Maybe
  20
  21 import Foreign
  22 import Foreign.C
  23 import System.IO.Unsafe
  24
  25 foreign import ccall unsafe "double_format"
  26     double_format :: CString -> CChar -> CInt -> CDouble -> IO CInt
  27
  28 -- | @binarySearch p low high@ computes the least integer on the interval
  29 -- [low, high] satisfying the given predicate, by binary search.  The
  30 -- predicate must partition the interval into two contiguous regions.
  31 binarySearch :: Integral a => (a -> Bool) -> a -> a -> a
  32 binarySearch p l u
  33     | l >  u    = error "empty interval"
  34     | l == u    = m
  35     | p m       = binarySearch p l m
  36     | otherwise = binarySearch p (m+1) u
  37     where
  38         m = l + div (u-l) 2
  39
  40 -- | Find a power of two such that the given rational number, when multiplied
  41 -- by that power and rounded to an integer, has exactly as many digits as the
  42 -- precision of the floating point type.  The search may stop for values with
  43 -- extremely large (or small) magnitude, in which case the result shall
  44 -- overflow (or underflow) when scaled to the floating type.
  45 scaleRational :: PrimFloat a => a -> Rational -> (Integer, Int)
  46 scaleRational t x = (fromJust . toIntegral . round . scale x $ e, e) where
  47     e         = binarySearch ((lbound <=) . scale x) l (u*2)
  48     (l, u)    = floatRange t
  49     lbound    = floatRadix t ^ (floatPrecision t - 1)
  50     scale x y = x * 2^^y
  51
  52 formatDouble :: Char -> Int -> Double -> String
  53 formatDouble c p x = unsafePerformIO $ do
  54     let format = castCharToCChar c
  55     size <- double_format nullPtr format (fromIntegral p) (toFloating x)
  56     allocaArray0 (fromIntegral size) $ \buf -> do
  57         double_format buf format (fromIntegral p) (toFloating x)
  58         peekCString buf