{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}

module QIO.QIORandom where

import Data.Monoid as Monoid
import QIO.QioSyn
import QIO.Qdata
import QIO.Qio
import Complex

-- complex numbers
type CC = Complex RR

rX :: RR -> Rotation
rX r (x,y) = if x==y then (cos (r/2):+0) else (0:+ (-(sin (r/2))))

rY :: RR -> Rotation
rY r (x,y) = if x==y then (cos (r/2):+0) else (s * sin (r/2):+0) where s = if x then 1 else -1

hadamards :: [Qbit] -> U
hadamards [] = mempty
hadamards (q:qs) = uhad q `mappend` hadamards qs

pow2 :: Int -> Int
pow2 x = pow2' x 0

pow2' :: Int -> Int -> Int
pow2' x y | 2^(y+1) > x = 2^y
	  | otherwise = pow2' x (y+1)

weightedU :: RR -> Qbit -> U
weightedU ps q | sqrt ps <= 1 = rot q (rX (2*(acos (sqrt ps))))
	       | otherwise = error ("weightedU: Invalid Probability: " ++ show ps) 		

weightedBool :: RR -> QIO Bool
weightedBool pf = do q <- mkQbit False
                     applyU (weightedU pf q)
                     measQ q

rlf :: [Bool] -> [Bool]
rlf (False:bs) = rlf bs
rlf bs = bs

rlf_l :: Int -> [Bool]
rlf_l x = rlf (reverse (int2bits x))

rlf_n :: Int -> Int
rlf_n x = length (rlf_l x)

trim :: Int -> [Qbit] -> [Qbit]
trim x qbs = drop ((length qbs)-(rlf_n x)) qbs

randomU :: Int -> [Qbit] -> U
randomU max qbs = randomU' max (trim max qbs)

randomU' :: Int -> [Qbit] -> U
randomU' _ [] = mempty
randomU' 0 _ = mempty
randomU' max (q:qbs) = weightedU (fromIntegral ((max+1)-p)/fromIntegral (max+1)) q
		       `mappend`
		       condQ q (\x -> if x then (randomU (max-p) qbs) 
			                  else (hadamards qbs))
		        where p = pow2 max

randomQInt :: Int -> QIO QInt
randomQInt max = do qbs <- mkQ (reverse (int2bits max))
		    applyU (randomU max qbs)
                    return (QInt (reverse qbs))

randomQIO :: (Int,Int) -> QIO Int
randomQIO (min,max) = do q <- randomInt (max-min)
                         return (q + min)

randomInt :: Int -> QIO Int
randomInt max = do q <- randomQInt max
	           measQ q

random :: Int -> QIO Int
random x = randomInt (x-1)
		    
		    
dice :: IO Int
dice = do x <- run (randomInt 5)
	  return (x+1)


inf_dice :: IO [Int]
inf_dice = do x <- dice
	      y <- inf_dice
              return (x:y)

dice_rolls :: Int -> IO [Int]
dice_rolls 0 = return []
dice_rolls y = do x <- dice
	          xs <- dice_rolls (y-1)
                  return (x:xs)


occs :: [Int] -> (Int,Int,Int,Int,Int,Int)
occs rs = (rs' rs 1,rs' rs 2,rs' rs 3,rs' rs 4,rs' rs 5,rs' rs 6)

rs' :: [Int] -> Int -> Int
rs' rs x = length ([y|y<-rs,y==x])

probs' :: Int -> IO (Int,Int,Int,Int,Int,Int)
probs' x = do xs <- dice_rolls x
	      return (occs xs)

probs :: Int -> IO (RR,RR,RR,RR,RR,RR)
probs x = do (a,b,c,d,e,f) <- probs' x
	     return (fromIntegral a/x',fromIntegral b/x',fromIntegral c/x',fromIntegral d/x',fromIntegral e/x',fromIntegral f/x')
		    where x' = fromIntegral x