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