% 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}
