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 #-}

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

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

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

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

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

98 word <- mword