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.Maybe
  19
  20 import Foreign
  21 import Foreign.C
  22
  23 foreign import ccall unsafe "double_format"
  24     double_format :: CString -> CChar -> CInt -> CDouble -> IO CInt
  25
  26 -- | @binarySearch p low high@ computes the least integer on the interval
  27 -- [low, high] satisfying the given predicate, by binary search.  The
  28 -- predicate must partition the interval into two contiguous regions.
  29 binarySearch :: Integral a => (a -> Bool) -> a -> a -> a
  30 binarySearch p l u
  31     | l >  u    = error "empty interval"
  32     | l == u    = m
  33     | p m       = binarySearch p l m
  34     | otherwise = binarySearch p (m+1) u
  35     where
  36         m = l + div (u-l) 2
  37
  38 -- | Find a power of two such that the given rational number, when multiplied
  39 -- by that power and rounded to an integer, has exactly as many digits as the
  40 -- precision of the floating point type.  The search may stop for values with
  41 -- extremely large (or small) magnitude, in which case the result shall
  42 -- overflow (or underflow) when scaled to the floating type.
  43 scaleRational :: PrimFloat a => a -> Rational -> (Integer, Int)
  44 scaleRational t x = (fromJust . toIntegral . round . scale x $ e, e) where
  45     e         = binarySearch ((lbound <=) . scale x) l (u*2)
  46     (l, u)    = floatRange t
  47     lbound    = floatRadix t ^ (floatPrecision t - 1)
  48     scale x y = x * 2^^y
  49
  50 formatDouble :: Char -> Int -> Double -> String
  51 formatDouble c p x = unsafePerformIO $ do
  52     let format = castCharToCChar c
  53     size <- double_format nullPtr format (fromIntegral p) (toFloating x)
  54     allocaArray0 (fromIntegral size) $ \buf -> do
  55         double_format buf format (fromIntegral p) (toFloating x)
  56         peekCString buf