Quantum.hs

{-# LANGUAGE NoMonomorphismRestriction, FlexibleContexts #-}
module Quantum where

import Data.Complex
import Control.Monad.Random
import Control.Lens

-- |Dot product of complex vectors
dotprod a b = sum (zipWith (*) a (map conjugate b))

-- |Multiply vector with matrix
vmmult v m = map (dotprod v) m

-- |Squared norm (length) of a vector
normSq v = dotprod v v

-- |Norm (length) of a vector
norm v = sqrt (normSq v)

-- |Scale a vector by a scalar
scale v s = map (*s) v

-- |Rescale a vector so that it's a quantum state
normalizeState v = scale v (1 / norm v)

-- |Born's rule: Given probability amplitudes corresponding
-- to a given result, the probability of getting that result is the
-- sum of the squares of the magnitudes of the amplitudes
born = sum . map ((^2) . magnitude)

-- |When measuring the bit at index bitIndex in a state of numBits
-- bits, the dimensionth index in the state vector contributes to the
-- probability of measuring 0
indexOfZero numBits bitIndex dimension =
  let l = 2^(numBits - bitIndex - 1) in
  (dimension `div` l) `mod` 2 == 0

-- |Measure the bit at index bitIndex of a quantum state v with
-- numBits bits, returning the measurement result and the residual
-- state. Measurement done in the computational basis.
measure v numBits bitIndex = do
  let zeroes = traversed . indices (indexOfZero numBits bitIndex)
      ones = traversed . indices (not . indexOfZero numBits bitIndex)
  roll <- getRandomR (0.0, 1.0 :: Double)
  let prob = born (toListOf ones v)
  let result = roll < prob
  return (result, normalizeState (set (if result then zeroes else ones) 0 v))

-- Some quantum gates and utilities therefor

cnot = [[1, 0, 0, 0],
        [0, 1, 0, 0],
        [0, 0, 0, 1],
        [0, 0, 1, 0]]

cswap = [[1, 0, 0, 0, 0, 0, 0, 0],
         [0, 1, 0, 0, 0, 0, 0, 0],
         [0, 0, 1, 0, 0, 0, 0, 0],
         [0, 0, 0, 1, 0, 0, 0, 0],
         [0, 0, 0, 0, 1, 0, 0, 0],
         [0, 0, 0, 0, 0, 0, 1, 0],
         [0, 0, 0, 0, 0, 1, 0, 0],
         [0, 0, 0, 0, 0, 0, 0, 1]]

-- |Scale a matrix by a scalar.
mscale scalar = map (map (*scalar))

-- |Hadamard transformation, without scaling
hadamardUnscaled 1 = [[1, 1], [1, -1]]
hadamardUnscaled n =
    let lower = hadamardUnscaled (n - 1) in
    zipWith (++) lower lower ++
    zipWith (++) lower (mscale (-1) lower)

-- |The (2^n)X(2^n) Hadamard gate.
hadamard n = mscale (2**(-(n/2))) $ hadamardUnscaled n

-- |Tensor product
tensor s t =
  concatMap (\rowS ->
               map (\rowT ->
                      concatMap (\fieldS ->
                                   map (\fieldT -> fieldS * fieldT)
                                   rowT)
                      rowS)
               t)
  s


-- example: Measuring a Bell state

i = 0 :+ 1

-- | (1/sqrt 2) * (|00> + |11>)
bell = normalizeState [1, 0, 0, 1]

-- | Always either (True, True) or (False, False) (unless bugs)!
measureBell = do
  (bit1, state1) <- measure bell 2 0
  (bit2, state2) <- measure state1 2 1
  return (bit1, bit2)