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