Randomish.hs

{-# LANGUAGE BangPatterns #-}

module Geometry.Randomish
        ( randomishPoints
        , randomishInts
        , randomishDoubles)
where
import Data.Word
import qualified Data.Vector.Generic            as G
import qualified Data.Vector.Unboxed.Mutable    as MV
import qualified Data.Vector.Unboxed            as V

-- | Some uniformly distributed points
randomishPoints
        :: Int                  -- ^ seed
        -> Int                  -- ^ number of points
        -> Float                -- ^ minimum coordinate
        -> Float                -- ^ maximum coordinate
        -> V.Vector (Float, Float)

randomishPoints seed' n pointMin pointMax
 = let  pts      = randomishFloats (n*2) pointMin pointMax seed'
        xs       = G.slice 0 n pts
        ys       = G.slice n n pts
   in   V.zip xs ys


-- | Use the "minimal standard" Lehmer generator to quickly generate some random
--   numbers with reasonable statistical properties. By "reasonable" we mean good
--   enough for games and test data, but not cryptography or anything where the
--   quality of the randomness really matters.
--
--   From "Random Number Generators: Good ones are hard to find"
--   Stephen K. Park and Keith W. Miller.
--   Communications of the ACM, Oct 1988, Volume 31, Number 10.
--
randomishInts
        :: Int                  -- Length of vector.
        -> Int                  -- Minumum value in output.
        -> Int                  -- Maximum value in output.
        -> Int                  -- Random seed.
        -> V.Vector Int         -- Vector of random numbers.

randomishInts !len !valMin' !valMax' !seed'

 = let  -- a magic number (don't change it)
        multiplier :: Word64
        multiplier = 16807

        -- a merzenne prime (don't change it)
        modulus :: Word64
        modulus = 2^(31 :: Integer) - 1

        -- if the seed is 0 all the numbers in the sequence are also 0.
        seed
         | seed' == 0   = 1
         | otherwise    = seed'

        !valMin = fromIntegral valMin'
        !valMax = fromIntegral valMax' + 1
        !range  = valMax - valMin

        {-# INLINE f #-}
        f x             = multiplier * x `mod` modulus
 in G.create
     $ do
        vec             <- MV.new len

        let go !ix !x
                | ix == len     = return ()
                | otherwise
                = do    let x'  = f x
                        MV.write vec ix $ fromIntegral $ (x `mod` range) + valMin
                        go (ix + 1) x'

        go 0 (f $ f $ f $ fromIntegral seed)
        return vec


-- | Generate some randomish doubles with terrible statistical properties.
--   This is good enough for test data, but not much else.
randomishDoubles
        :: Int                  -- Length of vector
        -> Double               -- Minimum value in output
        -> Double               -- Maximum value in output
        -> Int                  -- Random seed.
        -> V.Vector Double      -- Vector of randomish doubles.

randomishDoubles !len !valMin !valMax !seed
 = let  range   = valMax - valMin

        mx      = 2^(30 :: Integer) - 1
        mxf     = fromIntegral mx
        ints    = randomishInts len 0 mx seed

   in   V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints


-- | Generate some randomish doubles with terrible statistical properties.
--   This is good enough for test data, but not much else.
randomishFloats
        :: Int                  -- Length of vector
        -> Float                -- Minimum value in output
        -> Float                -- Maximum value in output
        -> Int                  -- Random seed.
        -> V.Vector Float       -- Vector of randomish doubles.

randomishFloats !len !valMin !valMax !seed
 = let  range   = valMax - valMin

        mx      = 2^(30 :: Integer) - 1
        mxf     = fromIntegral mx
        ints    = randomishInts len 0 mx seed

   in   V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints