content: % Copyright (C) 2009 John Millikin <jmillikin@gmail.com>
  % 
  % This program is free software: you can redistribute it and/or modify
  % it under the terms of the GNU General Public License as published by
  % the Free Software Foundation, either version 3 of the License, or
  % any later version.
  % 
  % This program is distributed in the hope that it will be useful,
  % but WITHOUT ANY WARRANTY; without even the implied warranty of
  % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  % GNU General Public License for more details.
  % 
  % You should have received a copy of the GNU General Public License
  % along with this program.  If not, see <http://www.gnu.org/licenses/>.
  
  \ignore{
  \begin{code}
  {-# LANGUAGE GeneralizedNewtypeDeriving #-}
  module Data.Binary.IEEE754 (
      -- * Parsing
        getFloat16be, getFloat16le
      , getFloat32be, getFloat32le
      , getFloat64be, getFloat64le
  
      -- * Serializing
      , putFloat32be, putFloat32le
      , putFloat64be, putFloat64le
  ) where
  
  import Data.Bits ((.&.), (.|.), shiftL, shiftR, Bits)
  import Data.Word (Word8)
  import Data.List (foldl')
  
  import qualified Data.ByteString as B
  import Data.Binary.Get (Get, getByteString)
  import Data.Binary.Put (Put, putByteString)
  \end{code}
  }
  
  \section{Parsing}
  
  \subsection{Public interface}
  
  \begin{code}
  getFloat16be :: Get Float
  getFloat16be = getFloat (ByteCount 2) splitBytes
  \end{code}
  
  \begin{code}
  getFloat16le :: Get Float
  getFloat16le = getFloat (ByteCount 2) $ splitBytes . reverse
  \end{code}
  
  \begin{code}
  getFloat32be :: Get Float
  getFloat32be = getFloat (ByteCount 4) splitBytes
  \end{code}
  
  \begin{code}
  getFloat32le :: Get Float
  getFloat32le = getFloat (ByteCount 4) $ splitBytes . reverse
  \end{code}
  
  \begin{code}
  getFloat64be :: Get Double
  getFloat64be = getFloat (ByteCount 8) splitBytes
  \end{code}
  
  \begin{code}
  getFloat64le :: Get Double
  getFloat64le = getFloat (ByteCount 8) $ splitBytes . reverse
  \end{code}
  
  \subsection{Implementation}
  
  Split the raw byte array into (sign, exponent, significand) components.
  The exponent and signifcand are drawn directly from the bits in the
  original float, and have not been unbiased or otherwise modified.
  
  \begin{code}
  splitBytes :: [Word8] -> RawFloat
  splitBytes bs = RawFloat width sign exp' sig expWidth sigWidth where
      width = ByteCount (length bs)
      nBits = bitsInWord8 bs
      sign = if head bs .&. 0x80 == 0x80
               then Negative
               else Positive
  
      expStart = 1
      expWidth = exponentWidth nBits
      expEnd = expStart + expWidth
      exp' = Exponent . fromIntegral $ bitSlice bs expStart expEnd
  
      sigWidth = nBits - expEnd
      sig  = Significand $ bitSlice bs expEnd nBits
  \end{code}
  
  \subsubsection{Encodings and special values}
  
  The next step depends on the value of the exponent $e$, size of the
  exponent field in bits $w$, and value of the significand.
  
  \begin{table}[h]
  \begin{center}
  \begin{tabular}{lrl}
  \toprule
  Exponent                & Significand & Format \\
  \midrule
  $0$                     & $0$           & Zero \\
  $0$                     & $> 0$         & Denormalised \\
  $1 \leq e \leq 2^w - 2$ & \textit{any} & Normalised \\
  $2^w-1$                 & $0$           & Infinity \\
  $2^w-1$                 & $> 0$         & NaN \\
  \bottomrule
  \end{tabular}
  \end{center}
  \end{table}
  
  There's no built-in literals for Infinity or NaN, so they
  are constructed using the {\tt Read} instances for {\tt Double} and
  {\tt Float}.
  
  \begin{code}
  merge :: (Read a, RealFloat a) => RawFloat -> a
  merge f@(RawFloat _ _ e sig eWidth _)
      | e == 0 = if sig == 0
                   then 0.0
                   else denormalised f
      | e == eMax - 1 = if sig == 0
                          then read "Infinity"
                          else read "NaN"
      | otherwise = normalised f
      where eMax = 2 `pow` eWidth
  \end{code}
  
  If a value is normalised, its significand has an implied {\tt 1} bit
  in its most-significant bit. The significand must be adjusted by
  this value before being passed to {\tt encodeField}.
  
  \begin{code}
  normalised :: RealFloat a => RawFloat -> a
  normalised f = encodeFloat fraction exp' where
      Significand sig = rawSignificand f
      Exponent exp' = unbiased - sigWidth
  
      fraction = sig + (1 `bitShiftL` rawSignificandWidth f)
  
      sigWidth = fromIntegral $ rawSignificandWidth f
      unbiased = unbias (rawExponent f) (rawExponentWidth f)
  \end{code}
  
  For denormalised values, the implied {\tt 1} bit is the least-significant
  bit of the exponent.
  
  \begin{code}
  denormalised :: RealFloat a => RawFloat -> a
  denormalised f = encodeFloat sig exp' where
      Significand sig = rawSignificand f
      Exponent exp' = unbiased - sigWidth + 1
  
      sigWidth = fromIntegral $ rawSignificandWidth f
      unbiased = unbias (rawExponent f) (rawExponentWidth f)
  \end{code}
  
  By composing {\tt splitBytes} and {\tt merge}, the absolute value of the
  float is calculated. Before being returned to the calling function, it
  must be signed appropriately.
  
  \begin{code}
  getFloat :: (Read a, RealFloat a) => ByteCount
              -> ([Word8] -> RawFloat) -> Get a
  getFloat (ByteCount width) parser = do
      raw <- fmap (parser . B.unpack) $ getByteString width
      let absFloat = merge raw
      return $ case rawSign raw of
                 Positive ->  absFloat
                 Negative -> -absFloat
  \end{code}
  
  \section{Serialising}
  
  \subsection{Public interface}
  
  \begin{code}
  putFloat32be :: Float -> Put
  putFloat32be = putFloat (ByteCount 4) id
  \end{code}
  
  \begin{code}
  putFloat32le :: Float -> Put
  putFloat32le = putFloat (ByteCount 4) reverse
  \end{code}
  
  \begin{code}
  putFloat64be :: Double -> Put
  putFloat64be = putFloat (ByteCount 8) id
  \end{code}
  
  \begin{code}
  putFloat64le :: Double -> Put
  putFloat64le = putFloat (ByteCount 8) reverse
  \end{code}
  
  \subsection{Implementation}
  
  Serialisation is similar to parsing. First, the float is converted to
  a structure representing raw bitfields. The values returned from
  {\tt decodeFloat} are clamped to their expected lengths before being
  stored in the {\tt RawFloat}.
  
  \begin{code}
  splitFloat :: RealFloat a => ByteCount -> a -> RawFloat
  splitFloat width x = raw where
      raw = RawFloat width sign clampedExp clampedSig expWidth sigWidth
      sign = if isNegativeNaN x || isNegativeZero x || x < 0
               then Negative
               else Positive
      clampedExp = clamp expWidth exp'
      clampedSig = clamp sigWidth sig
      (exp', sig) = case (dFraction, dExponent, biasedExp) of
                      (0, 0, _) -> (0, 0)
                      (_, _, 0) -> (0, Significand $ truncatedSig + 1)
                      _         -> (biasedExp, Significand truncatedSig)
      expWidth = exponentWidth $ bitCount width
      sigWidth = bitCount width - expWidth - 1 -- 1 for sign bit
  
      (dFraction, dExponent) = decodeFloat x
  
      rawExp = Exponent $ dExponent + fromIntegral sigWidth
      biasedExp = bias rawExp expWidth
      truncatedSig = abs dFraction - (1 `bitShiftL` sigWidth)
  \end{code}
  
  Then, the {\tt RawFloat} is converted to a list of bytes by mashing all
  the fields together into an {\tt Integer}, and chopping up that integer
  in 8-bit blocks.
  
  \begin{code}
  rawToBytes :: RawFloat -> [Word8]
  rawToBytes raw = integerToBytes mashed width where
      RawFloat width sign exp' sig expWidth sigWidth = raw
      sign' :: Word8
      sign' = case sign of
                Positive -> 0
                Negative -> 1
      mashed = mashBits sig sigWidth .
               mashBits exp' expWidth .
               mashBits sign' 1 $ 0
  \end{code}
  
  {\tt clamp}, given a maximum bit count and a value, will strip any 1-bits
  in positions above the count.
  
  \begin{code}
  clamp :: (Num a, Bits a) => BitCount -> a -> a
  clamp = (.&.) . mask where
      mask 1 = 1
      mask n | n > 1 = (mask (n - 1) `shiftL` 1) + 1
      mask _ = undefined
  \end{code}
  
  For merging the fields, just shift the starting integer over a bit and
  then \textsc{or} it with the next value. The weird parameter order allows
  easy composition.
  
  \begin{code}
  mashBits :: (Bits a, Integral a) => a -> BitCount -> Integer -> Integer
  mashBits _ 0 x = x
  mashBits y n x = (x `bitShiftL` n) .|. fromIntegral y
  \end{code}
  
  Given an integer, read it in 255-block increments starting from the LSB.
  Each increment is converted to a byte and added to the final list.
  
  \begin{code}
  integerToBytes :: Integer -> ByteCount -> [Word8]
  integerToBytes _ 0 = []
  integerToBytes x n = bytes where
      bytes = integerToBytes (x `shiftR` 8) (n - 1) ++ [step]
      step = fromIntegral x .&. 0xFF
  \end{code}
  
  Finally, the raw parsing is wrapped up in {\tt Put}. The second parameter
  allows the same code paths to be used for little- and big-endian
  serialisation.
  
  \begin{code}
  putFloat :: (RealFloat a) => ByteCount -> ([Word8] -> [Word8]) -> a -> Put
  putFloat width f x = putByteString $ B.pack bytes where
      bytes = f . rawToBytes . splitFloat width $ x
  \end{code}
  
  \section{Raw float components}
  
  Information about the raw bit patterns in the float is stored in
  {\tt RawFloat}, so they don't have to be passed around to the various
  format cases. The exponent should be biased, and the significand
  shouldn't have it's implied MSB (if applicable).
  
  \begin{code}
  data RawFloat = RawFloat
      { rawWidth            :: ByteCount
      , rawSign             :: Sign
      , rawExponent         :: Exponent
      , rawSignificand      :: Significand
      , rawExponentWidth    :: BitCount
      , rawSignificandWidth :: BitCount
      }
      deriving (Show)
  \end{code}
  
  \section{Exponents}
  
  Calculate the proper size of the exponent field, in bits, given the
  size of the full structure.
  
  \begin{code}
  exponentWidth :: BitCount -> BitCount
  exponentWidth k
      | k == 16         = 5
      | k == 32         = 8
      | k `mod` 32 == 0 = ceiling (4 * logBase 2 (fromIntegral k)) - 13
      | otherwise       = error "Invalid length of floating-point value"
  \end{code}
  
  \begin{code}
  bias :: Exponent -> BitCount -> Exponent
  bias e eWidth = e - (1 - (2 `pow` (eWidth - 1)))
  \end{code}
  
  \begin{code}
  unbias :: Exponent -> BitCount -> Exponent
  unbias e eWidth = e + 1 - (2 `pow` (eWidth - 1))
  \end{code}
  
  \section{Byte and bit counting}
  
  \begin{code}
  data Sign = Positive | Negative
      deriving (Show)
  
  newtype Exponent = Exponent Int
      deriving (Show, Eq, Num, Ord, Real, Enum, Integral, Bits)
  
  newtype Significand = Significand Integer
      deriving (Show, Eq, Num, Ord, Real, Enum, Integral, Bits)
  
  newtype BitCount = BitCount Int
      deriving (Show, Eq, Num, Ord, Real, Enum, Integral)
  
  newtype ByteCount = ByteCount Int
      deriving (Show, Eq, Num, Ord, Real, Enum, Integral)
  
  bitCount :: ByteCount -> BitCount
  bitCount (ByteCount x) = BitCount (x * 8)
  
  bitsInWord8 :: [Word8] -> BitCount
  bitsInWord8 = bitCount . ByteCount . length
  
  bitShiftL :: (Bits a) => a -> BitCount -> a
  bitShiftL x (BitCount n) = shiftL x n
  
  bitShiftR :: (Bits a) => a -> BitCount -> a
  bitShiftR x (BitCount n) = shiftR x n
  \end{code}
  
  \section{Utility}
  
  Considering a byte list as a sequence of bits, slice it from start
  inclusive to end exclusive, and return the resulting bit sequence as an
  integer.
  
  \begin{code}
  bitSlice :: [Word8] -> BitCount -> BitCount -> Integer
  bitSlice bs = sliceInt (foldl' step 0 bs) bitCount' where
      step acc w     = shiftL acc 8 + fromIntegral w
      bitCount'      = bitsInWord8 bs
  \end{code}
  
  Slice a single integer by start and end bit location
  
  \begin{code}
  sliceInt :: Integer -> BitCount -> BitCount -> BitCount -> Integer
  sliceInt x xBitCount s e = fromIntegral sliced where
      sliced = (x .&. startMask) `bitShiftR` (xBitCount - e)
      startMask = n1Bits (xBitCount - s)
      n1Bits n  = (2 `pow` n) - 1
  \end{code}
  
  Integral version of {\tt (**)}
  
  \begin{code}
  pow :: (Integral a, Integral b, Integral c) => a -> b -> c
  pow b e = floor $ fromIntegral b ** fromIntegral e
  \end{code}
  
  Detect whether a float is {\tt $-$NaN}
  
  \begin{code}
  isNegativeNaN :: RealFloat a => a -> Bool
  isNegativeNaN x = isNaN x && (floor x > 0)
  \end{code}
  