-- | 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 t -> t -> t
forall a. Num a => a -> a -> a
* t
x0 t -> t -> t
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 = (Word64 -> Word64) -> Word64 -> Word64 -> Word64 -> Word64
forall t. Num t => (t -> t) -> t -> t -> t -> t
lcg Word64 -> Word64
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 = (Int32 -> Int32) -> Int32 -> Int32 -> Int32 -> Int32
forall t. Num t => (t -> t) -> t -> t -> t -> t
lcg Int32 -> Int32
forall a. a -> a
id Int32
1103515245 Int32
12345

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

-- | Iterate lcgWord64Knuth using cpuTimeSeedWord64.
lcgWord64KnuthCpuTime :: IO [Word64]
lcgWord64KnuthCpuTime :: IO [Word64]
lcgWord64KnuthCpuTime = (Word64 -> [Word64]) -> IO Word64 -> IO [Word64]
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((Word64 -> Word64) -> Word64 -> [Word64]
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 = Word64 -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Word64 -> Int -> Word64
forall a. Bits a => a -> Int -> a
shiftR Word64
n Int
11) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Word64 -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Word64 -> Int -> Word64
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 = ([Word64] -> [Double]) -> IO [Word64] -> IO [Double]
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((Word64 -> Double) -> [Word64] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Word64 -> Double
word64ToUnitDouble) IO [Word64]
lcgWord64KnuthCpuTime