diff --git a/ChangeLog.md b/ChangeLog.md
new file mode 100644
--- /dev/null
+++ b/ChangeLog.md
@@ -0,0 +1,30 @@
+# Revision history for quil
+
+## 0.1.7.5 -- 2017-12-19
+
+* added example usage
+
+## 0.1.7.3 -- 2017-12-18
+
+* added measurement
+
+## 0.1.5.0 -- 2017-12-17
+
+* tested two- and three-qubit gates
+
+## 0.1.4.0 -- 2017-12-16
+
+* renamed product operators
+* tested one-qubit gates
+
+## 0.1.3.0 -- 2017-12-16
+
+* added Quil gates
+
+## 0.1.2.0 -- 2017-12-16
+
+* added qubit types and operations
+
+## 0.1.0.0 -- 2017-12-15
+
+* skeleton
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,20 @@
+Copyright (c) 2017 Brian W Bush
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be included
+in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
diff --git a/ReadMe.md b/ReadMe.md
new file mode 100644
--- /dev/null
+++ b/ReadMe.md
@@ -0,0 +1,38 @@
+quil: An instruction set for quantum computing
+==============================================
+
+This Haskell library implements the Quil language for quantum computing, as specified in "A Practical Quantum Instruction Set Architecture" <>, using the indexing conventions in <>.
+
+Please report issues at <>.
+
+
+Example
+-------
+
+This example creates a wavefunction for the [Bell state](https://en.wikipedia.org/wiki/Bell_state) and then performs a measurement on its highest qubit.
+
+ > import Control.Monad.Random (evalRandIO)
+ > import Data.Qubit ((^^*), groundState)
+ > import Data.Qubit.Gate (H, CNOT)
+
+ > -- Construct the Bell state.
+ > let bell = [h 0, cnot 0 1] ^^* groundState 2
+ > bell
+ 0.7071067811865475|00> + 0.7071067811865475|11> @ [1,0]
+
+ > -- Measure the Bell wavefunction.
+ > bell' <- evalRandIO $ measure [1] bell
+ > bell'
+ ([(1,0)],1.0|00> @ [1,0])
+
+ > -- Measure it again.
+ > evalRandIO $ measure [1] bell'
+ ([(1,0)],1.0|00> @ [1,0])
+
+ > -- Measure another Bell wavefunction.
+ > evalRandIO $ measure [1] bell
+ ([(1,0)],1.0|00> @ [1,0])
+
+ > -- Measure another Bell wavefunction.
+ > evalRandIO $ measure [1] bell
+ ([(1,1)],1.0|11> @ [1,0])
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/haquil.cabal b/haquil.cabal
new file mode 100644
--- /dev/null
+++ b/haquil.cabal
@@ -0,0 +1,55 @@
+name : haquil
+version : 0.1.7.5
+synopsis : A Haskell implementation of the Quil instruction set for quantum computing.
+description : This Haskell library implements the Quil language for quantum computing, as specified in "A Practical Quantum Instruction Set Architecture" \.
+
+license : MIT
+license-file : LICENSE
+copyright : (c) 2017-18 Brian W Bush
+
+author : Brian W Bush
+maintainer : Brian W Bush
+
+homepage : https://bitbucket.org/functionally/haquil
+bug-reports : https://bwbush.atlassian.net/projects/HQUIL/issues/
+package-url : https://bitbucket.org/functionally/quil/downloads/haquil-0.1.7.5.tar.gz
+
+category : Language
+stability : Unstable
+
+cabal-version : >= 1.10
+build-type : Simple
+extra-source-files : ReadMe.md
+ ChangeLog.md
+
+source-repository head
+ type : git
+ location : https://bitbucket.org/functionally/haquil
+
+library
+ exposed-modules : Data.Qubit
+ Data.Qubit.Gate
+ other-modules : Data.Int.Util
+ hs-source-dirs : src
+ build-depends : base >= 4.9 && < 4.10
+ , hTensor
+ , MonadRandom
+ , vector
+ exposed : True
+ buildable : True
+ default-language : Haskell2010
+ ghc-options : -Wall
+
+test-suite haquil-test
+ main-is : Test.hs
+ hs-source-dirs : src
+ build-depends : base
+ , hTensor
+ , MonadRandom
+ , QuickCheck
+ , template-haskell
+ , vector
+ type : exitcode-stdio-1.0
+ buildable : True
+ default-language : Haskell2010
+ ghc-options : -Wall
diff --git a/src/Data/Int/Util.hs b/src/Data/Int/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Int/Util.hs
@@ -0,0 +1,85 @@
+-----------------------------------------------------------------------------
+--
+-- Module : $Header$
+-- Copyright : (c) 2017 Brian W Bush
+-- License : MIT
+--
+-- Maintainer : Brian W Bush
+-- Stability : Stable
+-- Portability : Portable
+--
+-- | Integer functions.
+--
+-----------------------------------------------------------------------------
+
+
+{-# LANGUAGE Safe #-}
+
+
+module Data.Int.Util (
+-- * Functions
+ ilog2
+, isqrt
+, ilog2sqrt
+) where
+
+
+import Control.Arrow ((&&&))
+
+
+-- | Square root of an integer.
+--
+-- This is equivalent to:
+--
+-- @
+-- isqrt x == floor (sqrt $ from Integral x)
+-- @
+isqrt :: (Enum a, Num a, Ord a)
+ => a -- ^ The integer.
+ -> a -- ^ The integral square root.
+isqrt n
+ | n < fst (head table) = error "isqrt: illegal argument."
+ | otherwise = snd . head $ dropWhile ((< n) . fst) table
+ where
+ table = tabulate (^ (2 :: Int))
+
+
+-- | Base-two logarithm of an integer.
+--
+-- This is equivalent to:
+--
+-- @
+-- ilog2 x == floor (logBase 2 $ fromIntegral x)
+-- @
+ilog2 :: (Enum a, Integral a, Ord a)
+ => a -- ^ The integer.
+ -> a -- ^ The integral base-two logarithm.
+ilog2 n
+ | n < fst (head table) = error "ilog2: illegal argument."
+ | otherwise = snd . head $ dropWhile ((< n) . fst) table
+ where
+ table = tabulate (2^)
+
+
+-- | Base-two logarithm of the square root of an integer.
+--
+-- This is equivalent to:
+--
+-- @
+-- ilog2sqrt x == floor (logBase 2 . sqrt $ fromIntegral x)
+-- @
+ilog2sqrt :: (Enum a, Integral a, Num a, Ord a)
+ => a -- ^ The integer.
+ -> a -- ^ The integral base-two logarithm of the square root.
+ilog2sqrt n
+ | n < fst (head table) = error "ilog2sqrt: illegal argument."
+ | otherwise = snd . head $ dropWhile ((< n) . fst) table
+ where
+ table = tabulate (^ (2 :: Int))
+
+
+-- | Tabulate a function for later use.
+tabulate :: (Enum a, Num a)
+ => (a -> a) -- ^ The function.
+ -> [(a, a)] -- ^ Its tabulation.
+tabulate f = (f &&& id) <$> [0..]
diff --git a/src/Data/Qubit.hs b/src/Data/Qubit.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Qubit.hs
@@ -0,0 +1,449 @@
+-----------------------------------------------------------------------------
+--
+-- Module : $Header$
+-- Copyright : (c) 2017 Brian W Bush
+-- License : MIT
+--
+-- Maintainer : Brian W Bush
+-- Stability : Stable
+-- Portability : Portable
+--
+-- | Qubits and operations on them, using the indexing conventions of \<\>.
+--
+-----------------------------------------------------------------------------
+
+
+{-# LANGUAGE GeneralizedNewtypeDeriving #-}
+
+
+module Data.Qubit (
+-- * Types
+ QIndex
+, QState(..)
+, Amplitude
+, Wavefunction
+, Operator
+-- * Construction
+, qubit
+, pureQubit
+, qubits
+, groundState
+, pureState
+, qubitsOperator
+-- * Properties
+, wavefunctionOrder
+, wavefunctionIndices
+, wavefunctionAmplitudes
+, rawWavefunction
+, operatorOrder
+, operatorIndices
+, operatorAmplitudes
+, rawOperator
+-- * Operations
+, (^*^)
+, (^*)
+, (*^)
+, (^^*)
+, (*^^)
+, probabilities
+, project
+, measure
+, wavefunctionProbability
+) where
+
+
+import Control.Arrow (first, second)
+import Control.Monad (replicateM)
+import Control.Monad.Random.Class (fromList)
+import Control.Monad.Random.Lazy (Rand, RandomGen)
+import Data.Complex (Complex(..), magnitude)
+import Data.Function (on)
+import Data.Int.Util (ilog2, isqrt)
+import Data.List (elemIndex, groupBy, intersect, intercalate, sortBy)
+import Data.Maybe (catMaybes)
+import Numeric.LinearAlgebra.Array ((.*))
+import Numeric.LinearAlgebra.Array.Util (Idx(iName), asScalar, coords, dims, order, outers, renameExplicit, reorder)
+import Numeric.LinearAlgebra.Tensor (Tensor, Variant(..), listTensor, switch)
+
+import qualified Data.Vector.Storable as V (toList)
+
+
+-- | States of a quibit.
+data QState =
+ QState0
+ | QState1
+ deriving (Bounded, Enum, Eq, Ord)
+
+instance Read QState where
+ readsPrec = (fmap (first toEnum) .) . readsPrec
+
+instance Show QState where
+ show = show . fromEnum
+
+
+-- | Index for qubits in a wavefunction.
+type QIndex = Int
+
+
+-- | Dimension names in order of storage.
+names :: Tensor a -> [String]
+names = fmap iName . dims
+
+
+-- | Prefixes for 'Tensor' indices.
+indexPrefix :: Variant -> Char
+indexPrefix Contra = 'm'
+indexPrefix Co = 'n'
+
+
+-- | 'Tensor' index for a collection of qubits.
+indexLabels :: Variant -- ^ Whether contravariant or covariant.
+ -> [QIndex] -- ^ Indices of the qubits.
+ -> [(String, String)] -- ^ Pairs of default `Tensor` indices and qubit indices.
+indexLabels variant indices =
+ let
+ n = length indices
+ o =
+ case variant of
+ Contra -> 0
+ Co -> n
+ in
+ zipWith (\k v -> (show k, indexPrefix variant : show v)) [(1+o)..] $ reverse indices
+
+
+-- | 'Tensor' index for a wavefunction.
+wavefunctionLabels :: Int -- ^ The number of qubits.
+ -> [(String, String)] -- ^ Pairs of default 'Tensor' indices and qubit indices.
+wavefunctionLabels n = indexLabels Contra [0..(n-1)]
+
+
+-- | 'Tensor' index for an operator on qubits.
+operatorLabels :: [QIndex] -- ^ Indices of the qubits.
+ -> [(String, String)] -- ^ Pairs of default `Tensor` indices and qubit indices.
+operatorLabels indices = indexLabels Contra indices ++ indexLabels Co indices
+
+
+-- | Amplitude of a state in a wavefunction.
+type Amplitude = Complex Double
+
+
+-- | Compactly show an amplitude.
+showAmplitude :: Amplitude -> String
+showAmplitude (a :+ b)
+ | b == 0 = show a
+ | a == 0 = show b ++ "i"
+ | otherwise = "(" ++ show a ++ "+" ++ show b ++ "i)"
+
+
+-- | Show a state and its amplitude.
+type ShowState = Int -- ^ Number of qubits in the state.
+ -> ([QState], Amplitude) -- ^ The state and its amplitude.
+ -> String -- ^ The string representation.
+
+
+-- | Help to compactly show a tensor.
+showTensor :: (a -> Tensor Amplitude) -- ^ Function for extracting the tensor of amplitudes.
+ -> ShowState -- ^ Function for showing a state and its amplitude.
+ -> a -- ^ What to show.
+ -> String -- ^ The string respresentation.
+showTensor toTensor format x =
+ let
+ x' = canonicalOrder $ toTensor x
+ n = order x'
+ in
+ (++ (show $ tensorIndices x'))
+ . (++ " @ ")
+ . intercalate " + "
+ . fmap (format n)
+ . filter ((/= 0) . snd)
+ $ tensorAmplitudes x'
+
+
+-- | Indices in a tensor.
+tensorIndices :: Tensor Amplitude -- ^ The tensor.
+ -> [QIndex] -- ^ The qubit indices.
+tensorIndices =
+ fmap (read . tail)
+ . filter ((== indexPrefix Contra) . head)
+ . names
+
+
+-- | Transpose a tensor to canonical order.
+canonicalOrder :: Tensor Amplitude -> Tensor Amplitude
+canonicalOrder x =
+ let
+ f (v : i) (v' : i') = (v, - read i :: Int) `compare` (v', - read i')
+ f _ _ = undefined
+ in
+ reorder (sortBy f $ names x) x
+
+
+-- | Amplitudes in a tensor.
+tensorAmplitudes :: Tensor Amplitude -- ^ The tensor.
+ -> [([QState], Amplitude)] -- ^ List of qubit states and their amplitudes, where indices of states are ordered according to 'tensorIndices'.
+tensorAmplitudes x =
+ let
+ n = order x
+ in
+ zip (replicateM n [minBound..maxBound])
+ . V.toList
+ $ coords x
+
+
+-- | A wavefunction for qubits.
+newtype Wavefunction = Wavefunction {rawWavefunction' :: Tensor Amplitude}
+ deriving (Eq, Floating, Fractional, Num)
+
+instance Show Wavefunction where
+ show =
+ showTensor rawWavefunction'
+ $ \_ (k, v) ->
+ showAmplitude v ++ "|" ++ concatMap show k ++ ">"
+
+
+-- | The 'Tensor' encoding the wavefunction.
+rawWavefunction :: Wavefunction -> Tensor Amplitude
+rawWavefunction = canonicalOrder . rawWavefunction'
+
+
+-- | Number of qubits in a wavefunction.
+wavefunctionOrder :: Wavefunction -> Int
+wavefunctionOrder = order . rawWavefunction'
+
+
+-- | Qubit indices in a wavefunction.
+wavefunctionIndices :: Wavefunction -- ^ The wavefunction.
+ -> [QIndex] -- ^ List of qubit indices.
+wavefunctionIndices = tensorIndices . rawWavefunction'
+
+
+-- | Amplitudes of states in a qubit wavefunction.
+wavefunctionAmplitudes :: Wavefunction -- ^ The wavefunction.
+ -> [([QState], Amplitude)] -- ^ List of qubit states and their amplitudes, where indices of states are ordered according to 'wavefunctionIndices'.
+wavefunctionAmplitudes = tensorAmplitudes . rawWavefunction'
+
+
+-- | Construct a qubit from the amplitudes of its states.
+--
+-- The squares of the norms of the amplitudes must sum to one.
+qubit :: QIndex -- ^ The index of the qubit in the wavefunction.
+ -> (Amplitude, Amplitude) -- ^ The amplitude of the 0 and 1 states, respectively.
+ -> Wavefunction -- ^ The wavefunction for the qubit.
+qubit index (a0, a1) =
+ Wavefunction
+ . renameExplicit [("1", indexPrefix Contra : show index)]
+ $ listTensor [2] [a0, a1]
+
+
+-- | Construct a qubit with a pure state.
+pureQubit :: QIndex -- ^ Which qubit.
+ -> QState -- ^ The state of the qubit.
+ -> Wavefunction -- ^ The wavefunction.
+pureQubit index QState0 = qubit index (1, 0)
+pureQubit index QState1 = qubit index (0, 1)
+
+
+-- | Construct a wavefunction for the amplitudes of its qubit states.
+--
+-- Amplitudes ordered so that the 0 state appears before the 1 state and the lower qubit indices cycle faster than then higher qubit indices. For example, a two-qubit state has its amplitudes ordered |00>, |01>, |10>, |11>. This ordering can be generated as follows, where qubits are orderd from higher indices to lower ones:
+--
+-- >>> sequence $ replicate 3 [minBound..maxBound] :: [[QState]]
+-- [[0,0,0],[0,0,1],[0,1,0],[0,1,1],[1,0,0],[1,0,1],[1,1,0],[1,1,1]]
+--
+-- The squares of the norms of the amplitudes must sum to one.
+qubits :: [Amplitude] -- ^ The amplitudes.
+ -> Wavefunction -- ^ The wavefunction.
+qubits amplitudes =
+ let
+ n = ilog2 $ length amplitudes
+ in
+ Wavefunction
+ . renameExplicit (wavefunctionLabels n)
+ $ listTensor (replicate n 2) amplitudes
+
+
+-- | Construct the ground state where each qubit is in state 0.
+groundState :: Int -- ^ Number of qubits.
+ -> Wavefunction -- ^ The ground state wavefunction.
+groundState = pureState . flip replicate QState0
+
+
+-- | Constructa a pure state.
+pureState :: [QState] -- ^ The state of each qubit, ordered from higher index to lower index.
+ -> Wavefunction -- ^ The wavefunction.
+pureState =
+ Wavefunction
+ . outers
+ . (rawWavefunction' <$>)
+ . (uncurry pureQubit <$>)
+ . zip [0..]
+ . reverse
+
+
+-- | An operator on wavefunctions.
+newtype Operator = Operator {rawOperator' :: Tensor (Complex Double)}
+ deriving (Eq, Floating, Fractional, Num)
+
+instance Show Operator where
+ show =
+ showTensor rawOperator'
+ $ \n (k, v) ->
+ let
+ (kr, kc) = splitAt (n `div` 2) $ concatMap show k
+ in
+ showAmplitude v ++ "|" ++ kr ++ "><" ++ kc ++ "|"
+
+
+-- | The 'Tensor' encoding the operator.
+rawOperator :: Operator -> Tensor Amplitude
+rawOperator = canonicalOrder . rawOperator'
+
+
+-- | Number of qubits for an operator.
+operatorOrder :: Operator -> Int
+operatorOrder = isqrt . order . rawOperator'
+
+
+-- | Qubit indices in an operator.
+operatorIndices :: Operator -- ^ The operator.
+ -> [QIndex] -- ^ List of qubit indices.
+operatorIndices = tensorIndices . rawOperator'
+
+
+-- | Amplitudes of state transitions in a qubit operator.
+operatorAmplitudes :: Operator -- ^ The wavefunction.
+ -> [(([QState], [QState]), Amplitude)] -- ^ List of qubit state transitions and their amplitudes, in row-major order with the states ordered according to 'operatorIndices'.
+operatorAmplitudes (Operator x) =
+ let
+ x' = tensorAmplitudes x
+ n = length . fst $ head x'
+ in
+ first (splitAt $ n `div` 2) <$> x'
+
+
+-- | Construct an operator on qubit wavefunctions.
+--
+-- Amplitudes in row-major order where amplitudes are ordered so that the 0 state appears before the 1 state and the lower qubit indices cycle faster than then higher qubit indices. For example, a three-qubit operator has its amplitudes ordered \<00|00>, \<00|01>, \<00|10>, \<00|11>, \<01|00>, \<01|01>, \<01|10>, \<01|11>, \<10|00>, \<10|01>, \<10|10>, \<10|11>, \<11|00>, \<11|01>, \<11|10>, \<11|11>, where states in the bras and kets are correspond to the order of the first argument to `qubitsOperator`. This ordering can be generated as follows:
+--
+-- >>> fmap (splitAt 2) . sequence $ replicate (2 * 2) [minBound..maxBound] :: [([QState], [QState])]
+-- [([0,0],[0,0]),([0,0],[0,1]),([0,0],[1,0]),([0,0],[1,1]),([0,1],[0,0]),([0,1],[0,1]),([0,1],[1,0]),([0,1],[1,1]),([1,0],[0,0]),([1,0],[0,1]),([1,0],[1,0]),([1,0],[1,1]),([1,1],[0,0]),([1,1],[0,1]),([1,1],[1,0]),([1,1],[1,1])]
+--
+-- The operator must be unitary.
+qubitsOperator :: [QIndex] -- ^ The qubit indices for which the operator applies, in descending order according to \<\>.
+ -> [Amplitude] -- ^ The amplitudes of the operator matrix, in row-major order with the states ordered from higher indices to lower ones.
+ -> Operator -- ^ The wavefunction operator.
+qubitsOperator indices =
+ let
+ n = length indices
+ in
+ Operator
+ . renameExplicit (operatorLabels $ reverse indices)
+ . listTensor (replicate n 2 ++ replicate n (-2))
+
+
+-- | Multiply two wavefunction or operator tensors.
+mult :: Tensor Amplitude -> Tensor Amplitude -> Tensor Amplitude
+mult x y =
+ let
+ lft = filter ((== indexPrefix Co ) . head) $ names x
+ rght = filter ((== indexPrefix Contra) . head) $ names y
+ cmmn = fmap tail lft `intersect` fmap tail rght
+ x' = renameExplicit [(i, '@' : tail i) | i <- filter ((`elem` cmmn) . tail) lft ] x
+ y' = renameExplicit [(i, '@' : tail i) | i <- filter ((`elem` cmmn) . tail) rght] y
+ in
+ x' * y'
+
+
+-- | Apply two operators in sequence.
+(^*^) :: Operator -> Operator -> Operator
+Operator x ^*^ Operator y = Operator $ x `mult` y
+infixr 7 ^*^
+
+
+-- | Apply an operator to a wavefunction.
+(^*) :: Operator -> Wavefunction -> Wavefunction
+Operator x ^* Wavefunction y = Wavefunction $ x `mult` y
+infixr 6 ^*
+
+
+-- | Apply an operator to a wavefunction.
+(*^) :: Wavefunction -> Operator -> Wavefunction
+(*^) = flip (^*)
+infixl 6 *^
+
+
+-- | Apply a sequence of operators to a wavefunction.
+(^^*) :: Foldable t => t Operator -> Wavefunction -> Wavefunction
+(^^*) = flip . foldl $ flip (^*)
+infixr 6 ^^*
+
+
+-- | Apply a sequence of operators to a wavefunction.
+(*^^) :: Foldable t => Wavefunction -> t Operator -> Wavefunction
+(*^^) = flip (^^*)
+infixl 6 *^^
+
+
+-- | Probabilities of a selection of qubits.
+probabilities :: [QIndex] -- ^ Which qubits.
+ -> Wavefunction -- ^ The wavefunciton.
+ -> [([(QIndex, QState)], Double)] -- ^ The probabilities for the combinations of qubit states.
+probabilities indices x =
+ let
+ indices' = wavefunctionIndices x
+ positions = catMaybes $ (`elemIndex` indices') <$> indices
+ in
+ filter ((/= 0) . snd)
+ . fmap (\kvs -> (fst $ head kvs, sum $ snd <$> kvs))
+ . groupBy ((==) `on` fst)
+ $ sortBy (compare `on` fst)
+ [
+ (
+ fmap snd . filter ((`elem` positions) . fst) . zip [0..] $ zip indices' states
+ , magnitude amplitude ^ (2::Int)
+ )
+ |
+ (states, amplitude) <- wavefunctionAmplitudes x
+ ]
+
+
+-- | Project a wavefunction onto a particular state.
+project :: [(QIndex, QState)] -- ^ The qubits for the state.
+ -> Wavefunction -- ^ The wavefunction.
+ -> Wavefunction -- ^ The projected wavefunction.
+project states x =
+ let
+ y = canonicalOrder $ rawWavefunction' x
+ p =
+ canonicalOrder
+ . outers
+ $ fmap rawWavefunction'
+ [
+ case i `lookup` states of
+ Nothing -> qubit i (1, 1)
+ Just s -> pureQubit i s
+ |
+ i <- [0..(order y - 1)]
+ ]
+ z = p .* y
+ in
+ Wavefunction $ z / sqrt (z * switch z)
+
+
+-- | Measure qubits in a wavefunction.
+measure :: RandomGen g
+ => [QIndex] -- ^ Which qubits to measure.
+ -> Wavefunction -- ^ The wavefunction.
+ -> Rand g ([(QIndex, QState)], Wavefunction) -- ^ Action for the resulting measurement and wavefunction.
+measure indices x =
+ do
+ let
+ candidates = probabilities indices x
+ collapse <- fromList $ second toRational <$> candidates
+ return (collapse, project collapse x)
+
+
+-- | The total probability for the wave function, which should be 1.
+wavefunctionProbability :: Wavefunction
+ -> Amplitude
+wavefunctionProbability (Wavefunction x) = asScalar $ x * switch x
diff --git a/src/Data/Qubit/Gate.hs b/src/Data/Qubit/Gate.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Qubit/Gate.hs
@@ -0,0 +1,340 @@
+-----------------------------------------------------------------------------
+--
+-- Module : $Header$
+-- Copyright : (c) 2017 Brian W Bush
+-- License : MIT
+--
+-- Maintainer : Brian W Bush
+-- Stability : Stable
+-- Portability : Portable
+--
+-- | Gates for qubit wavefunctions, mostly as in Quil \<> and using the conventions of \<\>.
+--
+-----------------------------------------------------------------------------
+
+
+module Data.Qubit.Gate (
+-- * One-Cubit Gates
+ i
+, x
+, y
+, z
+, h
+, s
+, t
+-- * Parameterized One-Cubit Gates
+, phase
+, rx
+, ry
+, rz
+-- * Two-Cubit Gates
+, cnot
+, swap
+, iswap
+-- * Parameterized Two-Cubit Gates
+, cphase00
+, cphase01
+, cphase10
+, cphase
+, pswap
+, cz
+-- * Three-Cubit Gates
+, ccnot
+, cswap
+) where
+
+
+import Data.Complex (Complex(..), cis)
+import Data.Qubit (Operator, QIndex, qubitsOperator)
+
+
+-- | The identity gate. This corresponds to I in Quil \<>, with the indexing convention of \<\>.
+i :: QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+i i0 =
+ qubitsOperator [i0]
+ [
+ 1, 0
+ , 0, 1
+ ]
+
+
+-- | The Pauli-x gate. This corresponds to X in Quil \<>, with the indexing convention of \<\>.
+x :: QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+x i0 =
+ qubitsOperator [i0]
+ [
+ 0, 1
+ , 1, 0
+ ]
+
+
+-- | The Pauli-y gate. This corresponds to Y in Quil \<>, with the indexing convention of \<\>.
+y :: QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+y i0 =
+ qubitsOperator [i0]
+ [
+ 0 , 0 :+ (-1)
+ , 0 :+ 1, 0
+ ]
+
+
+-- | The Pauli-z gate. This corresponds to Z in Quil \<>, with the indexing convention of \<\>.
+z :: QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+z i0 =
+ qubitsOperator [i0]
+ [
+ 1, 0
+ , 0, -1
+ ]
+
+
+-- | The Hadamard gate. This corresponds to H in Quil \<>, with the indexing convention of \<\>.
+h :: QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+h i0 =
+ qubitsOperator [i0]
+ [
+ 1, 1
+ , 1, -1
+ ] / sqrt 2
+
+
+-- | The phase-shift gate. This corresponds to PHASE in Quil \<>, with the indexing convention of \<\>.
+phase :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+phase theta i0 =
+ qubitsOperator [i0]
+ [
+ 1, 0
+ , 0, cis theta
+ ]
+
+
+-- | The one-qubit 90-degree phase-shift gate. This corresponds to S in Quil \<>, with the indexing convention of \<\>.
+s :: QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+s i0 =
+ qubitsOperator [i0]
+ [
+ 1, 0
+ , 0, 0 :+ 1
+ ]
+
+
+-- | The one-qubit 45-degree phase-shift gate. This corresponds to T in Quil \<>, with the indexing convention of \<\>.
+t :: QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+t i0 =
+ qubitsOperator [i0]
+ [
+ 1, 0
+ , 0, (1 :+ 1) / sqrt 2
+ ]
+
+
+-- | A controlled phase gate on |00>. This corresponds to CPHASE00 in Quil \<>, with the indexing convention of \<\>.
+cphase00 :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+cphase00 theta i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ cis theta, 0, 0, 0
+ , 0 , 1, 0, 0
+ , 0 , 0, 1, 0
+ , 0 , 0, 0, 1
+ ]
+
+
+-- | A controlled phase gate on |01>. This corresponds to CPHASE01 in Quil \<>, with the indexing convention of \<\>.
+cphase01 :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+cphase01 theta i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ 1, 0 , 0, 0
+ , 0, cis theta, 0, 0
+ , 0, 0 , 1, 0
+ , 0, 0 , 0, 1
+ ]
+
+
+-- | A conrolled pahse gate on |10>. This corresponds to CPHASE10 in Quil \<>, with the indexing convention of \<\>.
+cphase10 :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+cphase10 theta i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ 1, 0, 0 , 0
+ , 0, 1, 0 , 0
+ , 0, 0, cis theta, 0
+ , 0, 0, 0 , 1
+ ]
+
+
+-- | A controlled phase gate on |11>. This corresponds to CPHASE in Quil \<>, with the indexing convention of \<\>.
+cphase :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+cphase theta i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ 1, 0, 0, 0
+ , 0, 1, 0, 0
+ , 0, 0, 1, 0
+ , 0, 0, 0, cis theta
+ ]
+
+
+-- | The x-rotation gate. This corresponds to RX in Quil \<>, with the indexing convention of \<\>.
+rx :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+rx theta i0 =
+ qubitsOperator [i0]
+ [
+ cos (theta / 2) :+ 0 , 0 :+ (- sin (theta / 2))
+ , 0 :+ (- sin (theta / 2)), cos (theta / 2) :+ 0
+ ]
+
+
+-- | The y-rotation gate. This corresponds to RY in Quil \<>, with the indexing convention of \<\>.
+ry :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+ry theta i0 =
+ qubitsOperator [i0]
+ [
+ cos (theta / 2) :+ 0, (- sin (theta / 2)) :+ 0
+ , sin (theta / 2) :+ 0, cos (theta / 2) :+ 0
+ ]
+
+
+-- | The z-rotation gate. This corresponds to RZ in Quil \<>, with the indexing convention of \<\>.
+rz :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+rz theta i0 =
+ qubitsOperator [i0]
+ [
+ cis (- theta / 2), 0
+ , 0 , cis (theta / 2)
+ ]
+
+
+-- | The controlled-not gate for two qubits. This corresponds to CNOT in Quil \<>, with the indexing convention of \<\>.
+cnot :: QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+cnot i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ 1, 0, 0, 0
+ , 0, 1, 0, 0
+ , 0, 0, 0, 1
+ , 0, 0, 1, 0
+ ]
+
+
+-- | The controlled-not gate by three qubits. This corresponds to CCNOT in Quil \<>, with the indexing convention of \<\>.
+ccnot :: QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the middle qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+ccnot i0 i1 i2 =
+ qubitsOperator [i0, i1, i2]
+ [
+ 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, 1, 0, 0
+ , 0, 0, 0, 0, 0, 0, 0, 1
+ , 0, 0, 0, 0, 0, 0, 1, 0
+ ]
+
+
+-- | The phase-swap gate. This corresponds to PSWAP in Quil \<>, with the indexing convention of \<\>.
+pswap :: Double -- ^ The phase angle [radians].
+ -> QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+pswap theta i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ 1, 0 , 0 , 0
+ , 0, 0 , cis theta, 0
+ , 0, cis theta, 0 , 0
+ , 0, 0 , 0 , 1
+ ]
+
+
+-- | The swap gate with 0-degree phase. This corresponds to SWAP in Quil \<>, with the indexing convention of \<\>.
+swap :: QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+swap i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ 1, 0, 0, 0
+ , 0, 0, 1, 0
+ , 0, 1, 0, 0
+ , 0, 0, 0, 1
+ ]
+
+-- | The swap gate with 90-degree phase. This corresponds to ISWAP in Quil \<>, with the indexing convention of \<\>.
+iswap :: QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+iswap i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ 1, 0 , 0 , 0
+ , 0, 0 , 0 :+ 1, 0
+ , 0, 0 :+ 1, 0 , 0
+ , 0, 0 , 0 , 1
+ ]
+
+-- | The three-qubit controlled-swap gate. This corresponds to CSWAP in Quil \<>, with the indexing convention of \<\>.
+cswap :: QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the middle qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+cswap i0 i1 i2 =
+ qubitsOperator [i0, i1, i2]
+ [
+ 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
+ ]
+
+
+-- | The controlled-z gate, with the indexing convention of \<\>.
+cz :: QIndex -- ^ The index of the higher qubit in the wavefunction.
+ -> QIndex -- ^ The index of the lower qubit in the wavefunction.
+ -> Operator -- ^ The operator.
+cz i0 i1 =
+ qubitsOperator [i0, i1]
+ [
+ 1, 0, 0, 0
+ , 0, 1, 0, 0
+ , 0, 0, 1, 0
+ , 0, 0, 0, -1
+ ]
diff --git a/src/Test.hs b/src/Test.hs
new file mode 100644
--- /dev/null
+++ b/src/Test.hs
@@ -0,0 +1,337 @@
+-----------------------------------------------------------------------------
+--
+-- Module : $Header$
+-- Copyright : (c) 2017 Brian W Bush
+-- License : MIT
+--
+-- Maintainer : Brian W Bush
+-- Stability : Stable
+-- Portability : Portable
+--
+-- | Testing qbits and operations on them.
+--
+-- The test data is derived from pyQuil \<\ and recorded in \<\>.
+-----------------------------------------------------------------------------
+
+
+{-# LANGUAGE TemplateHaskell #-}
+
+{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
+
+
+module Main (
+ main
+) where
+
+
+import Control.Monad (replicateM)
+import Data.Complex (Complex(..), cis, magnitude)
+import Data.Int.Util (ilog2)
+import Data.Qubit (QState(..), (^*), groundState, pureQubit, pureState, qubit, qubits, rawWavefunction, wavefunctionAmplitudes, wavefunctionIndices, wavefunctionOrder)
+import Numeric.LinearAlgebra.Array ((.*))
+import Numeric.LinearAlgebra.Array.Util (coords, scalar)
+import Test.QuickCheck.All (quickCheckAll)
+import System.Exit (exitFailure, exitSuccess)
+
+import qualified Data.Qubit.Gate as G
+import qualified Data.Vector.Storable as V (toList)
+
+
+-- Helper functions.
+
+i = scalar $ 0 :+ 1
+
+cos' = scalar . (:+ 0) . cos
+sin' = scalar . (:+ 0) . sin
+cis' = scalar . cis
+
+applyGate = (rawWavefunction .) . (. (qubits . V.toList . coords)) . (^*)
+
+epsilon = 1e-5
+
+xs ~~ ys = V.toList (coords xs) ~~~ V.toList (coords ys)
+infix 4 ~~
+
+(~~~) = (and .) . zipWith (\x y -> magnitude (x - y) <= epsilon)
+
+
+-- Wavefunction tests.
+
+prop_qubit n a0 a1 =
+ let
+ q = qubit n (a0, a1)
+ in
+ wavefunctionOrder q == 1
+ && wavefunctionIndices q == [n]
+ && wavefunctionAmplitudes q == [([QState0], a0), ([QState1], a1)]
+
+prop_pure_qubit n b =
+ let
+ z = if b then QState1 else QState0
+ q = pureQubit n z
+ in
+ wavefunctionOrder q == 1
+ && wavefunctionIndices q == [n]
+ && fmap fst (wavefunctionAmplitudes q) == [[QState0], [QState1]]
+ && fmap snd (wavefunctionAmplitudes q) ~~~ (if b then [0, 1] else [1, 0])
+
+prop_qubits as =
+ let
+ n = ilog2 $ length as
+ as' = take (2^n) as
+ q = qubits as'
+ in
+ length as < 2
+ || n > 15
+ || wavefunctionOrder q == n
+ && wavefunctionIndices q == reverse [0..(n-1)]
+ && fmap fst (wavefunctionAmplitudes q) == replicateM n [minBound..maxBound]
+ && fmap snd (wavefunctionAmplitudes q) ~~~ as'
+
+prop_ground n =
+ let
+ q = groundState n
+ in
+ n < 2
+ || n > 15
+ || wavefunctionOrder q == n
+ && wavefunctionIndices q == [0..(n-1)]
+ && fmap fst (wavefunctionAmplitudes q) == replicateM n [minBound..maxBound]
+ && fmap snd (wavefunctionAmplitudes q) ~~~ (1 : replicate (2^n - 1) 0)
+
+prop_pure_state bs =
+ let
+ n = length bs
+ z = [if b then QState1 else QState0 | b <- bs]
+ q = pureState z
+ in
+ n < 2
+ || n > 15
+ || wavefunctionOrder q == n
+ && wavefunctionIndices q == [0..(n-1)]
+ && fmap fst (wavefunctionAmplitudes q) == replicateM n [minBound..maxBound]
+ && fmap snd (wavefunctionAmplitudes q) ~~~ [if z' == reverse z then 1 else 0 | z' <- replicateM n [minBound..maxBound]]
+
+
+-- Tests for one-qubit gates.
+
+b0 = rawWavefunction $ pureState [QState0]
+b1 = rawWavefunction $ pureState [QState1]
+
+prop_gate_i =
+ let
+ g = applyGate $ G.i 0
+ in
+ g b0 ~~ b0
+ && g b1 ~~ b1
+
+prop_gate_x =
+ let
+ g = applyGate $ G.x 0
+ in
+ g b0 ~~ b1
+ && g b1 ~~ b0
+
+prop_gate_y =
+ let
+ g = applyGate $ G.y 0
+ in
+ g b0 ~~ i .* b1
+ && g b1 ~~ - i .* b0
+
+prop_gate_z =
+ let
+ g = applyGate $ G.z 0
+ in
+ g b0 ~~ b0
+ && g b1 ~~ - b1
+
+prop_gate_h =
+ let
+ g = applyGate $ G.h 0
+ in
+ g b0 ~~ (b0 + b1) / sqrt 2
+ && g b1 ~~ (b0 - b1) / sqrt 2
+
+prop_gate_phase theta =
+ let
+ g = applyGate $ G.phase theta 0
+ in
+ g b0 ~~ b0
+ && g b1 ~~ cis' theta .* b1
+
+prop_gate_s =
+ let
+ g = applyGate $ G.s 0
+ in
+ g b0 ~~ b0
+ && g b1 ~~ i .* b1
+
+prop_gate_t =
+ let
+ g = applyGate $ G.t 0
+ in
+ g b0 ~~ b0
+ && g b1 ~~ scalar (1 :+ 1) .* b1 / sqrt 2
+
+prop_gate_rx theta =
+ let
+ g = applyGate $ G.rx theta 0
+ theta' = theta / 2
+ in
+ g b0 ~~ cos' theta' .* b0 - sin' theta' .* i .* b1
+ && g b1 ~~ - sin' theta' .* i .* b0 + cos' theta' .* b1
+
+prop_gate_ry theta =
+ let
+ g = applyGate $ G.ry theta 0
+ theta' = theta / 2
+ in
+ g b0 ~~ cos' theta' .* b0 + sin' theta' .* b1
+ && g b1 ~~ - sin' theta' .* b0 + cos' theta' .* b1
+
+prop_gate_rz theta =
+ let
+ g = applyGate $ G.rz theta 0
+ theta' = theta / 2
+ in
+ g b0 ~~ cis' (- theta') .* b0
+ && g b1 ~~ cis' theta' .* b1
+
+
+-- Tests for two-cubit gates.
+
+b00 = rawWavefunction $ pureState [QState0, QState0]
+b01 = rawWavefunction $ pureState [QState0, QState1]
+b10 = rawWavefunction $ pureState [QState1, QState0]
+b11 = rawWavefunction $ pureState [QState1, QState1]
+
+prop_gate_cphase00 theta =
+ let
+ g = applyGate $ G.cphase00 theta 1 0
+ in
+ g b00 ~~ cis' theta .* b00
+ && g b01 ~~ b01
+ && g b10 ~~ b10
+ && g b11 ~~ b11
+
+prop_gate_cphase01 theta =
+ let
+ g = applyGate $ G.cphase01 theta 1 0
+ in
+ g b00 ~~ b00
+ && g b01 ~~ cis' theta .* b01
+ && g b10 ~~ b10
+ && g b11 ~~ b11
+
+prop_gate_cphase10 theta =
+ let
+ g = applyGate $ G.cphase10 theta 1 0
+ in
+ g b00 ~~ b00
+ && g b01 ~~ b01
+ && g b10 ~~ cis' theta .* b10
+ && g b11 ~~ b11
+
+prop_gate_cphase theta =
+ let
+ g = applyGate $ G.cphase theta 1 0
+ in
+ g b00 ~~ b00
+ && g b01 ~~ b01
+ && g b10 ~~ b10
+ && g b11 ~~ cis' theta .* b11
+
+prop_gate_cnot =
+ let
+ g = applyGate $ G.cnot 1 0
+ in
+ g b00 ~~ b00
+ && g b01 ~~ b01
+ && g b10 ~~ b11
+ && g b11 ~~ b10
+
+prop_gate_pswap theta =
+ let
+ g = applyGate $ G.pswap theta 1 0
+ in
+ g b00 ~~ b00
+ && g b01 ~~ cis' theta .* b10
+ && g b10 ~~ cis' theta .* b01
+ && g b11 ~~ b11
+
+prop_gate_swap =
+ let
+ g = applyGate $ G.swap 1 0
+ in
+ g b00 ~~ b00
+ && g b01 ~~ b10
+ && g b10 ~~ b01
+ && g b11 ~~ b11
+
+prop_gate_iswap =
+ let
+ g = applyGate $ G.iswap 1 0
+ in
+ g b00 ~~ b00
+ && g b01 ~~ i .* b10
+ && g b10 ~~ i .* b01
+ && g b11 ~~ b11
+
+prop_gate_cz =
+ let
+ g = applyGate $ G.cz 1 0
+ in
+ g b00 ~~ b00
+ && g b01 ~~ b01
+ && g b10 ~~ b10
+ && g b11 ~~ - b11
+
+
+-- Tests for three-cubit gates.
+
+b000 = rawWavefunction $ pureState [QState0, QState0, QState0]
+b001 = rawWavefunction $ pureState [QState0, QState0, QState1]
+b010 = rawWavefunction $ pureState [QState0, QState1, QState0]
+b011 = rawWavefunction $ pureState [QState0, QState1, QState1]
+b100 = rawWavefunction $ pureState [QState1, QState0, QState0]
+b101 = rawWavefunction $ pureState [QState1, QState0, QState1]
+b110 = rawWavefunction $ pureState [QState1, QState1, QState0]
+b111 = rawWavefunction $ pureState [QState1, QState1, QState1]
+
+prop_gate_ccnot =
+ let
+ g = applyGate $ G.ccnot 2 1 0
+ in
+ g b000 ~~ b000
+ && g b001 ~~ b001
+ && g b010 ~~ b010
+ && g b011 ~~ b011
+ && g b100 ~~ b100
+ && g b101 ~~ b101
+ && g b110 ~~ b111
+ && g b111 ~~ b110
+
+prop_gate_cswap =
+ let
+ g = applyGate $ G.cswap 2 1 0
+ in
+ g b000 ~~ b000
+ && g b001 ~~ b001
+ && g b010 ~~ b010
+ && g b011 ~~ b011
+ && g b100 ~~ b100
+ && g b101 ~~ b110
+ && g b110 ~~ b101
+ && g b111 ~~ b111
+
+
+-- Run tests.
+
+return []
+main =
+ do
+ success <- $quickCheckAll
+ if success
+ then exitSuccess
+ else exitFailure