diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,25 @@ # Revision history for libBF-hs +## 0.6 -- 2021-01-29 + +* Fix a bug with `frint` and `bfRoundInt` that was causing incorrect +rounding modes to be selected. + +* Implement additional operations on `BigFloat` values, including +some missing predicates, conversions to and from binary form +and fused-multiply-add. + +* Fix a compile issue on 32-bit systems. + +## 0.5.2 -- 2021-01-12 + +* Use `install-includes` over `c-sources` for header files to avoid linker + issues. + +## 0.5.1 -- 2020-07-13 + +* Add header files to `c-sources` field to include them in `sdist`. + ## 0.5.0 -- 2020-07-01 * First version. Released on an unsuspecting world. diff --git a/libBF.cabal b/libBF.cabal --- a/libBF.cabal +++ b/libBF.cabal @@ -1,7 +1,7 @@ cabal-version: 2.2 name: libBF -version: 0.5.1 +version: 0.6 synopsis: A binding to the libBF library. description: LibBF is a C library for working with arbitray precision IEEE 754 floating point numbers. @@ -18,7 +18,9 @@ type: git location: https://github.com/GaloisInc/libBF-hs.git - +flag system-libbf + default: False + description: Use system libbf instead library exposed-modules: @@ -28,25 +30,33 @@ build-depends: base >=4.12.0.0 && < 5, - deepseq + deepseq, + hashable >= 1.3 build-tool-depends: hsc2hs:hsc2hs hs-source-dirs: src - include-dirs: - libbf-2020-01-19 + if flag(system-libbf) + extra-libraries: bf + c-sources: + cbits/libbf-hs.c + else + include-dirs: + libbf-2020-01-19 - includes: - libbf-2020-01-19/libbf.h + includes: + libbf-2020-01-19/libbf.h - c-sources: - libbf-2020-01-19/libbf.h - libbf-2020-01-19/cutils.h - libbf-2020-01-19/cutils.c - libbf-2020-01-19/libbf.c - cbits/libbf-hs.c + install-includes: + libbf-2020-01-19/libbf.h + libbf-2020-01-19/cutils.h + + c-sources: + libbf-2020-01-19/cutils.c + libbf-2020-01-19/libbf.c + cbits/libbf-hs.c ghc-options: -Wall default-language: Haskell2010 diff --git a/src/LibBF.hs b/src/LibBF.hs --- a/src/LibBF.hs +++ b/src/LibBF.hs @@ -1,3 +1,4 @@ +{-# Language BangPatterns #-} {-# Language BlockArguments #-} {-# Language Trustworthy #-} -- | Computation with high-precision floats. @@ -20,20 +21,27 @@ , bfToRep , BFRep(..) , BFNum(..) + , bfFromBits + , bfToBits -- * Predicates , bfIsFinite + , bfIsInf , bfIsZero , bfIsNaN + , bfIsNormal + , bfIsSubnormal , bfCompare , bfSign , bfExponent + , bfIsPos + , bfIsNeg , Sign(..) -- * Arithmetic - , bfNeg - , bfAdd, bfSub, bfMul, bfDiv - , bfMulWord, bfMulInt, bfMul2Exp + , bfNeg, bfAbs + , bfAdd, bfSub, bfMul, bfDiv, bfRem + , bfFMA, bfMulWord, bfMulInt, bfMul2Exp , bfSqrt , bfPow @@ -52,6 +60,8 @@ ) where +import Data.Bits +import Data.Hashable import Data.Word import Data.Int import System.IO.Unsafe @@ -60,7 +70,6 @@ import LibBF.Opts import Control.DeepSeq - -- | Arbitrary precision floating point numbers. newtype BigFloat = BigFloat BF @@ -141,6 +150,10 @@ BigFloat x <= BigFloat y = unsafe (cmpLEQ x y) +instance Hashable BigFloat where + hashWithSalt s x = hashWithSalt s (bfToRep x) + + {-| Compare the two numbers. The special values are ordered like this: * -0 < 0 @@ -153,7 +166,7 @@ bfCompare (BigFloat x) (BigFloat y) = unsafe (cmp x y) --- | Is this a "normal" (i.e., non-infinite, non NaN) number. +-- | Is this a finite (i.e., non-infinite, non NaN) number. bfIsFinite :: BigFloat -> Bool bfIsFinite (BigFloat x) = unsafe (isFinite x) @@ -161,10 +174,50 @@ bfIsNaN :: BigFloat -> Bool bfIsNaN (BigFloat x) = unsafe (M.isNaN x) --- | Get the sign of a number. Assumes the input is not NaN. +-- | Is this value infinite +bfIsInf :: BigFloat -> Bool +bfIsInf (BigFloat x) = unsafe (isInf x) + +-- | This is a "normal" number, which means it is not +-- a NaN, not a zero, not infinite, and not subnormal. +bfIsNormal :: BFOpts -> BigFloat -> Bool +bfIsNormal opts bf = + case bfToRep bf of + rep@(BFRep _sgn (Num _ _)) -> not (repIsSubnormal opts rep) + _ -> False + +-- | This number is "subnormal", which means it is among the smallest +-- representable numbers for the given precision and exponent bits. +-- These numbers differ from "normal" numbers in that they do not use +-- an implicit leading 1 bit in the binary representation. +bfIsSubnormal :: BFOpts -> BigFloat -> Bool +bfIsSubnormal opts bf = repIsSubnormal opts (bfToRep bf) + +-- | Get the sign of a number. Returns 'Nothing' if the number is `NaN`. bfSign :: BigFloat -> Maybe Sign bfSign (BigFloat x) = unsafe (getSign x) +-- | Compute the absolute value of a number. +bfAbs :: BigFloat -> BigFloat +bfAbs bf = + case bfSign bf of + Just Neg -> bfNeg bf + _ -> bf + +-- | Is this value positive +bfIsPos :: BigFloat -> Bool +bfIsPos bf = + case bfSign bf of + Just Pos -> True + _ -> False + +-- | Is this value negative +bfIsNeg :: BigFloat -> Bool +bfIsNeg bf = + case bfSign bf of + Just Neg -> True + _ -> False + -- | Get the exponent for the given number. -- Infinity, zero and NaN do not have an exponent. bfExponent :: BigFloat -> Maybe Int64 @@ -208,6 +261,15 @@ bfDiv :: BFOpts -> BigFloat -> BigFloat -> (BigFloat,Status) bfDiv opt (BigFloat x) (BigFloat y) = newBigFloat' (fdiv opt x y) +-- | Compute the remainder @x - y * n@ where @n@ is the integer +-- nearest to @x/y@ (with ties broken to even values of @n@). +bfRem :: BFOpts -> BigFloat -> BigFloat -> (BigFloat, Status) +bfRem opt (BigFloat x) (BigFloat y) = newBigFloat' (frem opt x y) + +-- | Compute the fused-multiply-add @(x*y)+z@ +bfFMA :: BFOpts -> BigFloat -> BigFloat -> BigFloat -> (BigFloat, Status) +bfFMA opt (BigFloat x) (BigFloat y) (BigFloat z) = newBigFloat' (ffma opt x y z) + -- | Square root of two numbers useing the given options. bfSqrt :: BFOpts -> BigFloat -> (BigFloat,Status) bfSqrt opt (BigFloat x) = newBigFloat' (fsqrt opt x) @@ -219,11 +281,11 @@ fround opt bf ) --- | Round to an integer using the given parameters. -bfRoundInt :: BFOpts -> BigFloat -> (BigFloat,Status) -bfRoundInt opt (BigFloat x) = newBigFloat' (\bf -> +-- | Round to an integer using the given rounding mode. +bfRoundInt :: RoundMode -> BigFloat -> (BigFloat,Status) +bfRoundInt r (BigFloat x) = newBigFloat' (\bf -> do setBF x bf - frint opt bf + frint r bf ) -- | Exponentiate a word to a positive integer power. @@ -268,6 +330,109 @@ -------------------------------------------------------------------------------- +-- | Make a float using "raw" bits representing the bitvector +-- representation of a floating-point value with the +-- exponent and precision bits given by the options. +bfFromBits :: + BFOpts -> + Integer {- ^ Raw bits -} -> + BigFloat +bfFromBits opts bits + | expoBiased == 0 && mant == 0 = -- zero + if isNeg then bfNegZero else bfPosZero + | expoBiased == eMask && mant == 0 = -- infinity + if isNeg then bfNegInf else bfPosInf + | expoBiased == eMask = bfNaN -- NaN + + | expoBiased == 0 = -- Subnormal + case bfMul2Exp opts' (bfFromInteger mant) (expoVal + 1) of + (num,Ok) -> if isNeg then bfNeg num else num + (_,s) -> error $ unwords ["bfFromBits", "subnormal case", "Unexpected status:", show s, show bits, show mant, show expoVal, show e, show p ] + + | otherwise = -- Normal + case bfMul2Exp opts' (bfFromInteger mantVal) expoVal of + (num,Ok) -> if isNeg then bfNeg num else num + (_,s) -> error $ unwords ["bfFromBits", "normal case", "Unexpected status:", show s, show bits, show mantVal, show expoVal, show e, show p ] + + where + e = getExpBits opts + p = getPrecBits opts + + opts' = opts <> allowSubnormal + + p' = fromIntegral p - 1 :: Int + eMask = (1 `shiftL` e) - 1 :: Int64 + pMask = (1 `shiftL` p') - 1 :: Integer + + isNeg = testBit bits (e + p') + + mant = pMask .&. bits :: Integer + mantVal = mant `setBit` p' :: Integer + -- accounts for the implicit 1 bit + + expoBiased = eMask .&. fromInteger (bits `shiftR` p') :: Int64 + bias = eMask `shiftR` 1 :: Int64 + expoVal = expoBiased - bias - fromIntegral p' :: Int64 + + +-- | Turn a float into raw bits. +-- @NaN@ is represented as a positive "quiet" @NaN@ +-- (most significant bit in the significand is set, the rest of it is 0). +bfToBits :: BFOpts -> BigFloat -> Integer +bfToBits opts bf = res + where + res = (isNeg `shiftL` (e+p')) + .|. (expBiased `shiftL` p') + .|. (mant `shiftL` 0) + + e = getExpBits opts + p = getPrecBits opts + + p' = fromIntegral p - 1 :: Int + + eMask = (1 `shiftL` e) - 1 :: Integer + pMask = (1 `shiftL` p') - 1 :: Integer + + (isNeg, expBiased, mant) = + case bfToRep bf of + BFNaN -> (0, eMask, 1 `shiftL` (p' - 1)) + BFRep s num -> (sign, be, ma) + where + sign = case s of + Neg -> 1 + Pos -> 0 + + (be,ma) = + case num of + Zero -> (0,0) + Num i ev + | ex <= 0 -> + (0, i `shiftL` (p'-m-1+fromInteger ex)) -- subnormal case + | otherwise -> + (ex, (i `shiftL` (p' - m)) .&. pMask) -- normal case + where + m = msb 0 i - 1 + bias = eMask `shiftR` 1 + ex = toInteger ev + bias + toInteger m + + Inf -> (eMask,0) + + msb !n j = if j == 0 then n else msb (n+1) (j `shiftR` 1) + +-- | test if a given big float representation is subnormal +repIsSubnormal :: BFOpts -> BFRep -> Bool +repIsSubnormal opts (BFRep _s (Num i ev)) = ex <= 0 + where + e = getExpBits opts + eMask = (1 `shiftL` e) - 1 :: Integer + bias = eMask `shiftR` 1 + + m = msb (0 :: Int) i - 1 + ex = toInteger ev + bias + toInteger m + + msb !n j = if j == 0 then n else msb (n+1) (j `shiftR` 1) + +repIsSubnormal _opts _rep = False diff --git a/src/LibBF/Mutable.hsc b/src/LibBF/Mutable.hsc --- a/src/LibBF/Mutable.hsc +++ b/src/LibBF/Mutable.hsc @@ -30,6 +30,7 @@ , getExp , isFinite + , isInf , LibBF.Mutable.isNaN , isZero @@ -42,12 +43,15 @@ , fmulInt , fmulWord , fmul2Exp + , ffma , fdiv + , frem , fsqrt , fpow , fround , frint + -- * Convert from a number , toDouble , toString @@ -68,6 +72,7 @@ import Data.Word import Data.Int import Data.Bits +import Data.Hashable import Data.List(unfoldr) import Control.Monad(foldM,when) import Control.Exception(bracket) @@ -171,7 +176,7 @@ -- | Indicates if a number is positive or negative. -data Sign = Neg {-^ Negative -} | Pos {-^ Positive -} +data Sign = Neg {-^ Negative -} | Pos {-^ Positive -} deriving (Eq,Ord,Show) @@ -344,12 +349,6 @@ isZero :: BF -> IO Bool isZero = bfQuery bf_is_zero - - - - - - foreign import capi "libbf.h bf_neg" bf_neg :: Ptr BF -> IO () @@ -377,6 +376,8 @@ foreign import ccall "bf_div" bf_div :: Ptr BF -> Ptr BF -> Ptr BF -> LimbT -> FlagsT -> IO Status +foreign import ccall "bf_rem" + bf_rem :: Ptr BF -> Ptr BF -> Ptr BF -> LimbT -> FlagsT -> CInt -> IO Status foreign import ccall "bf_pow" bf_pow :: Ptr BF -> Ptr BF -> Ptr BF -> LimbT -> FlagsT -> IO Status @@ -385,7 +386,7 @@ bf_round :: Ptr BF -> LimbT -> FlagsT -> IO Status foreign import ccall "bf_rint" - bf_rint :: Ptr BF -> LimbT -> FlagsT -> IO Status + bf_rint :: Ptr BF -> CInt -> IO Status foreign import ccall "bf_sqrt" bf_sqrt :: Ptr BF -> Ptr BF -> LimbT -> FlagsT -> IO Status @@ -401,8 +402,6 @@ fun r a b prec flags - - -- | Negate the number. fneg :: BF -> IO () fneg = bf1 bf_neg @@ -426,6 +425,21 @@ fmul :: BFOpts -> BF -> BF -> BF -> IO Status fmul = bfArith bf_mul +-- | Compute the fused-multiply-add. +-- @ffma opts x y z r@ computes @r := (x*y)+z@. +ffma :: BFOpts -> BF -> BF -> BF -> BF -> IO Status +ffma (BFOpts prec f) (BF x) (BF y) (BF z) (BF r) = + withForeignPtr x \xp -> + withForeignPtr y \yp -> + withForeignPtr z \zp -> + withForeignPtr r \out -> + do s1 <- bf_mul out xp yp #{const BF_PREC_INF} #{const BF_RNDN} + case s1 of + MemError -> return s1 + _ -> + do s2 <- bf_add out out zp prec f + pure (s1 <> s2) + -- | Multiply the number by the given word, and store the result -- in the second number. fmulWord :: BFOpts -> BF -> Word64 -> BF -> IO Status @@ -445,6 +459,16 @@ fdiv :: BFOpts -> BF -> BF -> BF -> IO Status fdiv = bfArith bf_div +-- | Compute the remainder @x - y * n@ where @n@ is the integer +-- nearest to @x/y@ (with ties broken to even values of @n@). +-- Output is written into the final argument. +frem :: BFOpts -> BF -> BF -> BF -> IO Status +frem (BFOpts p f) (BF fin1) (BF fin2) (BF fout) = + withForeignPtr fin1 \in1 -> + withForeignPtr fin2 \in2 -> + withForeignPtr fout \out -> + bf_rem out in1 in2 p f #{const BF_RNDN} + -- | Compute the square root of the first number and store the result -- in the second. fsqrt :: BFOpts -> BF -> BF -> IO Status @@ -455,8 +479,8 @@ fround (BFOpts p f) = bf1 (\ptr -> bf_round ptr p f) -- | Round to the neareset integer. -frint :: BFOpts -> BF -> IO Status -frint (BFOpts p f) = bf1 (\ptr -> bf_rint ptr p f) +frint :: RoundMode -> BF -> IO Status +frint (RoundMode r) = bf1 (\ptr -> bf_rint ptr (fromIntegral r :: CInt)) -- | Exponentiate the first number by the second, -- and store the result in the third number. @@ -536,12 +560,22 @@ | BFNaN -- ^ Not a number deriving (Eq,Ord,Show) --- | Representations for unsign floating point numbers. +instance Hashable BFRep where + hashWithSalt s BFNaN = s `hashWithSalt` (0::Int) + hashWithSalt s (BFRep Pos num) = s `hashWithSalt` (1::Int) `hashWithSalt` num + hashWithSalt s (BFRep Neg num) = s `hashWithSalt` (2::Int) `hashWithSalt` num + +-- | Representations for unsigned floating point numbers. data BFNum = Zero -- ^ zero | Num Integer !Int64 -- ^ @x * 2 ^ y@ | Inf -- ^ infinity deriving (Eq,Ord,Show) +instance Hashable BFNum where + hashWithSalt s Zero = s `hashWithSalt` (0::Int) + hashWithSalt s (Num mag ex) = s `hashWithSalt` (1::Int) `hashWithSalt` mag `hashWithSalt` ex + hashWithSalt s Inf = s `hashWithSalt` (2::Int) + -- | Returns 'Nothing' for @NaN@. getSign :: BF -> IO (Maybe Sign) getSign = bf1 (\ptr -> @@ -559,8 +593,14 @@ e > #{const BF_EXP_ZERO} then Just (fromIntegral e) else Nothing) +{-| Check if the given numer is infinite. -} +isInf :: BF -> IO Bool +isInf = bf1 (\ptr -> + do e <- #{peek bf_t, expn} ptr + if | (e :: SLimbT) == #{const BF_EXP_INF} -> pure True + | otherwise -> pure False) --- | Get the represnetation of the number. +-- | Get the representation of the number. toRep :: BF -> IO BFRep toRep = bf1 (\ptr -> do s <- #{peek bf_t, sign} ptr @@ -587,4 +627,3 @@ pure (norm base bias) -- (BFRep sgn (Num base (e - bias))) ) - diff --git a/src/LibBF/Opts.hsc b/src/LibBF/Opts.hsc --- a/src/LibBF/Opts.hsc +++ b/src/LibBF/Opts.hsc @@ -14,12 +14,14 @@ -- ** Precision , precBits + , getPrecBits , precBitsMin , precBitsMax , infPrec -- ** Exponent Size , expBits + , getExpBits , expBitsMin , expBitsMax @@ -94,9 +96,13 @@ -- | Use this many bits to represent the mantissa in the computation. -- The input should be in the interval defined by 'precMin' and 'precMax' -precBits :: Int -> BFOpts +precBits :: Word -> BFOpts precBits n = BFOpts (fromIntegral n) 0 +-- | Retrieve how many bits to represent the mantissa in the computation. +getPrecBits :: BFOpts -> Word +getPrecBits (BFOpts n _) = fromIntegral n + -- | Use the given rounding mode. -- If none is specified, then the default is 'NearEven'. rnd :: RoundMode -> BFOpts @@ -119,11 +125,18 @@ foreign import capi "libbf.h bf_set_exp_bits" bf_set_exp_bits :: CInt -> FlagsT +foreign import capi "libbf.h bf_get_exp_bits" + bf_get_exp_bits :: FlagsT -> CInt + -- | Set how many bits to use to represent the exponent. -- Should fit in the range defined by 'expBitsMin' and 'expBitsMax'. expBits :: Int -> BFOpts expBits n = BFOpts 0 (bf_set_exp_bits (fromIntegral n)) +-- | Get the number of exponent bits from a @BFOpts@ value. +getExpBits :: BFOpts -> Int +getExpBits (BFOpts _ f) = fromIntegral (bf_get_exp_bits f) + {-| The smallest supported number of bits in the exponent. -} foreign import capi "libbf.h value BF_EXP_BITS_MIN" expBitsMin :: Int @@ -170,12 +183,12 @@ ShowFmt a x <> ShowFmt b y = ShowFmt (max a b) (x .|. y) {-| Show this many significant digits total . -} -showFixed :: Word64 -> ShowFmt -showFixed n = ShowFmt n #{const BF_FTOA_FORMAT_FIXED} +showFixed :: Word -> ShowFmt +showFixed n = ShowFmt (fromIntegral n) #{const BF_FTOA_FORMAT_FIXED} {-| Show this many digits after the decimal point. -} -showFrac :: Word64 -> ShowFmt -showFrac n = ShowFmt n #{const BF_FTOA_FORMAT_FRAC} +showFrac :: Word -> ShowFmt +showFrac n = ShowFmt (fromIntegral n) #{const BF_FTOA_FORMAT_FRAC} {-| Use as many digits as necessary to match the required precision rounding to nearest and the subnormal+exponent configuration of 'FlagsT'. @@ -184,20 +197,20 @@ Infinite precision, indicated by giving 'Nothing' for the precision is supported when the radix is a power of two. -} -showFree :: Maybe Word64 -> ShowFmt +showFree :: Maybe Word -> ShowFmt showFree mb = ShowFmt prec #{const BF_FTOA_FORMAT_FREE} where prec = case mb of Nothing -> #{const BF_PREC_INF} - Just n -> n + Just n -> fromIntegral n {-| same as 'showFree' but uses the minimum number of digits (takes more computation time). -} -showFreeMin :: Maybe Word64 -> ShowFmt +showFreeMin :: Maybe Word -> ShowFmt showFreeMin mb = ShowFmt prec #{const BF_FTOA_FORMAT_FREE_MIN} where prec = case mb of Nothing -> #{const BF_PREC_INF} - Just n -> n + Just n -> fromIntegral n @@ -258,6 +271,13 @@ -- | A set of flags indicating things that might go wrong. newtype Status = Status CInt deriving (Eq,Ord) + +instance Semigroup Status where + Status a <> Status b = Status (a .|. b) + +instance Monoid Status where + mempty = Ok + mappend = (<>) checkStatus :: CInt -> Status -> Bool checkStatus n (Status x) = (x .&. n) > 0