{-|
Module      : TimeSeries.RPG64
Copyright   : (C) 2013 Parallel Scientific Labs, LLC
License     : GPL-2
Stability   : experimental
Portability : portable

Computing pseudorandom 64-bit words.
We use RC5 as a basis for our PRG.
-}

module TimeSeries.PRG64 where

import Data.Bits
import Data.List (intersperse, unfoldr)
import Data.Word

-- | Random unit vector and control vector, with specified size.
--
-- Same random sequence will returned with same seed, so that sketch
-- distance could be calculated with same random sequence.
--
rnd ::
  Word64     -- ^ /bw/, size of unit vector.
  -> Word64  -- ^ /nb/, size of control vector.
  -> Integer -- ^ Seed.
  -> ([Double], [Double]) -- ^ (/unit vector/, /control vector/).
rnd bw nb seed = (r,b) where
  r = take (fromIntegral bw) $ unfoldr f (prg64Init seed1)
  b = take (fromIntegral nb) $ unfoldr f (prg64Init seed2)
  f g = let (g',w64) = prg64Next g in Just (fromIntegral w64 / max64, g')
  seed1 = seed
  seed2 = seed + 0x9387429
  max64 = fromIntegral (maxBound :: Word64)


data PRG64 = PRG64 {
                prg64A, prg64B :: !Word64
        }
        deriving Show

-- |Key pairs for rounds a-la RC5. We have to have at lest 6 rounds for
-- avalanche effect.
prg64Keys :: [(Word64, Word64)]
prg64Keys =
        [ (0xb7e151628aed2a6b, 0x9e3779b97f4a7c15)      -- P and Q numbers for RC5-64 (from paper)
        -- all other values are invented.
        , (0x1234567890abcdef, 0xfedcba9876543210)      -- two dummy constants.
        , (0x0000000000000001, 0x0000000000000002)      -- also two dumm constants.
        , (0xC045000000000000, 0x34f5b9d947cc1a58)      -- 42 in floating point double, part of commit hash.
        , (0x1bbf46f5c26e05c7, 0xd44cae3c25bd8080)      -- parts of commit hashes.
        , (0xe593d49150265452, 0xdb1953f8241b620e)      -- also parts of commit hashes.
        , (0x22bf42c19c5591b8, 0xb1a3180c5032e4c2)      -- parts of commit hashes.
        ]

-- |Compute the next value.
prg64Next :: PRG64 -> (PRG64, Word64)
prg64Next (PRG64 a0 b0) = (PRG64 ae be, ae `xor` b0)
        where
                (ae,be) = foldr rounding (a0,b0) prg64Keys
                rounding (a,b) (s1,s2) = (a', b')
                        where
                                a' = rotateL (xor a b) (fromIntegral b) + s1
                                b' = rotateL (xor a' b) (fromIntegral a') + s2


prg64Init :: Integer -> PRG64
prg64Init seed = result
        where
                (initValue, nRounds') = seed `divMod` 11
                nRounds = fromIntegral (nRounds' + 3)
                a = fromIntegral $ shiftR initValue 64
                b = fromIntegral $ initValue
                result = fst $ head $ drop nRounds $ iterate (prg64Next . fst) (PRG64 a b, 0)

prg64Bits :: Num a => PRG64 -> [[a]]
prg64Bits prg = map randomBits randoms
        where
                randoms = map snd $ tail $ iterate (prg64Next . fst) (prg,0)
                bitsMasks = map (shiftL 1) [0..63]
                randomBit r mask = if (r .&. mask) == 0 then -1 else 1
                randomBits r = map (randomBit r) bitsMasks


generateRandomVectorsCSV :: Integer -> Int -> Int -> String
generateRandomVectorsCSV seed nbits nvalues
        | nbits > 64 = error "Too many bits."
        | otherwise = text
        where
                bits :: [[Int]]
                bits = prg64Bits (prg64Init seed)
                text = unlines $
                        map (concat . intersperse "\t" . map show . take nbits) $ take nvalues bits