basic polar method rejection sampler for normal distribution
[packages/random.git] / src / Data / Distribution / FloatingInterval.hs
1
2 {- | This module provides both complete and fast (sloppy)
3 random samples for the unit interval [+0,1)
4
5
6 complete in this context means all representable normalized floats are reachable
7
8 fast means
9 -}
10
11 {-# LANGUAGE ScopedTypeVariables #-}
12 module Data.Distribution.FloatingInterval where
13
14 import Data.Bits
15 import Data.Word
16 import Data.Random.Utils
17
18 {-
19 for more info, read the IEEE 2008 Floating point spec or wikipedia.
20 single precision float is also called binary32,
21 double precision is also called binary64
22
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
27
28 in both cases, for normalized form (exponent != 0, theres an implicit leading 1)
29
30
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
35
36 binary32 bias == 127
37 binary64 bias == 1023
38 -}
39
40 {- | sampleUnitIntervalDoubleM uniformly samples over the [+0,1) interval of
41 representable floating point numbers
42
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
51 leading <- return $ countLeadingZeros wd
52 if leading == 64
53 then computeMoreMantissa 64
54 else
55 computeNormalizedSignificandWith leading wd
56
57 computeNormalizedSignificandWith:: Int -> Word64 -> m Double
58 computeNormalizedSignificandWith leadingZeros rawWord =
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 -}
67
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
74
75
76 extracted docs from the original site:
77 """
78 #include <stdint.h>
79
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.
82
83 An alternative, faster multiplication-free operation is
84
85 #include <stdint.h>
86
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 """
93
94
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