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 module Data.Floating.Helpers where
  11
  12 import Prelude hiding (RealFloat(..), RealFrac(..))
  13 import Data.Floating.Classes
  14 import Data.Floating.Instances
  15 import Data.Ratio
  16 import Data.Maybe
  17
  18 -- | @binarySearch p low high@ computes the least integer on the interval
  19 -- [low, high] satisfying the given predicate, by binary search.  The
  20 -- predicate must partition the interval into two contiguous regions.
  21 binarySearch :: Integral a => (a -> Bool) -> a -> a -> a
  22 binarySearch p l u
  23     | l >  u    = error "empty interval"
  24     | l == u    = m
  25     | p m       = binarySearch p l m
  26     | otherwise = binarySearch p (m+1) u
  27     where
  28         m = l + div (u-l) 2
  29
  30 -- | Find a power of two such that the given rational number, when multiplied
  31 -- by that power and rounded to an integer, has exactly as many digits as the
  32 -- precision of the floating point type.  The search may stop for values with
  33 -- extremely large (or small) magnitude, in which case the result shall
  34 -- overflow (or underflow) when scaled to the floating type.
  35 scaleRational :: PrimFloat a => a -> Rational -> (Integer, Int)
  36 scaleRational t x = (fromJust . toIntegral . round . scale x $ e, e) where
  37     e         = binarySearch ((lbound <) . scale x) l (u*2)
  38     (l, u)    = floatRange t
  39     lbound    = floatRadix t ^ (floatPrecision t - 1)
  40     scale x y = x * 2^^y