-- | Linear congruential generator and related utilities.  Ordinarily use System.Random.
module Sound.Sc3.Common.Random where

import Data.Bits {- base -}
import Data.Int {- base -}
import Data.Word {- base -}
import System.CPUTime {- base -}

-- | Linear congruential generator given modulo function for type.
--
-- See <http://en.wikipedia.org/wiki/Linear_congruential_generator> for possible parameters.
lcg :: Num t => (t -> t) -> t -> t -> t -> t
lcg :: forall t. Num t => (t -> t) -> t -> t -> t -> t
lcg t -> t
modFunc t
a t
c t
x0 = t -> t
modFunc (t
a forall a. Num a => a -> a -> a
* t
x0 forall a. Num a => a -> a -> a
+ t
c)

{- | 'lcg' 6364136223846793005 1442695040888963407, so in (0, 18446744073709551615)

> take 5 (iterate lcgWord64Knuth 147092873413)
> (maxBound :: Word64) == (2 ^ 64 - 1)
-}
lcgWord64Knuth :: Word64 -> Word64
lcgWord64Knuth :: Word64 -> Word64
lcgWord64Knuth = forall t. Num t => (t -> t) -> t -> t -> t -> t
lcg forall a. a -> a
id Word64
6364136223846793005 Word64
1442695040888963407

{- | 'lcg' 1103515245 12345, so in (-2147483648, 2147483647)

> take 5 (iterate lcgInt32Glibc 873413)
> (minBound :: Int32,maxBound :: Int32) == (-2147483648, 2147483647)
-}
lcgInt32Glibc :: Int32 -> Int32
lcgInt32Glibc :: Int32 -> Int32
lcgInt32Glibc = forall t. Num t => (t -> t) -> t -> t -> t -> t
lcg forall a. a -> a
id Int32
1103515245 Int32
12345

-- | Run getCPUTime and convert to Word64
cpuTimeSeedWord64 :: IO Word64
cpuTimeSeedWord64 :: IO Word64
cpuTimeSeedWord64 = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. Integral a => a -> a -> a
div Integer
cpuTimePrecision) IO Integer
getCPUTime

-- | Run getCPUTime and convert to Int32
cpuTimeSeedInt32 :: IO Int32
cpuTimeSeedInt32 :: IO Int32
cpuTimeSeedInt32 = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c. (a -> b -> c) -> b -> a -> c
flip forall a. Integral a => a -> a -> a
div Integer
cpuTimePrecision) IO Integer
getCPUTime

-- | Iterate lcgWord64Knuth using cpuTimeSeedWord64.
lcgWord64KnuthCpuTime :: IO [Word64]
lcgWord64KnuthCpuTime :: IO [Word64]
lcgWord64KnuthCpuTime = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a. (a -> a) -> a -> [a]
iterate Word64 -> Word64
lcgWord64Knuth) IO Word64
cpuTimeSeedWord64

{- | Convert Word64 to Double in range (0, 1).
     Shifts input right 11 places (ie. discards 11 least significant bits) then divide by 2^53.
-}
word64ToUnitDouble :: Word64 -> Double
word64ToUnitDouble :: Word64 -> Double
word64ToUnitDouble Word64
n = forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a. Bits a => a -> Int -> a
shiftR Word64
n Int
11) forall a. Fractional a => a -> a -> a
/ forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a. Bits a => a -> Int -> a
shiftL (Word64
1 :: Word64) Int
53)

{- | word64ToUnitDouble of lcgWord64KnuthCpuTime

> x <- fmap (take 256) lcgDoubleKnuthCpuTime
> Sound.Sc3.Plot.plot_p1_ln [x]
-}
lcgDoubleKnuthCpuTime :: IO [Double]
lcgDoubleKnuthCpuTime :: IO [Double]
lcgDoubleKnuthCpuTime = forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (forall a b. (a -> b) -> [a] -> [b]
map Word64 -> Double
word64ToUnitDouble) IO [Word64]
lcgWord64KnuthCpuTime