{-# LANGUAGE BangPatterns #-}

module Data.Array.Repa.Algorithms.Randomish
        ( randomishIntArray
        , randomishIntVector
        , randomishDoubleArray
        , randomishDoubleVector)
where
import Data.Word
import Data.Vector.Unboxed                      (Vector)
import Data.Array.Repa                          as R
import qualified Data.Vector.Unboxed.Mutable    as MV
import qualified Data.Vector.Unboxed            as V
import qualified Data.Vector.Generic            as G


-- | 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. 
--
--   By nature of the algorithm, the maximum value in the output is clipped
--   to (valMin + 2^31 - 1)
-- 
--   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.
--
randomishIntArray
        :: Shape sh
        => sh                   -- ^ Shape of array
        -> Int                  -- ^ Minumum value in output.
        -> Int                  -- ^ Maximum value in output.
        -> Int                  -- ^ Random seed.       
        -> Array U sh Int       -- ^ Array of randomish numbers.

randomishIntArray :: sh -> Int -> Int -> Int -> Array U sh Int
randomishIntArray !sh
sh !Int
valMin !Int
valMax !Int
seed
        = sh -> Vector Int -> Array U sh Int
forall sh e. sh -> Vector e -> Array U sh e
fromUnboxed sh
sh (Vector Int -> Array U sh Int) -> Vector Int -> Array U sh Int
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Int -> Int -> Vector Int
randomishIntVector (sh -> Int
forall sh. Shape sh => sh -> Int
R.size sh
sh) Int
valMin Int
valMax Int
seed


randomishIntVector 
        :: Int                  -- ^ Length of vector.
        -> Int                  -- ^ Minumum value in output.
        -> Int                  -- ^ Maximum value in output.
        -> Int                  -- ^ Random seed.       
        -> Vector Int           -- ^ Vector of randomish numbers.

randomishIntVector :: Int -> Int -> Int -> Int -> Vector Int
randomishIntVector !Int
len !Int
valMin' !Int
valMax' !Int
seed'
 = let  -- a magic number
        -- (don't change it, the randomness depends on this specific number).
        multiplier :: Word64
        multiplier :: Word64
multiplier = Word64
16807

        -- a merzenne prime
        -- (don't change it, the randomness depends on this specific number).
        modulus :: Word64
        modulus :: Word64
modulus = Word64
2Word64 -> Integer -> Word64
forall a b. (Num a, Integral b) => a -> b -> a
^(Integer
31 :: Integer) Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
- Word64
1

        -- if the seed is 0 all the numbers in the sequence are also 0.
        seed :: Int
seed    
         | Int
seed' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0   = Int
1
         | Bool
otherwise    = Int
seed'

        !valMin :: Word64
valMin = Int -> Word64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
valMin'
        !valMax :: Word64
valMax = Int -> Word64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
valMax' Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
+ Word64
1
        !range :: Word64
range  = Word64
valMax Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
- Word64
valMin

        {-# INLINE f #-}
        f :: Word64 -> Word64
f Word64
x             = Word64
multiplier Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
* Word64
x Word64 -> Word64 -> Word64
forall a. Integral a => a -> a -> a
`mod` Word64
modulus
 in (forall s. ST s (Mutable Vector s Int)) -> Vector Int
forall (v :: * -> *) a.
Vector v a =>
(forall s. ST s (Mutable v s a)) -> v a
G.create 
     ((forall s. ST s (Mutable Vector s Int)) -> Vector Int)
-> (forall s. ST s (Mutable Vector s Int)) -> Vector Int
forall a b. (a -> b) -> a -> b
$ do       
        MVector s Int
vec             <- Int -> ST s (MVector (PrimState (ST s)) Int)
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
Int -> m (MVector (PrimState m) a)
MV.new Int
len

        let go :: Int -> Word64 -> ST s ()
go !Int
ix !Word64
x 
                | Int
ix Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
len     = () -> ST s ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
                | Bool
otherwise
                = do    let x' :: Word64
x'  = Word64 -> Word64
f Word64
x
                        MVector (PrimState (ST s)) Int -> Int -> Int -> ST s ()
forall (m :: * -> *) a.
(PrimMonad m, Unbox a) =>
MVector (PrimState m) a -> Int -> a -> m ()
MV.write MVector s Int
MVector (PrimState (ST s)) Int
vec Int
ix (Int -> ST s ()) -> Int -> ST s ()
forall a b. (a -> b) -> a -> b
$ Word64 -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Word64 -> Int) -> Word64 -> Int
forall a b. (a -> b) -> a -> b
$ (Word64
x Word64 -> Word64 -> Word64
forall a. Integral a => a -> a -> a
`mod` Word64
range) Word64 -> Word64 -> Word64
forall a. Num a => a -> a -> a
+ Word64
valMin
                        Int -> Word64 -> ST s ()
go (Int
ix Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Word64
x'

        Int -> Word64 -> ST s ()
go Int
0 (Word64 -> Word64
f (Word64 -> Word64) -> Word64 -> Word64
forall a b. (a -> b) -> a -> b
$ Word64 -> Word64
f (Word64 -> Word64) -> Word64 -> Word64
forall a b. (a -> b) -> a -> b
$ Word64 -> Word64
f (Word64 -> Word64) -> Word64 -> Word64
forall a b. (a -> b) -> a -> b
$ Int -> Word64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
seed)
        MVector s Int -> ST s (MVector s Int)
forall (m :: * -> *) a. Monad m => a -> m a
return MVector s Int
vec


-- | Generate some randomish doubles with terrible statistical properties.
--   This just takes randomish ints then scales them, so there's not much randomness in low-order bits.
randomishDoubleArray
        :: Shape sh
        => sh                   -- ^ Shape of array
        -> Double               -- ^ Minumum value in output.
        -> Double               -- ^ Maximum value in output.
        -> Int                  -- ^ Random seed.       
        -> Array U sh Double    -- ^ Array of randomish numbers.

randomishDoubleArray :: sh -> Double -> Double -> Int -> Array U sh Double
randomishDoubleArray !sh
sh !Double
valMin !Double
valMax !Int
seed
        = sh -> Vector Double -> Array U sh Double
forall sh e. sh -> Vector e -> Array U sh e
fromUnboxed sh
sh (Vector Double -> Array U sh Double)
-> Vector Double -> Array U sh Double
forall a b. (a -> b) -> a -> b
$ Int -> Double -> Double -> Int -> Vector Double
randomishDoubleVector (sh -> Int
forall sh. Shape sh => sh -> Int
R.size sh
sh) Double
valMin Double
valMax Int
seed


-- | Generate some randomish doubles with terrible statistical properties.
--   This just takes randmish ints then scales them, so there's not much randomness in low-order bits.
randomishDoubleVector
        :: Int                  -- ^ Length of vector
        -> Double               -- ^ Minimum value in output
        -> Double               -- ^ Maximum value in output
        -> Int                  -- ^ Random seed.
        -> Vector Double        -- ^ Vector of randomish doubles.

randomishDoubleVector :: Int -> Double -> Double -> Int -> Vector Double
randomishDoubleVector !Int
len !Double
valMin !Double
valMax !Int
seed
 = let  range :: Double
range   = Double
valMax Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
valMin

        mx :: a
mx      = a
2a -> Integer -> a
forall a b. (Num a, Integral b) => a -> b -> a
^(Integer
30 :: Integer) a -> a -> a
forall a. Num a => a -> a -> a
- a
1
        mxf :: b
mxf     = Integer -> b
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer
forall a. Num a => a
mx :: Integer)
        ints :: Vector Int
ints    = Int -> Int -> Int -> Int -> Vector Int
randomishIntVector Int
len Int
0 Int
forall a. Num a => a
mx Int
seed
        
   in   (Int -> Double) -> Vector Int -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map (\Int
n -> Double
valMin Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
forall a. Num a => a
mxf) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
range) Vector Int
ints