2 {- | This module provides both complete and fast (sloppy)
3 random samples for the unit interval [+0,1)
6 complete in this context means all representable normalized floats are reachable
8 fast means
9 -}
11 {-# LANGUAGE ScopedTypeVariables #-}
12 module Data.Distribution.FloatingInterval where
14 import Data.Bits
15 import Data.Word
16 import Data.Random.Utils
18 {-
20 single precision float is also called binary32,
21 double precision is also called binary64
23 the greatest negative exponent for double precision normalized floats is -1023, and
24 53 bits after the implicit MSB of the significand
25 for single precision its -126,
26 and
28 in both cases, for normalized form (exponent != 0, theres an implicit leading 1)
31 likewise, the exponent component of a floating point number is "biased"
32 in the sense of being an unsigned integer, with an additive positive correction factor
33 to map the actual exponent to its representation. To map from representation to value,
34 subtract the bias constant
36 binary32 bias == 127
37 binary64 bias == 1023
38 -}
40 {- | sampleUnitIntervalDoubleM uniformly samples over the [+0,1) interval of
41 representable floating point numbers
43 TODO: so for now the
44 -}
45 sampleUnitIntervalDoubleM :: forall m . Monad m => (m Word64) -> m Double
46 sampleUnitIntervalDoubleM mword = error "finish this" \$ mword {-computeMantissa
47 where
48 computeMantissa :: m Double
49 computeMantissa = do
50 wd <- mword
53 then computeMoreMantissa 64
54 else
57 computeNormalizedSignificandWith:: Int -> Word64 -> m Double
59 error "finish me" mkUnitDouble leadingZeros rawWord
60 computeMoreMantissa :: Int -> m Double
61 computeMoreMantissa = error "finish this too"
62 --- mkDouble takes the positive version of the negative exponent
63 --- and the normalized significand (which ellides the leading 1 digit)
64 mkUnitDouble :: Word64 -> Word64 -> Double
65 mkUnitDouble negUnBiasedExponent normalSignificand = toIEEE \$ undefined (negUnBiasedExponent )
66 -}
68 {- | sampleUnitIntervalDoubleReallySloppyM, using the same algorithm as in
69 http://xoroshiro.di.unimi.it/#remarks, which is also used by the rand package
70 in rust. It has issues, but its super fast. Generates all the representable floats
71 the correspond to dyadic (binary) rationals of the form k/2^{−53}. Note that
72 the lowest order bit will 0. Which is why the lowest order bit of the random word
73 is then xor'd against the corrected unit interval number in this specific
76 extracted docs from the original site:
77 """
78 #include <stdint.h>
80 (x >> 11) * (1. / (UINT64_C(1) << 53))
81 This conversion guarantees that all dyadic rationals of the form k / 2−53 will be equally likely. Note that this conversion prefers the high bits of x, but you can alternatively use the lowest bits.
83 An alternative, faster multiplication-free operation is
85 #include <stdint.h>
87 static inline double to_double(uint64_t x) {
88 const union { uint64_t i; double d; } u = { .i = UINT64_C(0x3FF) << 52 | x >> 12 };
89 return u.d - 1.0;
90 }
91 The code above cooks up by bit manipulation a real number in the interval [1..2), and then subtracts one to obtain a real number in the interval [0..1). If x is chosen uniformly among 64-bit integers, d is chosen uniformly among dyadic rationals of the form k / 2−52.
92 """
95 -}
96 sampleUnitIntervalDoubleReallySloppyM :: forall m . Monad m => (m Word64) -> m Double
97 sampleUnitIntervalDoubleReallySloppyM mword = do
98 word <- mword
99 return \$ toIEEE \$ (\x -> ( x `xor` (1 .&. word) )) \$ fromIEEE \$
100 toIEEE (0x3FF `unsafeShiftL` 52 .|. (word `unsafeShiftR` 12)) - 1