libBF 0.5.1 → 0.6

raw patch · 5 files changed

+301/−47 lines, 5 filesdep +hashablenew-uploaderPVP ok

version bump matches the API change (PVP)

Dependencies added: hashable

API changes (from Hackage documentation)

Files

• CHANGELOG.md +20/−0
• libBF.cabal +23/−13
• src/LibBF.hs +175/−10
• src/LibBF/Mutable.hsc +54/−15
• src/LibBF/Opts.hsc +29/−9

@@ -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.

@@ -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

@@ -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

@@ -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)))   )-

@@ -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