1{-# Language BangPatterns #-} 2{-# Language BlockArguments #-} 3{-# Language Trustworthy #-} 4-- | Computation with high-precision floats. 5module LibBF 6 ( 7 -- * Constants 8 BigFloat 9 , bfPosZero, bfNegZero 10 , bfPosInf, bfNegInf 11 , bfNaN 12 13 -- * Conversions 14 , bfFromWord 15 , bfFromInt 16 , bfFromDouble 17 , bfFromInteger 18 , bfFromString 19 , bfToDouble 20 , bfToString 21 , bfToRep 22 , BFRep(..) 23 , BFNum(..) 24 , bfFromBits 25 , bfToBits 26 27 -- * Predicates 28 , bfIsFinite 29 , bfIsInf 30 , bfIsZero 31 , bfIsNaN 32 , bfIsNormal 33 , bfIsSubnormal 34 , bfCompare 35 , bfSign 36 , bfExponent 37 , bfIsPos 38 , bfIsNeg 39 , Sign(..) 40 41 -- * Arithmetic 42 , bfNeg, bfAbs 43 , bfAdd, bfSub, bfMul, bfDiv, bfRem 44 , bfFMA, bfMulWord, bfMulInt, bfMul2Exp 45 , bfSqrt 46 , bfPow 47 48 -- * Rounding 49 , bfRoundFloat, bfRoundInt 50 51 -- * Mutability 52 , bfUnsafeThaw 53 , bfUnsafeFreeze 54 55 -- * Limits 56 57 58 -- * Configuration 59 , module LibBF.Opts 60 ) where 61 62 63import Data.Bits 64import Data.Hashable 65import Data.Word 66import Data.Int 67import System.IO.Unsafe 68 69import LibBF.Mutable as M 70import LibBF.Opts 71import Control.DeepSeq 72 73-- | Arbitrary precision floating point numbers. 74newtype BigFloat = BigFloat BF 75 76instance NFData BigFloat where 77 rnf x = x `seq` () 78 79 80instance Show BigFloat where 81 show = bfToString 16 (showFreeMin Nothing <> addPrefix) 82 83-------------------------------------------------------------------------------- 84{-# NOINLINE ctxt #-} 85{-# OPTIONS_GHC -fno-cse #-} 86ctxt :: BFContext 87ctxt = unsafePerformIO newContext 88 89newBigFloat :: (BF -> IO ()) -> BigFloat 90newBigFloat f = unsafe $ 91 do bf <- new ctxt 92 f bf 93 pure (BigFloat bf) 94 95newBigFloat' :: (BF -> IO a) -> (BigFloat,a) 96newBigFloat' f = unsafe $ 97 do bf <- new ctxt 98 a <- f bf 99 pure (BigFloat bf, a) 100 101unsafe :: IO a -> a 102unsafe = unsafePerformIO 103 104-------------------------------------------------------------------------------- 105-- Constants 106 107-- | Positive zero. 108bfPosZero :: BigFloat 109bfPosZero = newBigFloat (setZero Pos) 110 111-- | Negative zero. 112bfNegZero :: BigFloat 113bfNegZero = newBigFloat (setZero Neg) 114 115-- | Positive infinity. 116bfPosInf :: BigFloat 117bfPosInf = newBigFloat (setInf Pos) 118 119-- | Negative infinity. 120bfNegInf :: BigFloat 121bfNegInf = newBigFloat (setInf Neg) 122 123-- | Not-a-number. 124bfNaN :: BigFloat 125bfNaN = newBigFloat setNaN 126 127-- | A floating point number corresponding to the given word. 128bfFromWord :: Word64 -> BigFloat 129bfFromWord = newBigFloat . setWord 130 131-- | A floating point number corresponding to the given int. 132bfFromInt :: Int64 -> BigFloat 133bfFromInt = newBigFloat . setInt 134 135-- | A floating point number corresponding to the given double. 136bfFromDouble :: Double -> BigFloat 137bfFromDouble = newBigFloat . setDouble 138 139-- | A floating point number corresponding to the given integer. 140bfFromInteger :: Integer -> BigFloat 141bfFromInteger = newBigFloat . setInteger 142 143-- | IEEE 754 equality 144instance Eq BigFloat where 145 BigFloat x == BigFloat y = unsafe (cmpEq x y) 146 147-- | IEEE 754 comparisons 148instance Ord BigFloat where 149 BigFloat x < BigFloat y = unsafe (cmpLT x y) 150 BigFloat x <= BigFloat y = unsafe (cmpLEQ x y) 151 152 153instance Hashable BigFloat where 154 hashWithSalt s x = hashWithSalt s (bfToRep x) 155 156 157{-| Compare the two numbers. The special values are ordered like this: 158 159 * -0 < 0 160 * NaN == NaN 161 * NaN is larger than all other numbers 162 163Note that this differs from `(<=)` 164-} 165bfCompare :: BigFloat -> BigFloat -> Ordering 166bfCompare (BigFloat x) (BigFloat y) = unsafe (cmp x y) 167 168 169-- | Is this a finite (i.e., non-infinite, non NaN) number. 170bfIsFinite :: BigFloat -> Bool 171bfIsFinite (BigFloat x) = unsafe (isFinite x) 172 173-- | Is this value NaN. 174bfIsNaN :: BigFloat -> Bool 175bfIsNaN (BigFloat x) = unsafe (M.isNaN x) 176 177-- | Is this value infinite 178bfIsInf :: BigFloat -> Bool 179bfIsInf (BigFloat x) = unsafe (isInf x) 180 181-- | This is a "normal" number, which means it is not 182-- a NaN, not a zero, not infinite, and not subnormal. 183bfIsNormal :: BFOpts -> BigFloat -> Bool 184bfIsNormal opts bf = 185 case bfToRep bf of 186 rep@(BFRep _sgn (Num _ _)) -> not (repIsSubnormal opts rep) 187 _ -> False 188 189-- | This number is "subnormal", which means it is among the smallest 190-- representable numbers for the given precision and exponent bits. 191-- These numbers differ from "normal" numbers in that they do not use 192-- an implicit leading 1 bit in the binary representation. 193bfIsSubnormal :: BFOpts -> BigFloat -> Bool 194bfIsSubnormal opts bf = repIsSubnormal opts (bfToRep bf) 195 196-- | Get the sign of a number. Returns 'Nothing' if the number is `NaN`. 197bfSign :: BigFloat -> Maybe Sign 198bfSign (BigFloat x) = unsafe (getSign x) 199 200-- | Compute the absolute value of a number. 201bfAbs :: BigFloat -> BigFloat 202bfAbs bf = 203 case bfSign bf of 204 Just Neg -> bfNeg bf 205 _ -> bf 206 207-- | Is this value positive 208bfIsPos :: BigFloat -> Bool 209bfIsPos bf = 210 case bfSign bf of 211 Just Pos -> True 212 _ -> False 213 214-- | Is this value negative 215bfIsNeg :: BigFloat -> Bool 216bfIsNeg bf = 217 case bfSign bf of 218 Just Neg -> True 219 _ -> False 220 221-- | Get the exponent for the given number. 222-- Infinity, zero and NaN do not have an exponent. 223bfExponent :: BigFloat -> Maybe Int64 224bfExponent (BigFloat x) = unsafe (getExp x) 225 226-- | Is this value a zero. 227bfIsZero :: BigFloat -> Bool 228bfIsZero (BigFloat x) = unsafe (isZero x) 229 230-- | Negate a floating point number. 231bfNeg :: BigFloat -> BigFloat 232bfNeg (BigFloat x) = newBigFloat (\bf -> setBF x bf >> fneg bf) 233 234-- | Add two numbers useing the given options. 235bfAdd :: BFOpts -> BigFloat -> BigFloat -> (BigFloat,Status) 236bfAdd opt (BigFloat x) (BigFloat y) = newBigFloat' (fadd opt x y) 237 238-- | Subtract two numbers useing the given options. 239bfSub :: BFOpts -> BigFloat -> BigFloat -> (BigFloat,Status) 240bfSub opt (BigFloat x) (BigFloat y) = newBigFloat' (fsub opt x y) 241 242-- | Multiply two numbers using the given options. 243bfMul :: BFOpts -> BigFloat -> BigFloat -> (BigFloat,Status) 244bfMul opt (BigFloat x) (BigFloat y) = newBigFloat' (fmul opt x y) 245 246-- | Multiply a number and a word, using the given options. 247bfMulWord :: BFOpts -> BigFloat -> Word64 -> (BigFloat,Status) 248bfMulWord opt (BigFloat x) y = newBigFloat' (fmulWord opt x y) 249 250-- | Multiply a number and an int, using the given options. 251bfMulInt :: BFOpts -> BigFloat -> Int64 -> (BigFloat,Status) 252bfMulInt opt (BigFloat x) y = newBigFloat' (fmulInt opt x y) 253 254-- | Multiply a number by @2^e@. 255bfMul2Exp :: BFOpts -> BigFloat -> Int64 -> (BigFloat,Status) 256bfMul2Exp opt (BigFloat x) e = newBigFloat' (\p -> 257 do setBF x p 258 fmul2Exp opt e p) 259 260-- | Divide two numbers useing the given options. 261bfDiv :: BFOpts -> BigFloat -> BigFloat -> (BigFloat,Status) 262bfDiv opt (BigFloat x) (BigFloat y) = newBigFloat' (fdiv opt x y) 263 264-- | Compute the remainder @x - y * n@ where @n@ is the integer 265-- nearest to @x/y@ (with ties broken to even values of @n@). 266bfRem :: BFOpts -> BigFloat -> BigFloat -> (BigFloat, Status) 267bfRem opt (BigFloat x) (BigFloat y) = newBigFloat' (frem opt x y) 268 269-- | Compute the fused-multiply-add @(x*y)+z@ 270bfFMA :: BFOpts -> BigFloat -> BigFloat -> BigFloat -> (BigFloat, Status) 271bfFMA opt (BigFloat x) (BigFloat y) (BigFloat z) = newBigFloat' (ffma opt x y z) 272 273-- | Square root of two numbers useing the given options. 274bfSqrt :: BFOpts -> BigFloat -> (BigFloat,Status) 275bfSqrt opt (BigFloat x) = newBigFloat' (fsqrt opt x) 276 277-- | Round to a float matching the input parameters. 278bfRoundFloat :: BFOpts -> BigFloat -> (BigFloat,Status) 279bfRoundFloat opt (BigFloat x) = newBigFloat' (\bf -> 280 do setBF x bf 281 fround opt bf 282 ) 283 284-- | Round to an integer using the given rounding mode. 285bfRoundInt :: RoundMode -> BigFloat -> (BigFloat,Status) 286bfRoundInt r (BigFloat x) = newBigFloat' (\bf -> 287 do setBF x bf 288 frint r bf 289 ) 290 291-- | Exponentiate a word to a positive integer power. 292bfPow :: BFOpts -> BigFloat -> BigFloat -> (BigFloat, Status) 293bfPow opts (BigFloat x) (BigFloat y) = newBigFloat' (fpow opts x y) 294 295-- | Constant to a 'Double' 296bfToDouble :: RoundMode -> BigFloat -> (Double, Status) 297bfToDouble r (BigFloat x) = unsafe (toDouble r x) 298 299-- | Render as a 'String', using the given settings. 300bfToString :: Int {- ^ Base -} -> ShowFmt -> BigFloat -> String 301bfToString radix opts (BigFloat x) = 302 unsafe (toString radix opts x) 303 304-- | Parse a number from the given string. 305-- Returns @NaN` if the string does not correspond to a number we recognize. 306bfFromString :: Int {- ^ Base -} -> BFOpts -> String -> (BigFloat,Status) 307bfFromString radix opts str = 308 newBigFloat' \bf -> 309 do (status,_,usedAll) <- setString radix opts str bf 310 if usedAll 311 then pure status 312 else do setNaN bf 313 pure Ok 314 315-- | The float as an exponentiated 'Integer'. 316bfToRep :: BigFloat -> BFRep 317bfToRep (BigFloat x) = unsafe (toRep x) 318 319-- | Make a number mutable. 320-- WARNING: This does not copy the number, 321-- so it could break referential transperancy. 322bfUnsafeThaw :: BigFloat -> BF 323bfUnsafeThaw (BigFloat x) = x 324 325-- | Make a number immutable. 326-- WARNING: This does not copy the number, 327-- so it could break referential transperancy. 328bfUnsafeFreeze :: BF -> BigFloat 329bfUnsafeFreeze = BigFloat 330 331-------------------------------------------------------------------------------- 332 333-- | Make a float using "raw" bits representing the bitvector 334-- representation of a floating-point value with the 335-- exponent and precision bits given by the options. 336bfFromBits :: 337 BFOpts -> 338 Integer {- ^ Raw bits -} -> 339 BigFloat 340 341bfFromBits opts bits 342 | expoBiased == 0 && mant == 0 = -- zero 343 if isNeg then bfNegZero else bfPosZero 344 345 | expoBiased == eMask && mant == 0 = -- infinity 346 if isNeg then bfNegInf else bfPosInf 347 348 | expoBiased == eMask = bfNaN -- NaN 349 350 | expoBiased == 0 = -- Subnormal 351 case bfMul2Exp opts' (bfFromInteger mant) (expoVal + 1) of 352 (num,Ok) -> if isNeg then bfNeg num else num 353 (_,s) -> error $ unwords ["bfFromBits", "subnormal case", "Unexpected status:", show s, show bits, show mant, show expoVal, show e, show p ] 354 355 | otherwise = -- Normal 356 case bfMul2Exp opts' (bfFromInteger mantVal) expoVal of 357 (num,Ok) -> if isNeg then bfNeg num else num 358 (_,s) -> error $ unwords ["bfFromBits", "normal case", "Unexpected status:", show s, show bits, show mantVal, show expoVal, show e, show p ] 359 360 where 361 e = getExpBits opts 362 p = getPrecBits opts 363 364 opts' = opts <> allowSubnormal 365 366 p' = fromIntegral p - 1 :: Int 367 eMask = (1 `shiftL` e) - 1 :: Int64 368 pMask = (1 `shiftL` p') - 1 :: Integer 369 370 isNeg = testBit bits (e + p') 371 372 mant = pMask .&. bits :: Integer 373 mantVal = mant `setBit` p' :: Integer 374 -- accounts for the implicit 1 bit 375 376 expoBiased = eMask .&. fromInteger (bits `shiftR` p') :: Int64 377 bias = eMask `shiftR` 1 :: Int64 378 expoVal = expoBiased - bias - fromIntegral p' :: Int64 379 380 381-- | Turn a float into raw bits. 382-- @NaN@ is represented as a positive "quiet" @NaN@ 383-- (most significant bit in the significand is set, the rest of it is 0). 384bfToBits :: BFOpts -> BigFloat -> Integer 385bfToBits opts bf = res 386 where 387 res = (isNeg `shiftL` (e+p')) 388 .|. (expBiased `shiftL` p') 389 .|. (mant `shiftL` 0) 390 391 e = getExpBits opts 392 p = getPrecBits opts 393 394 p' = fromIntegral p - 1 :: Int 395 396 eMask = (1 `shiftL` e) - 1 :: Integer 397 pMask = (1 `shiftL` p') - 1 :: Integer 398 399 (isNeg, expBiased, mant) = 400 case bfToRep bf of 401 BFNaN -> (0, eMask, 1 `shiftL` (p' - 1)) 402 BFRep s num -> (sign, be, ma) 403 where 404 sign = case s of 405 Neg -> 1 406 Pos -> 0 407 408 (be,ma) = 409 case num of 410 Zero -> (0,0) 411 Num i ev 412 | ex <= 0 -> 413 (0, i `shiftL` (p'-m-1+fromInteger ex)) -- subnormal case 414 | otherwise -> 415 (ex, (i `shiftL` (p' - m)) .&. pMask) -- normal case 416 where 417 m = msb 0 i - 1 418 bias = eMask `shiftR` 1 419 ex = toInteger ev + bias + toInteger m 420 421 Inf -> (eMask,0) 422 423 msb !n j = if j == 0 then n else msb (n+1) (j `shiftR` 1) 424 425-- | test if a given big float representation is subnormal 426repIsSubnormal :: BFOpts -> BFRep -> Bool 427repIsSubnormal opts (BFRep _s (Num i ev)) = ex <= 0 428 where 429 e = getExpBits opts 430 eMask = (1 `shiftL` e) - 1 :: Integer 431 bias = eMask `shiftR` 1 432 433 m = msb (0 :: Int) i - 1 434 ex = toInteger ev + bias + toInteger m 435 436 msb !n j = if j == 0 then n else msb (n+1) (j `shiftR` 1) 437 438repIsSubnormal _opts _rep = False 439