packages feed

FixedPoint-simple 0.4.2 → 0.5

raw patch · 3 files changed

+135/−32 lines, 3 filesdep +deepseqdep +template-haskell

Dependencies added: deepseq, template-haskell

Files

Data/FixedPoint.lhs view
@@ -31,18 +31,23 @@ >       , Int8192 >       -- * Big Word Types >       , Word128(..)+>       , Word72 >       , Word256 >       , Word512 >       , Word576+>       , Word584 >       , Word1024 >       , Word1280 >       , Word2048+>       , Word2632 >       , Word4096 >       , Word8192 >       -- * Type Constructors >       , GenericFixedPoint(..) >       , BigInt(..) >       , BigWord(..)+>       -- * Helpers+>       , fromInternal, toInternal, fracBits >       ) where > import Data.Bits > import Data.Word@@ -50,6 +55,7 @@ > import Foreign.Storable > import Foreign.Ptr > import Numeric+> import Control.DeepSeq  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@@ -60,7 +66,7 @@ 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. -> -- | GenericFixedPoitn is a type constructor for arbitrarily-sized fixed point+> -- | GenericFixedPoint is a type constructor for arbitrarily-sized fixed point > -- tyes. Take note the first type variable, @flat@, should be a signed int > -- equal to the size of the fixed point integral plus fractional bits. > -- The second type variable, @internal@, should be unsigned and twice@@ -69,10 +75,9 @@ > -- fractional bits in the fixed point type.  See the existing type aliases, > -- such as @FixedPoint4816@, for examples. -> data GenericFixedPoint flat internal fracBitRepr = FixedPoint flat+> data GenericFixedPoint flat internal fracBitRepr = FixedPoint { toFlat :: flat } >               deriving (Eq, Ord) >-> toFlat (FixedPoint x) = x > fromFlat = FixedPoint > > type FixedPoint20482048 = GenericFixedPoint Int4096 Word8192 Word2048@@ -87,9 +92,11 @@ > toInternal :: (Integral a, Num b) => GenericFixedPoint a b c -> b > toInternal (FixedPoint a) = fromIntegral a >+> -- Convert a fixed point into its internal (high precision) representation for use in computations. > fromInternal :: (Integral b, Num a) =>  b -> GenericFixedPoint a b c > fromInternal w = FixedPoint (fromIntegral w) >+> -- | Obtain the number of bits used to represent the fractional component of this fixed point. > fracBits :: (Bits c) => GenericFixedPoint a b c -> Int > fracBits = bitSize . getC >@@ -164,30 +171,47 @@ >         GenericFixedPoint a b c > pi' = realToFrac pi -> -- | The square root operation uses Newton's method but is parameterized by the number-> -- of iterations and stops early if we have arrived at a fixed point (no pun intended).-> -- Suggested iterations: 500 (But it increases with the size of the input!)-> 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)++> -- | The square root operation converges in O(bitSize input).+> sqrt' :: (Ord b, Integral b, Bits b, Integral a, Num a, Bits a, Bits c) =>+>          GenericFixedPoint a b c -> GenericFixedPoint a b c+> sqrt' x = fromInternal . fpSqrtRaw .  toInternal $ x >  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 (i-1) gNew+> -- Note: Using 'internal' instead of an unsigned version of 'flat'+> -- is unnecessary but preferable to adding yet another type variable or a type function.+>  fpSqrtRaw n0 | n0 < 0  = error "fpSqrt of a negative number"+>  fpSqrtRaw n0 = case iterate (step (bitSize n0)) (0,0,1,n0) !! bitSize n0 of+>                  (_,a,_,_) -> a `shiftR` ((bitSize n0 - fracBits x) `div` 2)+>  step sz (!s,!a,!t,!n) =+>    let s0 = (s `shiftL` 2) .|. (n `shiftR` (sz - 2))+>        n1 = n `shiftL` 2+>        a0 = a `shiftL` 1+>        (a1,s1,t0) = if s0 < t then (a0,       s0    , t - 1)+>                               else (a0 .|. 1, s0 - t, t + 1)+>        t1 = (t0 `shiftL` 1) .|. 1+>    in (s1,a1,t1,n1) -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+The below exp function is a lookup table augmented with a taylor series.+The taylor functuion 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.+e^x = (e^(x/2))^2 to break the problem down (if the table hasn't already). -This could probably be improved using a lookup table.+> expTable :: [(Double,Double)]+> expTable = let ys = [b*2**x | b <- [-1,1], x <- [8,7..(-10)]] in zip ys (map exp ys)+>+> exp' :: (Show a, Ord a, Fractional a, Eq a) => Int -> a -> a+> exp' n a =+>       let table = map (\(a,b) -> (realToFrac a, realToFrac b)) expTable+>       in (\(r,x) -> r * expTaylor n x) $ foldl op (1,a) table+>  where+>  op (r,x) (p,ep)+>     | signum p == signum x && abs p <= abs x = (r*ep, x - p)+>     | otherwise = (r,x) -> 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+> expTaylor :: (Ord a, Fractional a, Eq a) => Int -> a -> a+> expTaylor 0 a = 1+> expTaylor n a+>    | not (a > (-1) && a < 1) = let t = expTaylor n (a/2) in t*t >    | otherwise               = go 1 1 a >  where >  go !i !total !term@@ -199,7 +223,7 @@ 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' :: (Show a, Eq a, Ord a, Num a, Fractional a) => Int -> a -> a > erf' n x = >       let a1 = 0.254829592 >           a2 =  -0.284496736@@ -210,7 +234,7 @@ >           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');+>           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 @@ -290,7 +314,7 @@ >               | 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@@ -338,6 +362,9 @@  Larger word aliases follow. +> -- |A 72 bit unsigned word+> type Word72 = BigWord Word8 Word64+> > -- |A 256 bit unsigned word > type Word256 = BigWord Word128 Word128 >@@ -347,6 +374,9 @@ > -- |A 576 bit unsigned word > type Word576 = BigWord Word512 Word64 >+> -- |A 584 bit unsigned word+> type Word584 = BigWord Word72 Word512+> > -- |A 1024 bit unsigned word > type Word1024 = BigWord Word512 Word512 >@@ -356,6 +386,9 @@ > -- |A 2048 bit unsigned word > type Word2048 = BigWord Word1024 Word1024 >+> -- |A 2632 bit unsigned word+> type Word2632 = BigWord Word584 Word2048+> > -- |A 4096 bit unsigned word > type Word4096 = BigWord Word2048 Word2048 >@@ -435,14 +468,10 @@ >                  . (`shiftL` i) >                  . (id :: Integer -> Integer) >                  . fromIntegral $ b->               -- | i > bitSize l = shiftL (BigWord (fromIntegral l) 0) (i - bitSize l)->               -- | otherwise     = BigWord ((h `shiftL` i) .|. (fromIntegral (l `shiftR` (bitSize l - i)))) (l `shiftL` i) >       shiftR b i = fromIntegral >                  . (`shiftR` i) >                  . (id :: Integer -> Integer) >                  . fromIntegral $ b->               -- | i > bitSize h = shiftR (BigWord 0 h) (i - bitSize h)->               -- | otherwise     = BigWord (h `shiftR` i) ((l `shiftR` i) .|. fromIntegral (h `shiftL` (bitSize h - i))) >       isSigned _ = False >       testBit (BigWord h l) i >               | i >= bitSize l = testBit h (i - bitSize l)@@ -639,6 +668,18 @@ >       alignment ~(BigInt a) = alignment a >       peekElemOff ptr i = fmap BigInt (peekElemOff (castPtr ptr) i) >       pokeElemOff ptr i (BigInt a) = pokeElemOff (castPtr ptr) i a+>+> instance NFData a => NFData (BigInt a) where+>    rnf (BigInt a) = rnf a+>+> instance (NFData a, NFData b) => NFData (BigWord a b) where+>    rnf (BigWord a b) = rnf a `seq` rnf b+>+> instance NFData Word128 where+>    rnf (W128 a b) = rnf a `seq` rnf b+>+> instance NFData flat => NFData (GenericFixedPoint flat i r) where+>    rnf (FixedPoint f) = rnf f > > instance (Storable a, Storable b) => Storable (BigWord a b) where >       sizeOf ~(BigWord a b) = sizeOf a + sizeOf b
+ Data/FixedPoint/TH.hs view
@@ -0,0 +1,62 @@+module Data.FixedPoint.TH +    ( mkWord+    , mkInt+    , mkFixedPoint+    ) where++import Language.Haskell.TH+import Data.Maybe++-- |@$(mkWord X)@ Makes a type alias named @WordX@ for a word of @X@ bits.+-- Notice @X@ must be a multiple of 8, 'Data.Word.Word8' must be in scope,+-- 'Data.FixedPoint.BigWord' must be in scope, and this splice will add+-- all smaller @WordY@ type aliases needed that aren't already in scope.+mkWord :: Int -> DecsQ+mkWord i+  | i `rem` 8 /= 0 = error ("Can not build a word of bit size " ++ show i)+  | otherwise = do+        info <- lookupTypeName (mkS i)+        let b = isNothing info+        if b then do+                let (h,l) = getParts i+                hD <- mkWord h+                lD <- mkWord l+                a <- tySynD (mkW i) [] (appT (appT (conT $ mkName "BigWord") (conT $ mkW h)) (conT $ mkW l))+                return $ a:(hD++lD)+             else return []++mkS :: Int -> String+mkS = ("Word" ++) . show++mkW,mkI :: Int -> Name+mkW = mkName . mkS++mkI = mkName . ("Int" ++) . show++getParts i =+    let l = 2^(floor (logBase 2 (fromIntegral i)))+        h = i - l+    in (h,l)++-- |@$(mkInt X)@ Makes a type alias named @IntX@ for an int of X bits.+-- See the requirements under 'mkWord' for additional information.+mkInt :: Int -> DecsQ+mkInt i = do+    d <- mkWord i+    e <- tySynD (mkName . ("Int" ++) . show $ i) [] (appT (conT $ mkName "BigInt") (conT $ mkW i))+    return (e:d)++-- @mkFixedPoint X Y@ Builds a fixed point alias named @FixedPointX_Y@. See+-- the requirements under 'mkWord' for additional information.+mkFixedPoint :: Int -> Int -> DecsQ+mkFixedPoint int frac+  | (int + frac) `rem` 8 /= 0 = error "For fixed points, The sum of the integral and fractional bits must be a multiple of 8."+  | frac `rem` 8 /= 0 = error "For fixed points, the fractional representation must be a multiple of 8."+  | otherwise = do+      let flat = int + frac+      f <- mkInt flat+      i <- mkWord (flat*2)+      r <- mkWord frac+      x <- tySynD (mkName $ "FixedPoint" ++ show int ++ "_" ++ show frac)+                  [] (appT (appT (appT (conT $ mkName "GenericFixedPoint") (conT $ mkI flat)) (conT $ mkW $ flat*2)) (conT $ mkW frac))+      return (x : r ++ i ++ f)
FixedPoint-simple.cabal view
@@ -1,5 +1,5 @@ Name:                FixedPoint-simple-Version:             0.4.2+Version:             0.5 Synopsis:            Fixed point, large word, and large int numerical representations (types and common class instances) Description:         This library uses elementary techniques to implement fixed point types in terms                      of basic integrals such as Word64.  All mathematical operations are implemented@@ -20,8 +20,8 @@   Library-  Exposed-modules:     Data.FixedPoint-  Build-depends:       base >= 4 && < 5+  Exposed-modules:     Data.FixedPoint, Data.FixedPoint.TH+  Build-depends:       base >= 4 && < 5, deepseq, template-haskell >= 2.8   ghc-options:         -O2 -funbox-strict-fields   -- Other-modules:          -- Build-tools: