provisional fast and really sloppy unit interval generator. For basic uses should...
authorCarter Tazio Schonwald <carter.schonwald@gmail.com>
Fri, 26 May 2017 15:50:44 +0000 (11:50 -0400)
committerCarter Tazio Schonwald <carter.schonwald@gmail.com>
Fri, 26 May 2017 15:50:44 +0000 (11:50 -0400)
cmmbits/floatsAndBits.cmm [new file with mode: 0644]
random.cabal
src/Data/Distribution/FloatingInterval.hs [new file with mode: 0644]
src/Data/Distribution/Interval.hs [deleted file]
src/Data/Distribution/Normal.hs
src/Data/Random/Utils.hs [new file with mode: 0644]
testCast/WordFloat.hs [new file with mode: 0644]

diff --git a/cmmbits/floatsAndBits.cmm b/cmmbits/floatsAndBits.cmm
new file mode 100644 (file)
index 0000000..d610ca3
--- /dev/null
@@ -0,0 +1,68 @@
+#include "Cmm.h"
+#include "MachDeps.h"
+
+#if WORD_SIZE_IN_BITS == 64
+#define DOUBLE_SIZE_WDS   1
+#else
+#define DOUBLE_SIZE_WDS   2
+#endif
+
+stg_word64ToDoublezhPrivate(I64 w)
+{
+    D_ d;
+    P_ ptr;
+
+    STK_CHK_GEN_N (DOUBLE_SIZE_WDS);
+
+    reserve DOUBLE_SIZE_WDS = ptr {
+        I64[ptr] = w;
+        d = D_[ptr];
+    }
+
+    return (d);
+}
+
+stg_doubleToWord64zhPrivate(D_ d)
+{
+    I64 w;
+    P_ ptr;
+
+    STK_CHK_GEN_N (DOUBLE_SIZE_WDS);
+
+    reserve DOUBLE_SIZE_WDS = ptr {
+        D_[ptr] = d;
+        w = I64[ptr];
+    }
+
+    return (w);
+}
+
+stg_word32ToFloatzhPrivate(W_ w)
+{
+    F_ f;
+    P_ ptr;
+
+    STK_CHK_GEN_N (1);
+
+    reserve 1 = ptr {
+        I32[ptr] = %lobits32(w);
+        f = F_[ptr];
+    }
+
+    return (f);
+}
+
+stg_floatToWord32zhPrivate(F_ f)
+{
+    W_ w;
+    P_ ptr;
+
+    STK_CHK_GEN_N (1);
+
+    reserve 1 = ptr {
+        F_[ptr] = f;
+        w = TO_W_(I32[ptr]);
+    }
+
+    return (w);
+}
index 06b4dec..31944b7 100644 (file)
@@ -54,7 +54,12 @@ library
   exposed-modules:
     System.Random.SplitMix.Internal
     System.Random.PCG.Internal
+    Data.Distribution.FloatingInterval
+    Data.Distribution.Normal
+    Data.Distribution.Permutation
+    Data.Random.Utils
 
+  c-sources: cmmbits/floatsAndBits.cmm
   -- Modules included in this library but not exported.
   -- other-modules:
 
@@ -75,3 +80,17 @@ library
   -- Base language which the package is written in.
   default-language:    Haskell2010
 
+
+test-suite word_and_float
+  type: exitcode-stdio-1.0
+  main-is: WordFloat.hs
+  hs-source-dirs:
+      testCast
+  ghc-options: -Wall
+  build-depends:
+      base == 4.*
+    , random
+    , tasty == 0.11.*
+    , tasty-hunit == 0.9.*
+  other-modules:
+  default-language: Haskell2010
diff --git a/src/Data/Distribution/FloatingInterval.hs b/src/Data/Distribution/FloatingInterval.hs
new file mode 100644 (file)
index 0000000..194e60b
--- /dev/null
@@ -0,0 +1,96 @@
+
+{- | This module provides both complete and fast (sloppy)
+random samples for the unit interval [+0,1)
+
+
+complete in this context means all representable normalized floats are reachable
+
+fast means
+-}
+
+{-# LANGUAGE ScopedTypeVariables #-}
+module Data.Distribution.FloatingInterval where
+
+import Data.Bits
+import Data.Word
+import Data.Random.Utils
+
+{-
+for more info, read the IEEE 2008 Floating point spec or wikipedia.
+single precision float is also called binary32,
+double precision is also called binary64
+
+the greatest negative exponent for double precision normalized floats is -1023, and
+    53 bits after the implicit MSB of the significand
+for single precision its -126,
+  and
+
+in both cases, for normalized form (exponent != 0, theres an implicit leading 1)
+
+
+likewise, the exponent component of a floating point number is "biased"
+in the sense of being an unsigned integer, with an additive positive correction factor
+to map the actual exponent to its representation. To map from representation to value,
+subtract the bias constant
+
+binary32 bias == 127
+binary64 bias == 1023
+-}
+
+{- |  sampleUnitIntervalDoubleM uniformly samples over the [+0,1) interval  of
+  representable floating point numbers
+
+  TODO: so for now the
+-}
+sampleUnitIntervalDoubleM :: forall m . Monad m => (m Word64) -> m Double
+sampleUnitIntervalDoubleM mword = error "finish this" $ mword  {-computeMantissa
+  where
+    computeMantissa :: m Double
+    computeMantissa = do
+        wd <- mword
+        leading <- return $ countLeadingZeros wd
+        if leading == 64
+          then computeMoreMantissa 64
+        else
+          computeNormalizedSignificandWith leading wd
+
+    computeNormalizedSignificandWith:: Int -> Word64 -> m Double
+    computeNormalizedSignificandWith leadingZeros rawWord =
+        error "finish me"  mkUnitDouble  leadingZeros rawWord
+    computeMoreMantissa :: Int  -> m Double
+    computeMoreMantissa = error "finish this too"
+    --- mkDouble takes the positive version of the negative exponent
+    --- and the normalized significand (which ellides the leading 1 digit)
+    mkUnitDouble :: Word64 -> Word64 -> Double
+    mkUnitDouble negUnBiasedExponent normalSignificand = toIEEE $ undefined (negUnBiasedExponent )
+-}
+
+{- | sampleUnitIntervalDoubleReallySloppyM, using the same algorithm as in
+http://xoroshiro.di.unimi.it/#remarks, which is also used by the rand package
+in rust. It has issues, but its super fast. Generates all the representable floats
+the correspond to dyadic (binary) rationals of the form k/2^{−53}. Note that
+the lowest order bit will 0. Which is why the lowest order bit  of the random word
+is then xor'd against the corrected unit interval number in this specific
+
+
+extracted docs from the original site:
+   #include <stdint.h>
+
+    (x >> 11) * (1. / (UINT64_C(1) << 53))
+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.
+
+An alternative, faster multiplication-free operation is
+
+    #include <stdint.h>
+
+    static inline double to_double(uint64_t x) {
+       const union { uint64_t i; double d; } u = { .i = UINT64_C(0x3FF) << 52 | x >> 12 };
+       return u.d - 1.0;
+    }
+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.
+ -}
+sampleUnitIntervalDoubleReallySloppyM ::  forall m . Monad m => (m Word64) -> m Double
+sampleUnitIntervalDoubleReallySloppyM mword = do
+        word <- mword
+        return $ toIEEE $ (\x -> (  x `xor` (1 .&. word) )) $ fromIEEE $
+          toIEEE (0x3FF `unsafeShiftL` 52 .|. (word `unsafeShiftR` 12)) - 1
diff --git a/src/Data/Distribution/Interval.hs b/src/Data/Distribution/Interval.hs
deleted file mode 100644 (file)
index f164b6e..0000000
+++ /dev/null
@@ -1,4 +0,0 @@
-module Data.Distribution.Interval where
-
-
-unit
index a7f9406..de28105 100644 (file)
@@ -1,2 +1,5 @@
 module Data.Distribution.Normal where
 
+{-
+
+-}
diff --git a/src/Data/Random/Utils.hs b/src/Data/Random/Utils.hs
new file mode 100644 (file)
index 0000000..efb64d7
--- /dev/null
@@ -0,0 +1,132 @@
+{-# LANGUAGE Trustworthy #-}
+{-# LANGUAGE CPP
+           , GHCForeignImportPrim
+           , MagicHash
+           , UnboxedTuples
+           , UnliftedFFITypes
+  #-}
+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}
+
+#include "MachDeps.h"
+module Data.Random.Utils(
+    castWord64ToDouble
+    ,castDoubleToWord64
+    ,castWord32ToFloat
+    ,castFloatToWord32
+    ,CastIEEE(..)) where
+import GHC.Word(Word32(..),Word64(..))
+import GHC.Prim (Word#,Float#,Double#)
+import GHC.Types
+{-
+from commit
+aa206346e6f12c9f88fdf051185741761ea88fbb
+of the ghc git repo, due for inclusion in ghc 8.4
+
+this should be move out of random into its own micro package for pre ghc 8.4 compat
+-}
+
+
+
+
+class CastIEEE word ieee | word -> ieee , ieee -> word where
+    toIEEE :: word -> ieee
+    fromIEEE :: ieee -> word
+
+instance CastIEEE Word32 Float where
+  {-# INLINE toIEEE #-}
+  {-# INLINE fromIEEE #-}
+  toIEEE = castWord32ToFloat
+  fromIEEE = castFloatToWord32
+instance CastIEEE Word64 Double where
+  {-# INLINE toIEEE #-}
+  {-# INLINE fromIEEE #-}
+  toIEEE = castWord64ToDouble
+  fromIEEE = castDoubleToWord64
+
+
+
+
+{-
+Note [Casting from integral to floating point types]
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+To implement something like `reinterpret_cast` from C++ to go from a
+floating-point type to an integral type one might niavely think that the
+following should work:
+      cast :: Float -> Word32
+      cast (F# f#) = W32# (unsafeCoerce# f#)
+Unfortunately that is not the case, because all the `unsafeCoerce#` does is tell
+the compiler that the types have changed. When one does the above cast and
+tries to operate on the resulting `Word32` the code generator will generate code
+that performs an integer/word operation on a floating-point register, which
+results in a compile error.
+The correct way of implementing `reinterpret_cast` to implement a primpop, but
+that requires a unique implementation for all supported archetectures. The next
+best solution is to write the value from the source register to memory and then
+read it from memory into the destination register and the best way to do that
+is using CMM.
+-}
+
+
+
+
+
+
+-- | @'castWord32ToFloat' w@ does a bit-for-bit copy from an integral value
+-- to a floating-point value.
+--
+-- @since 4.10.0.0
+
+{-# INLINE castWord32ToFloat #-}
+castWord32ToFloat :: Word32 -> Float
+castWord32ToFloat (W32# w#) = F# (stgWord32ToFloat w#)
+
+foreign import prim "stg_word32ToFloatzhPrivate"
+    stgWord32ToFloat :: Word# -> Float#
+
+
+-- | @'castFloatToWord32' f@ does a bit-for-bit copy from a floating-point value
+-- to an integral value.
+--
+-- @since 4.10.0.0
+
+{-# INLINE castFloatToWord32 #-}
+castFloatToWord32 :: Float -> Word32
+castFloatToWord32 (F# f#) = W32# (stgFloatToWord32 f#)
+
+foreign import prim "stg_floatToWord32zhPrivate"
+    stgFloatToWord32 :: Float# -> Word#
+
+
+
+-- | @'castWord64ToDouble' w@ does a bit-for-bit copy from an integral value
+-- to a floating-point value.
+--
+-- @since 4.10.0.0
+
+{-# INLINE castWord64ToDouble #-}
+castWord64ToDouble :: Word64 -> Double
+castWord64ToDouble (W64# w) = D# (stgWord64ToDouble w)
+
+foreign import prim "stg_word64ToDoublezhPrivate"
+#if WORD_SIZE_IN_BITS == 64
+    stgWord64ToDouble :: Word# -> Double#
+#else
+    stgWord64ToDouble :: Word64# -> Double#
+#endif
+
+
+-- | @'castFloatToWord32' f@ does a bit-for-bit copy from a floating-point value
+-- to an integral value.
+--
+-- @since 4.10.0.0
+
+{-# INLINE castDoubleToWord64 #-}
+castDoubleToWord64 :: Double -> Word64
+castDoubleToWord64 (D# d#) = W64# (stgDoubleToWord64 d#)
+
+foreign import prim "stg_doubleToWord64zhPrivate"
+#if WORD_SIZE_IN_BITS == 64
+    stgDoubleToWord64 :: Double# -> Word#
+#else
+    stgDoubleToWord64 :: Double# -> Word64#
+#endif
diff --git a/testCast/WordFloat.hs b/testCast/WordFloat.hs
new file mode 100644 (file)
index 0000000..b4dea21
--- /dev/null
@@ -0,0 +1,31 @@
+
+module Main where
+
+
+import Data.Random.Utils
+import Data.Word(Word32,Word64)
+import Test.Tasty
+import Test.Tasty.HUnit
+
+
+main :: IO ()
+main = defaultMain $ testGroup "zero Rep tests"
+          [testGroup "single precision"
+                    [testCase "signed zeros representation" $
+                        assertBool "impossible, binary rep of +/-0 cannot agree"
+                        (castFloatToWord32 (negate 0)  /= castFloatToWord32 0)
+
+                    ,testCase  "zero word32 is +zero" $
+                          assertBool "word32 zero is +zero, but this failure says nope"
+                            (castWord32ToFloat (0 :: Word32) == (0 :: Float ))]
+          ,testGroup "double  precision zero test "
+                    [testCase "signed zeros" $
+                      assertBool "zeros agree"
+                      (castDoubleToWord64 (negate 0)  /= castDoubleToWord64 0)
+
+                    ,testCase "zero word64 is +zero" $
+                      assertBool "word64 zero is +zero"
+                        (castWord64ToDouble (0 :: Word64) == (0 :: Double ))
+                        ]]
+
+