>
>
>
> module Data.FixedPoint
> (
> FixedPoint4816
> , FixedPoint3232
> , FixedPoint6464
> , FixedPoint128128
> , FixedPoint256256
> , FixedPoint512512
> , FixedPoint10241024
>
> , erf'
> , exp'
> , sqrt'
> , pi'
>
> , Int128
> , Int256
> , Int512
> , Int1024
> , Int2048
> , Int4096
> , Int8192
>
> , Word128(..)
> , Word256
> , Word512
> , Word1024
> , Word2048
> , Word4096
> , Word8192
>
> , GenericFixedPoint(..)
> , BigInt(..)
> , BigWord(..)
> ) where
> import Data.Bits
> import Data.Word
> import Data.Int
> import Debug.Trace
> import Numeric
This code implements n.m fixed point types allowing for a range from (2^(n-1),-2^(n-1)].
Given a type `GenericFixedPoint flat internal fracBitRepr` the values m and n
are:
m = bitSize fracBitRepr
n = bitSize flat - m
The 'Flat' representation is a signed n+m bit value. The 'Internal' representation should be a 2*(n+m)
unsigned value for use in division.
>
>
>
>
>
>
>
>
> data GenericFixedPoint flat internal fracBitRepr = FixedPoint flat
> deriving (Eq, Ord)
>
> toFlat (FixedPoint x) = x
> fromFlat = FixedPoint
>
> type FixedPoint20482048 = GenericFixedPoint Int4096 Word8192 Word2048
> type FixedPoint10241024 = GenericFixedPoint Int2048 Word4096 Word1024
> type FixedPoint512512 = GenericFixedPoint Int1024 Word2048 Word512
> type FixedPoint256256 = GenericFixedPoint Int512 Word1024 Word256
> type FixedPoint128128 = GenericFixedPoint Int256 Word512 Word128
> type FixedPoint6464 = GenericFixedPoint Int128 Word256 Word64
> type FixedPoint3232 = GenericFixedPoint Int64 Word128 Word32
> type FixedPoint4816 = GenericFixedPoint Int64 Word128 Word16
>
> toInternal :: (Integral a, Num b) => GenericFixedPoint a b c -> b
> toInternal (FixedPoint a) = fromIntegral a
>
> fromInternal :: (Integral b, Num a) => b -> GenericFixedPoint a b c
> fromInternal w = FixedPoint (fromIntegral w)
>
> fracBits :: (Bits c) => GenericFixedPoint a b c -> Int
> fracBits = bitSize . getC
>
> getC :: GenericFixedPoint a b c -> c
> getC = const undefined
>
> instance (Integral a, Integral b, Bits a, Bits b, Bits c) =>
> Show (GenericFixedPoint a b c) where
> show = (show :: Double -> String) . realToFrac
>
> instance (Enum a, Num a, Bits a, Bits c) =>
> Enum (GenericFixedPoint a b c) where
> succ (FixedPoint a) = FixedPoint (a + 1)
> pred (FixedPoint a) = FixedPoint (a 1)
> fromEnum f@(FixedPoint a) = fromEnum (a `shiftR` fracBits f)
> toEnum x = let r = FixedPoint (toEnum x `shiftL` fracBits r) in r
>
> xOR a b = a && not b || not a && b
>
> instance (Ord a, Num a, Bits a, Bits b, Integral a, Integral b, Bits c) =>
> Num (GenericFixedPoint a b c) where
>
>
>
>
> (FixedPoint a) + (FixedPoint b) = FixedPoint (a + b)
> aval@(FixedPoint a) * bval@(FixedPoint b) =
> let w = (toInternal $ abs aval) * (toInternal $ abs bval)
> op = if xOR (aval < 0) (bval < 0) then negate else id
> in op $ fromInternal (w `shiftR` fracBits aval)
> a b = a + negate b
> negate (FixedPoint a) = FixedPoint (negate a)
> signum (FixedPoint a) = FixedPoint $ signum a
> abs (FixedPoint a) = FixedPoint (abs a)
> fromInteger i = let r = FixedPoint (fromInteger i `shiftL` fracBits r) in r
>
> instance (Ord a, Integral a, Bits a, Num a, Bits b, Integral b, Bits c) =>
> Fractional (GenericFixedPoint a b c) where
> aval / bval =
> let wa = toInternal $ abs aval
> wb = toInternal $ abs bval
> wr = ((wa`shiftL` fracBits aval) `div` wb)
> op = if xOR (aval < 0) (bval < 0) then negate else id
> in (op $ fromInternal wr)
>
> fromRational a =
> let (r,k) = properFraction a
> (rf,_) = properFraction (abs $ k * 2^(fracBits res))
> signFix = if a < 0 then negate else id
> res = FixedPoint ((abs r `shiftL` fracBits res) .|. rf)
> in signFix res
>
> instance (Integral a, Ord a, Num a, Bits a, Bits b, Integral b, Bits c) =>
> Real (GenericFixedPoint a b c) where
> toRational f@(FixedPoint a)
> | a < 0 = negate (toRational $ negate f)
> | otherwise = fromIntegral a / (2^(fracBits f))
>
> instance (Integral a, Bits a, Integral b, Num a, Bits b, Bits c) =>
> Read (GenericFixedPoint a b c) where
> readsPrec n s = [ (realToFrac (r::Double), s) | (r,s) <- readsPrec n s]
Now we need some advanced functions (beyond +,-,*,/) on our fixed point type.
Specifically, we want 'exp' (exponentiation with a base of 'e' ~ 2.71), erf (the
"error function"), and square root. All of these will be implemented by some
form of approximation such as a Taylor series. We thus parameterize the number
of terms to allow testing / user control over cost and accuracy.
> pi' :: (Integral a, Bits a, Integral b, Num a, Bits b, Bits c) =>
> GenericFixedPoint a b c
> pi' = realToFrac pi
>
>
>
> sqrt' :: (Eq a, Integral a, Num a, Bits a, Integral b, Bits b, Bits c) =>
> Int -> GenericFixedPoint a b c -> GenericFixedPoint a b c
> sqrt' cnt input = fromFlat (go cnt 1) `shiftL` (fracBits input `div` 2)
> where
> a = toFlat input
> go 0 g = g
> go i g =
> let gNew = ((a`div`g) + g) `div` 2
> in if gNew == g then g else go (i1) gNew
The below exp function includes a taylor series (the 'go' function) but that
operation alone looses precision too quickly so we restrict it's use to an
acceptable range. Outside of that range we depend on the property
e^x = (e^(x/2))^2 to break the problem down.
This could probably be improved using a lookup table.
> exp' :: (Ord a, Fractional a, Eq a) => Int -> a -> a
> exp' 0 a = 1
> exp' n a
> | not (a > (1) && a < 1) = let t = exp' n (a/2) in t*t
> | otherwise = go 1 1 a
> where
> go !i !total !term
> | i <= n = let iNew = i + 1 in go iNew (total + term) (term*a/fromIntegral iNew)
> | otherwise = total
The first error function was a direct Taylor series implementation. It lost
accuracy too quickly due to the factorial term. A superior version from
picomath.org uses precomputed values (picomath released the code as public
domain).
> erf' :: (Eq a, Ord a, Num a, Fractional a) => Int -> a -> a
> erf' n x =
> let a1 = 0.254829592
> a2 = 0.284496736
> a3 = 1.421413741
> a4 = 1.453152027
> a5 = 1.061405429
> p = 0.3275911
> sign = if x < 0 then (1) else 1
> x' = abs x
> t = 1 / (1 + p * x')
> y = 1 (((((a5*t + a4)*t) + a3)*t + a2)*t + a1) * t * exp' n (x'*x');
> in sign * y
>
> flat1 :: (a -> Int -> a ) -> GenericFixedPoint a b c -> Int
> -> GenericFixedPoint a b c
> flat1 op a i = fromFlat . flip op i . toFlat $ a
> flat2 :: (a -> a -> a) -> GenericFixedPoint a b c
> -> GenericFixedPoint a b c -> GenericFixedPoint a b c
> flat2 op a b = fromFlat (op (toFlat a) (toFlat b))
>
> instance (Ord a, Bits a, Bits b, Integral a, Integral b, Bits c) =>
> Bits (GenericFixedPoint a b c) where
> setBit = flat1 setBit
> (.|.) = flat2 (.|.)
> xor = flat2 xor
> (.&.) = flat2 (.&.)
> complement = fromFlat . complement . toFlat
> bitSize a = fracBits a * 2
> isSigned _ = False
> shiftL = flat1 shiftL
> shiftR = flat1 shiftR
To implement FixedPoint division without loosing precision large words are
needed to represent our internal state.
Unlike other word sizes, Word128 is explicitly implemented and unpacked. It is
suspected that this will result in a performance difference. Benchmark results
would be interesting. If this proves to be a big win then implementing the
other Word sizes by generating code using TH, instead of using overloaded type
classes, would be a beneficial task.
> data Word128 = W128 !Word64 !Word64
> instance Num Word128 where
> W128 ah al + W128 bh bl =
> let rl = al + bl
> rh = ah + bh + if rl < al then 1 else 0
> in W128 rh rl
> W128 ah al W128 bh bl =
> let rl = al bl
> rh = ah bh if rl > al then 1 else 0
> in W128 rh rl
> a * b = go 0 0
> where
> go 64 r = r
> go i r
> | testBit b i = go (i+1) (r + (a `shiftL` i))
> | otherwise = go (i+1) r
> negate a = a
> abs a = a
> signum a = if a > 0 then 1 else 0
> fromInteger i = W128 (fromIntegral $ i `shiftR` 64) (fromIntegral i)
>
> pointwise op (W128 a b) = W128 (op a) (op b)
> pointwise2 op (W128 a b) (W128 c d) = W128 (op a c) (op b d)
>
> instance Eq Word128 where
> a == b = EQ == compare a b
>
> instance Bits Word128 where
> bit i | i >= 64 = W128 (bit $ i 64) 0
> | otherwise = W128 0 (bit i)
> complement = pointwise complement
> (.&.) = pointwise2 (.&.)
> (.|.) = pointwise2 (.|.)
> xor = pointwise2 xor
> setBit (W128 h l) i
> | i >= 64 = W128 (setBit h (i 64)) l
> | otherwise = W128 h (setBit l i)
> shiftL (W128 h l) i
> | i > bitSize l = shiftL (W128 l 0) (i bitSize l)
> | otherwise = W128 ((h `shiftL` i) .|. (l `shiftR` (bitSize l i))) (l `shiftL` i)
> shiftR (W128 h l) i
> | i > bitSize h = shiftR (W128 0 h) (i bitSize h)
> | otherwise = W128 (h `shiftR` i) ((l `shiftR` i) .|. h `shiftL` (bitSize h i))
> isSigned _ = False
> testBit (W128 h l) i
> | i >= bitSize l = testBit h (i bitSize l)
> | otherwise = testBit l i
> bitSize _ = 128
>
> instance Enum Word128 where
> toEnum i = W128 0 (toEnum i)
> fromEnum (W128 _ l) = fromEnum l
> pred (W128 h 0) = W128 (pred h) maxBound
> pred (W128 h l) = W128 h (pred l)
> succ (W128 h l) = if l == maxBound then W128 (succ h) 0 else W128 h (succ l)
>
> instance Ord Word128 where
> compare (W128 ah al) (W128 bh bl) = compare (ah,al) (bh,bl)
>
> instance Real Word128 where
> toRational w = toRational (fromIntegral w :: Integer)
>
> instance Integral Word128 where
> toInteger (W128 h l) = (fromIntegral h `shiftL` bitSize l) + fromIntegral l
> divMod = quotRem
> quotRem a@(W128 ah al) b@(W128 bh bl) =
> let r = a q*b
> q = go 0 (bitSize a) 0
> in (q,r)
> where
>
> go :: Word128 -> Int -> Word128 -> Word128
> go t 0 v = if v >= b then t+1 else t
> go t i v
> | v >= b = go (setBit t i) i' v2
> | otherwise = go t i' v1
> where
> i' = i 1
> newBit = if (testBit a i') then 1 else 0
> v1 = (v `shiftL` 1) .|. newBit
> v2 = ((v b) `shiftL` 1) .|. newBit
>
> instance Show Word128 where
> show = show . fromIntegral
>
> instance Read Word128 where
> readsPrec i s = let readsPrecI :: Int -> ReadS Integer
> readsPrecI = readsPrec
> in [(fromIntegral i, str) | (i,str) <- readsPrecI i s]
>
> instance Bounded Word128 where
> maxBound = W128 maxBound maxBound
> minBound = W128 minBound minBound
Larger word aliases follow.
>
> type Word256 = BigWord Word128
>
>
> type Word512 = BigWord Word256
>
>
> type Word1024 = BigWord Word512
>
>
> type Word2048 = BigWord Word1024
>
>
> type Word4096 = BigWord Word2048
>
>
> type Word8192 = BigWord Word4096
>
>
>
>
> data BigWord a = BigWord a a
> instance (Bits a, Num a, Ord a) => Num (BigWord a) where
> BigWord ah al + BigWord bh bl =
> let rl = al + bl
> rh = ah + bh + if rl < al then 1 else 0
> in BigWord rh rl
> BigWord ah al BigWord bh bl =
> let rl = al bl
> rh = ah bh if rl > al then 1 else 0
> in BigWord rh rl
> a * b = go 0 0
> where
> go i r
> | i == bitSize r = r
> | testBit b i = go (i+1) (r + (a `shiftL` i))
> | otherwise = go (i+1)r
> negate a = a
> abs a = a
> signum a = if a > 0 then 1 else 0
> fromInteger i =
> let r@(BigWord _ b) = BigWord (fromIntegral $ i `shiftR` (bitSize b)) (fromIntegral i)
> in r
>
> pointwiseBW op (BigWord a b) = BigWord (op a) (op b)
> pointwiseBW2 op (BigWord a b) (BigWord c d) = BigWord (op a c) (op b d)
>
> instance (Ord a) => Eq (BigWord a) where
> a == b = EQ == compare a b
>
> instance (Ord a, Bits a) => Bits (BigWord a) where
> bit i | i >= bitSize b = r1
> | otherwise = r2
> where r1@(BigWord _ b) = BigWord (bit $ i bitSize b) 0
> r2 = BigWord 0 (bit i)
> complement = pointwiseBW complement
> (.&.) = pointwiseBW2 (.&.)
> (.|.) = pointwiseBW2 (.|.)
> xor = pointwiseBW2 xor
> setBit (BigWord h l) i
> | i >= bitSize l = BigWord (setBit h (ibitSize l)) l
> | otherwise = BigWord h (setBit l i)
> shiftL (BigWord h l) i
> | i > bitSize l = shiftL (BigWord l 0) (i bitSize l)
> | otherwise = BigWord ((h `shiftL` i) .|. (l `shiftR` (bitSize l i))) (l `shiftL` i)
> shiftR (BigWord h l) i
> | i > bitSize h = shiftR (BigWord 0 h) (i bitSize h)
> | otherwise = BigWord (h `shiftR` i) ((l `shiftR` i) .|. h `shiftL` (bitSize h i))
> isSigned _ = False
> testBit (BigWord h l) i
> | i >= bitSize l = testBit h (i bitSize l)
> | otherwise = testBit l i
> bitSize ~(BigWord h l) = bitSize h + bitSize l
>
> instance (Bounded a,Eq a,Num a, Enum a) => Enum (BigWord a) where
> toEnum i = BigWord 0 (toEnum i)
> fromEnum (BigWord _ l) = fromEnum l
> pred (BigWord h 0) = BigWord (pred h) maxBound
> pred (BigWord h l) = BigWord h (pred l)
> succ (BigWord h l) = if l == maxBound then BigWord (succ h) 0 else BigWord h (succ l)
>
> instance Bounded a => Bounded (BigWord a) where
> maxBound = BigWord maxBound maxBound
> minBound = BigWord minBound minBound
>
> instance Ord a => Ord (BigWord a) where
> compare (BigWord a b) (BigWord c d) = compare (a,b) (c,d)
>
> instance (Bits a, Real a, Bounded a, Integral a) => Real (BigWord a) where
> toRational w = toRational (fromIntegral w :: Integer)
>
> instance (Bounded a, Integral a, Bits a) => Integral (BigWord a) where
> toInteger (BigWord h l) = (fromIntegral h `shiftL` bitSize l) + fromIntegral l
> divMod = quotRem
> quotRem a b =
> let r = a q * b
> q = go 0 (bitSize a) 0
> in (q, r)
> where
>
> go t 0 v = if v >= b then t + 1 else t
> go t i v
> | v >= b = go (setBit t i) i' v2
> | otherwise = go t i' v1
> where
> i' = i 1
> newBit = if (testBit a i') then 1 else 0
> v1 = (v `shiftL` 1) .|. newBit
> v2 = ((vb) `shiftL` 1) .|. newBit
>
> instance (Bounded a, Bits a, Integral a) => Show (BigWord a) where
> show = show . fromIntegral
>
> instance (Num a, Bits a, Ord a) => Read (BigWord a) where
> readsPrec i s = let readsPrecI :: Int -> ReadS Integer
> readsPrecI = readsPrec
> in [(fromIntegral i, str) | (i,str) <- readsPrecI i s]
>
For fixed point, the flat representation needs to be signed.
>
> type Int128 = BigInt Word128
>
> type Int256 = BigInt Word256
>
> type Int512 = BigInt Word512
>
> type Int1024 = BigInt Word1024
>
> type Int2048 = BigInt Word2048
>
> type Int4096 = BigInt Word4096
>
> type Int8192 = BigInt Word8192
>
>
>
>
>
> newtype BigInt a = BigInt { unBI :: a }
>
> instance (Ord a, Bits a) => Ord (BigInt a) where
> compare (BigInt a) (BigInt b)
> | testBit a (bitSize a 1) = if testBit b (bitSize b 1)
> then compare a b
> else LT
> | testBit b (bitSize b 1) = GT
> | otherwise = compare a b
>
> instance (Eq a) => Eq (BigInt a) where
> BigInt a == BigInt b = a == b
>
> instance (Show a, Num a, Bits a, Ord a) => Show (BigInt a) where
> show i@(BigInt a)
> | i < 0 = '-' : show (complement a + 1)
> | otherwise = show a
>
> instance (Num a, Bits a, Ord a) => Read (BigInt a) where
> readsPrec i s = let readsPrecI :: Int -> ReadS Integer
> readsPrecI = readsPrec
> in [(fromIntegral i, str) | (i,str) <- readsPrecI i s]
>
> instance (Num a, Bits a, Ord a) => Num (BigInt a) where
> (BigInt a) + (BigInt b) = BigInt (a+b)
> (BigInt a) (BigInt b) = BigInt (ab)
> (BigInt a) * (BigInt b) = BigInt (a*b)
> negate (BigInt a) = BigInt (complement a + 1)
> signum a = if a < 0 then 1 else if a > 0 then 1 else 0
> abs a = if a < 0 then negate a else a
> fromInteger i = if i < 0 then negate (BigInt $ fromInteger (abs i))
> else BigInt (fromInteger i)
>
> instance (Bits a, Ord a) => Bits (BigInt a) where
> (.&.) a b = BigInt (unBI a .&. unBI b)
> (.|.) a b = BigInt (unBI a .|. unBI b)
> xor a b = BigInt (unBI a `xor` unBI b)
> complement = BigInt . complement . unBI
> shiftL a i = BigInt . (`shiftL` i) . unBI $ a
> shiftR a i = (if a < 0 then \x -> foldl setBit x [bitSize a1, bitSize a 2 .. bitSize a i]
> else id)
> . BigInt
> . (`shiftR` i)
> . unBI
> $ a
> bit = BigInt . bit
> setBit a i = BigInt . (`setBit` i) . unBI $ a
> testBit a i = (`testBit` i) . unBI $ a
> bitSize (BigInt a) = bitSize a
> isSigned _ = True
>
> instance (Bits a, Ord a, Integral a, Bounded a, Num a) => Enum (BigInt a) where
> toEnum i = fromIntegral i
> fromEnum i = fromIntegral i
> pred a | a > minBound = (a 1)
> succ a | a < maxBound = (a + 1)
>
> instance (Integral a, Bits a, Bounded a) => Integral (BigInt a) where
> toInteger i@(BigInt h) = (if i < 0 then negate else id) (toInteger h)
> quotRem a b =
> let (BigInt ah) = abs a
> (BigInt bh) = abs b
> (q1,r1) = quotRem ah bh
> in if a < 0 && b < 0
> then (BigInt q1, negate $ BigInt r1)
> else if a < 0
> then (negate $ BigInt q1, negate $ BigInt r1)
> else if b < 0
> then (negate $ BigInt q1, BigInt r1)
> else (BigInt q1, BigInt r1)
>
> instance (Real a, Bounded a, Integral a, Bits a) => Real (BigInt a) where
> toRational = fromIntegral
>
>
> instance (Bounded a, Ord a, Bits a) => Bounded (BigInt a) where
> minBound = let r = fromIntegral (negate (2^ (bitSize r 1))) in r
> maxBound = let r = fromIntegral (2^(bitSize r 1) 1) in r