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