QIO 1.2 → 1.3

raw patch · 13 files changed

+504/−149 lines, 13 filesdep ~basedep ~containersdep ~mtlPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: base, containers, mtl, old-time, random

API changes (from Hackage documentation)

Files

• LICENSE +1/−1
• QIO.cabal +16/−14
• QIO/QArith.hs +37/−30
• QIO/QExamples.hs +22/−20
• QIO/QIORandom.hs +29/−23
• QIO/Qdata.hs +2/−2
• QIO/Qft.hs +8/−7
• QIO/Qio.hs +15/−6
• QIO/QioSyn.hs +11/−3
• QIO/QioSynAlt.hs +298/−0
• QIO/Shor.hs +34/−36
• QIO/Vec.hs +17/−2
• QIO/VecEq.hs +14/−5

@@ -1,4 +1,4 @@-Copyright (c) 2010-2012, Alexander S. Green+Copyright (c) 2010-2016, Alexander S. Green All rights reserved.  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

@@ -1,20 +1,22 @@-Name:		QIO-Version:	1.2-Cabal-Version:	>= 1.2-License:	BSD3-License-File:	LICENSE-Author:		Alexander S. Green-Maintainer: 	alexander.s.green@gmail.com-Homepage:	http://www.cs.nott.ac.uk/~asg/QIO/-Category:	Quantum-Synopsis:	The Quantum IO Monad is a library for defining quantum computations in Haskell-Description:	The Quantum IO Monad is a library for defining quantum computations in Haskell. It can be thought of as an embedded language within Haskell, and comes with functions for simulating the running of these quantum computations. The distribution contains many example computations written in QIO, including an implementation of Shor's algorithm.-Build-Type:	Simple+Name:			QIO+Version:		1.3+Cabal-Version:		>= 1.24+License:		BSD3+License-File:		LICENSE+Author:			Alexander S. Green+Maintainer: 		alexander.s.green@gmail.com+Homepage:		https://github.com/alexandersgreen/qio-haskell+Category:		Quantum+Synopsis:		The Quantum IO Monad is a library for defining quantum computations in Haskell+Description:		The Quantum IO Monad is a library for defining quantum computations in Haskell. It can be thought of as an embedded language within Haskell, and comes with functions for simulating the running of these quantum computations. The distribution contains many example computations written in QIO, including an implementation of Shor's algorithm.+Build-Type:		Simple  Library+  Default-language:+    Haskell2010   Build-Depends:		-    base >= 4 && < 5, containers, mtl, random, old-time+    base >= 4.9 && < 4.10, containers, mtl, random, old-time   Exposed-modules:     QIO.Heap, QIO.QArith, QIO.QExamples, QIO.QIORandom, QIO.Qdata, QIO.Qft,-    QIO.Qio, QIO.QioClass, QIO.QioSyn, QIO.Shor, QIO.Vec, QIO.VecEq+    QIO.Qio, QIO.QioClass, QIO.QioSyn, QIO.QioSynAlt, QIO.Shor, QIO.Vec, QIO.VecEq 

@@ -15,8 +15,9 @@ -- over the underlying qubits that make up a QInt. swapQInt :: QInt -> QInt -> U swapQInt (QInt xs) (QInt ys) = swapQInt' xs ys-			       where swapQInt' [] [] = mempty-				     swapQInt' (x:xs) (y:ys) = (swap x y) `mappend` swapQInt' xs ys+  where +   swapQInt' [] [] = mempty+   swapQInt' (x:xs) (y:ys) = (swap x y) `mappend` swapQInt' xs ys  -- | ifElseQ defines a quantum If statement, whereby depending on the state of -- the given (control) qubit, one of two unitaries are applied.@@ -94,11 +95,12 @@ -- | A small function to test applying the adder unitary, and then applying -- the reverse of the adder unitary, which should give the identity function. tBiAdder :: (Int,(Int,Bool)) -> QIO (Int,(Int,Bool))-tBiAdder xyc = do q @ (qx,(qy,qc)) <- mkQ xyc-		  applyU (adder qx qy qc)-		  applyU (urev (adder qx qy qc))-                  xyc <- measQ q-	          return xyc+tBiAdder xyc = do +  q @ (qx,(qy,qc)) <- mkQ xyc+  applyU (adder qx qy qc)+  applyU (urev (adder qx qy qc))+  xyc <- measQ q+  return xyc  -- | This unitary is for modular addition, and is done modulo some fixed -- classical modulus, given as the first Int argument.@@ -136,10 +138,12 @@ -- state 0 before the operation, then it is left in the sate a*x mod n. multMod :: Int -> Int -> QInt -> QInt -> U multMod n a (QInt x) y = multMod' n a x y 1-                         where multMod' _ _ [] _ _ = mempty-			       multMod' n a (x:xs) y p = cond x (\x -> (if x then (letU ((p*a) `mod` n) (\ qa -> (adderMod n qa y)) `mappend` (multMod' n a xs y (p*2)))-                                                                             else multMod' n a xs y (p*2)))-		               +  where +   multMod' _ _ [] _ _ = mempty+   multMod' n a (x:xs) y p = cond x (\x -> (+    if x then (letU ((p*a) `mod` n) (\ qa -> (adderMod n qa y)) `mappend` (multMod' n a xs y (p*2)))+         else multMod' n a xs y (p*2)))+                        -- | A small function for testing the modular multiplication unitary. This function -- initialises 'y' as zero, so the output is as expected. tmultMod :: Int -> Int -> Int -> QIO (Int,Int)@@ -171,19 +175,21 @@ -- | The unitary that represents modular exponentiation is constructed in terms -- of multiple \"steps\". This function defines those steps. modExpStep :: Qbit -> Int -> Int -> QInt -> Int -> U-modExpStep qc n a o p = letU 0 (\z ->                (condMultMod qc n p'                o z) -			             `mappend` (ifQ qc (swapQInt o z))-                                     `mappend` (urev (condMultMod qc n (inverseMod n p') o z)))-				  where p' | (a^(2^p)) == 0 = error "modExpStep: arguments too large"-					   | otherwise = (a^(2^p)) `mod` n+modExpStep qc n a o p = letU 0 (\z -> (condMultMod qc n p' o z) +                        `mappend` (ifQ qc (swapQInt o z)) +                        `mappend` (urev (condMultMod qc n (inverseMod n p') o z)))+  where +    p' | (a^(2^p)) == 0 = error "modExpStep: arguments too large"+       | otherwise = (a^(2^p)) `mod` n  -- | A QIO computation that forms a test of the 'modExpStep' unitary modExpStept :: Int -> Int -> Int -> Int -> QIO Int-modExpStept i n a p = do q <- mkQ True-		         one <- mkQ i-		         applyU (modExpStep q n a one p)	      -		         r <- measQ one	     -		         return r+modExpStept i n a p = do +  q <- mkQ True+  one <- mkQ i+  applyU (modExpStep q n a one p)          +  r <- measQ one      +  return r  -- | This function defines a unitary that implements modular exponentiation, as -- required in Shor's algorithm. Given classical arguments \n\ and \a\, a quantum@@ -191,21 +197,22 @@ -- leave the quantum register \o\ in the state \a\^\x\ mod \n\. modExp :: Int -> Int -> QInt -> QInt -> U modExp n a (QInt x) o = modExp' n a x o 0-                        where modExp' _ _ [] _ _ = mempty-			      modExp' n a (x:xs) o p =           modExpStep x n a o p -						      `mappend` (modExp' n a xs o (p+1))+  where +    modExp' _ _ [] _ _ = mempty+    modExp' n a (x:xs) o p = modExpStep x n a o p `mappend` (modExp' n a xs o (p+1))  -- | A QIO computation that forms a test of the modular exponentiation unitary. modExpt :: Int -> (Int,Int) -> QIO Int-modExpt n (a,x) = do qx <- mkQ x-		     one <- mkQ 1-                     applyU (modExp n a qx one)-                     r <- measQ one-		     return r+modExpt n (a,x) = do +  qx <- mkQ x+  one <- mkQ 1+  applyU (modExp n a qx one)+  r <- measQ one+  return r   -		       +                  

@@ -66,20 +66,22 @@ -- the qubit twice, and then the qubit is measured. This should correspond to -- the identity function on the given Boolean value. hadTwice :: Bool -> QIO Bool-hadTwice x = do q <- mkQ x-		applyU (uhad q `mappend` uhad q)-		b <- measQ q-		return b+hadTwice x = do +  q <- mkQ x+  applyU (uhad q `mappend` uhad q)+  b <- measQ q+  return b  -- | A different implementation of 'hadTwice' where QIO is used to apply two -- unitaries, each of which is a single Hadamard gate, as opposed to a single -- unitary, which is two Hadamard gates.  hadTwice' :: Bool -> QIO Bool-hadTwice' x = do q <- mkQ x-	 	 applyU (uhad q)-                 applyU (uhad q)-		 b <- measQ q-		 return b+hadTwice' x = do +  q <- mkQ x+  applyU (uhad q)+  applyU (uhad q)+  b <- measQ q+  return b  ---------------------------------------------- ---- Teleportation ---------------------------@@ -88,11 +90,11 @@ -- | The operations that Alice must perform in the classic quantum teleportation -- example. alice :: Qbit -> Qbit -> QIO (Bool,Bool)-alice aq eq  =  do  applyU (cond aq (\a -> if a then (unot eq)-                                              else (mempty)  )  )-                    applyU (uhad aq)-                    cd <- measQ (aq,eq)-                    return cd+alice aq eq  =  do  +  applyU (cond aq (\a -> if a then (unot eq) else (mempty)))+  applyU (uhad aq)+  cd <- measQ (aq,eq)+  return cd  -- | A definition of the Pauli-Z gate. uZZ :: Qbit -> U@@ -104,8 +106,7 @@ bobsU  (False,False)  eq  =  mempty bobsU  (False,True)  eq  =  (unot eq) bobsU  (True,False)  eq  =  (uZZ eq)-bobsU  (True,True)  eq  =	     ((unot eq) -       		    	   `mappend` (uZZ eq))+bobsU  (True,True)  eq  = ((unot eq) `mappend` (uZZ eq))  -- | The overall operations that Bob must perform in the classic quantum -- teleportation example@@ -123,10 +124,11 @@ -- | A small test function of quantum teleportation, which teleports a -- bell state, and then measures it. test_teleport :: QIO (Bool,Bool)-test_teleport = do (q1,q2) <- bell-                   tq2 <- teleportation q2-		   result <- measQ (q1,tq2)-		   return result+test_teleport = do +  (q1,q2) <- bell+  tq2 <- teleportation q2+  result <- measQ (q1,tq2)+  return result  -- | teleports a qubit in the state |1> teleport_true' :: QIO Qbit

@@ -26,13 +26,13 @@ pow2 :: Int -> Int pow2 x = pow2' 0   where pow2' y | 2^(y+1) > x = 2^y-	        | otherwise = pow2' (y+1)+            | otherwise = pow2' (y+1)  -- | A rotation that, given a qubit in state 0, leaves it in a super-position of -- 0 and 1, such that the probability of measuring as state 0 is \ps\. weightedU :: RR -> Qbit -> U weightedU ps q | sqrt ps <= 1 = rot q (rX (2*(acos (sqrt ps))))-	       | otherwise = error ("weightedU: Invalid Probability: " ++ show ps) 		+           | otherwise = error ("weightedU: Invalid Probability: " ++ show ps)       -- | A QIO computation that uses the "weightedU" unitary, to return a Bool that -- has a probablity of \pf\ of being False. weightedBool :: RR -> QIO Bool@@ -69,17 +69,18 @@   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+                 `mappend`+                 condQ q (\x -> if x then (randomU (max-p) qbs) +                                else (hadamards qbs))+                  where p = pow2 max  -- | A quantum computation that will return a quantum integer in a state that -- has equal probabilities of being measured in any of the state 0 to \max\. randomQInt :: Int -> QIO QInt-randomQInt max = do qbs <- mkQ (reverse (int2bits max))-		    applyU (randomU max qbs)-                    return (QInt (reverse qbs))+randomQInt max = do +  qbs <- mkQ (reverse (int2bits max))+  applyU (randomU max qbs)+  return (QInt (reverse qbs))  -- | A quantum computation that will return a quantum integer in a state that -- has equal probabilities of being measured in any of the state \min\ to \max\.@@ -89,25 +90,28 @@  -- | A quantum computation that measures the outcome of "randomQInt" randomInt :: Int -> QIO Int-randomInt max = do q <- randomQInt max-	           measQ q+randomInt max = do +  q <- randomQInt max+  measQ q  -- | A quantum computation that returns an integer that is equally likely to be -- any number in the range 0 to \x\-1 random :: Int -> QIO Int random x = randomInt (x-1)-		    +             -- | This function uses a Quantum computation to simulate the roll of a dice dice :: IO Int-dice = do x <- run (randomInt 5)-	  return (x+1)+dice = do +  x <- run (randomInt 5)+  return (x+1)  -- | This function simulates the given number of repitions of dice rolls dice_rolls :: Int -> IO [Int] dice_rolls 0 = return []-dice_rolls y = do x <- dice-	          xs <- dice_rolls (y-1)-                  return (x:xs)+dice_rolls y = do +  x <- dice+  xs <- dice_rolls (y-1)+  return (x:xs)  -- | Returns the number of occurences of 1 through 6 in the given list of Ints occs :: [Int] -> (Int,Int,Int,Int,Int,Int)@@ -118,14 +122,16 @@ -- | Returns the number of occurences of 1 through 6 in the given number of -- rolls of the dice. probs' :: Int -> IO (Int,Int,Int,Int,Int,Int)-probs' x = do xs <- dice_rolls x-	      return (occs xs)+probs' x = do +  xs <- dice_rolls x+  return (occs xs)  -- | Returns the percentage of occurences of 1 through 6, after the given number -- of rolls of the dice. 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+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 -		   +           

@@ -44,8 +44,8 @@         where letU' [] xs = xsu xs               letU' (a:as) xs = letU a (\ x -> letU' as (xs++[x]))     condQ qs qsu = condQ' qs []-	where condQ' [] xs = qsu xs-	      condQ' (a:as) xs = condQ a (\ x -> condQ' as (xs++[x]))+        where condQ' [] xs = qsu xs+              condQ' (a:as) xs = condQ a (\ x -> condQ' as (xs++[x]))  -- | A recursive conditional on a list of quantum data condQRec :: Qdata a qa => [qa] -> [(a -> U)] -> U

@@ -22,9 +22,9 @@ -- | The \"base\" step involved in a QFT is a series of controlled rotations. qftBase :: [Bool] -> Qbit -> U qftBase bs q =  f' bs q 2-	where f' [] q _ = uhad q-	      f' (b:bs) q x = if b then (rotK x q) `mappend` f' bs q (x+1) -			      else f' bs q (x+1)+    where f' [] q _ = uhad q+          f' (b:bs) q x = if b then (rotK x q) `mappend` f' bs q (x+1) +                  else f' bs q (x+1)  --need to change this into a conQRec??? -- e.g. qft [Qbit 0]@@ -38,7 +38,8 @@  -- | A test of the QFT unitary, over a quantum integer initialised to \n\. tryQft :: Int -> QIO Int-tryQft n = do QInt qs <- mkQ n-	      applyU(qft qs)-	      x <- measQ (QInt qs)-	      return x+tryQft n = do +  QInt qs <- mkQ n+  applyU(qft qs)+  x <- measQ (QInt qs)+  return x

@@ -12,6 +12,8 @@ import QIO.Vec import QIO.VecEq import QIO.Heap+import Control.Applicative (Applicative(..))+import Control.Monad       (liftM, ap)  -- | A "Pure" state can be thought of as a vector of classical basis states, stored -- as Heaps, along with complex amplitudes.@@ -46,10 +48,10 @@ -- a Unitary that applies the rotation to the given qubit. uMatrix :: Qbit -> (CC,CC,CC,CC) -> Unitary uMatrix q (m00,m01,m10,m11) = U (\ fv h -> (if (fromJust(h ? q)) -                                           then   (m01 <*> (unEmbed $ return (update h q False))) -                                                  <+> (m11 <*> (unEmbed $ return h)) -                                           else   (m00 <*> (unEmbed $ return h)) -                                                  <+> (m10 <*> (unEmbed $ return (update h q True)))))+                                           then   (m01 <.> (unEmbed $ return (update h q False))) +                                                  <+> (m11 <.> (unEmbed $ return h)) +                                           else   (m00 <.> (unEmbed $ return h)) +                                                  <+> (m10 <.> (unEmbed $ return (update h q True)))))  -- | A rotation can be converted into a "Unitary", using the 'uMatrix' function uRot :: Qbit -> Rotation -> Unitary@@ -85,7 +87,7 @@  -- | A quantum state is a defined as the next free qubit reference, along with the -- Pure state that represents the overall quantum state-data StateQ = StateQ { free :: Int, pure :: Pure }+data StateQ = StateQ { free :: Int, pureState :: Pure }  -- | The initial 'StateQ' initialStateQ :: StateQ@@ -129,9 +131,16 @@ instance Show a => Show (Prob a) where     show (Prob (Vec ps)) = show (filter (\ (a,p) -> p>0) ps) +instance Functor Prob where+    fmap = liftM+ +instance Applicative Prob where+    pure = Prob . return+    (<*>) = ap+ -- | Prob forms a Monad instance Monad Prob where-    return = Prob . return+    return = pure     (Prob ps) >>= f = Prob (ps >>= unProb . f)  -- | Prob is also a PMonad, where the result of both computations are combined into

@@ -6,6 +6,8 @@  import Data.Monoid as Monoid import Data.Complex+import Control.Applicative (Applicative(..))+import Control.Monad       (liftM, ap)  -- | For Real numbers, we simply use the built in Double type type RR = Double@@ -78,10 +80,16 @@ uphase :: Qbit -> RR -> U uphase x r = rot x (rphase r)  +instance Functor QIO where+    fmap = liftM+ +instance Applicative QIO  where+    pure  = QReturn+    (<*>) = ap  -- | The "QIO" type forms a Monad instance Monad QIO where-    return = QReturn+    return = pure     (QReturn a) >>= f = f a     (MkQbit b g) >>= f = MkQbit b (\ x -> g x >>= f)     (ApplyU u q) >>= f = ApplyU u (q >>= f)@@ -140,7 +148,7 @@  -- | We can display a qubit reference instance Show Qbit where-	show (Qbit q) = "(Qbit:" ++ show q ++ ")"+    show (Qbit q) = "(Qbit:" ++ show q ++ ")"  -- | We can display the matrix representation of a rotation instance Show Rotation where@@ -148,7 +156,7 @@  -- | We can display a representation of a unitary instance Show U where-	show u = show' u 0 (-1)+    show u = show' u 0 (-1)  -- | A helper function for the show instance of U show' :: U -> Int -> Int -> String

@@ -0,0 +1,298 @@+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}++-- | This module defines the Syntax of the Quantum IO Monad, which is+-- an embedded language for writing quantum computations. It is an+-- alternative definition using the approach of defining F-Algebras.+module QIO.QioSynAlt where++import Data.Monoid as Monoid+import Data.Complex+import Control.Applicative (Applicative(..))+import Control.Monad       (liftM, ap)++-- | For Real numbers, we simply use the built in Double type+type RR = Double++-- | For Complex numbers, we use the built in Complex numbers, over our Real +-- number type (i.e. Double)+type CC = Complex RR++-- | The amplitude of a complex number is the magnitude squared.+amp :: CC -> RR+amp k = (magnitude k)*(magnitude k)++-- | The type of Qubits in QIO are simply integer references.+newtype Qbit = Qbit Int deriving (Num, Enum, Eq, Ord)++-- | We can display a qubit reference+instance Show Qbit where+    show (Qbit q) = "(Qbit:" ++ show q ++ ")"++-- | A rotation is in essence a two-by-two complex valued matrix+type Rotation = ((Bool,Bool) -> CC)++-- | We can display the matrix representation of a rotation+instance Show Rotation where+    show f = "(" ++ (show (f (False,False))) ++ "," ++ (show (f (False,True))) ++ "," ++ (show (f (True,False))) ++ "," ++ (show (f (True,True))) ++ ")"++-- | The identity rotation+rid :: Rotation+rid (x,y) = if x==y then 1 else 0++-- | The not rotation+rnot :: Rotation+rnot (x,y) = if x==y then 0 else 1++-- | The hadamard rotation+rhad :: Rotation+rhad (x,y) = if x && y then -h else h where h = (1/sqrt 2)++-- | The phase rotation+rphase :: RR -> Rotation+rphase _ (False,False)  = 1+rphase r (True,True)    = exp(0:+r)+rphase _ (_,_)          = 0++-- | Returns the inverse (or reverse) of the given rotation+rrev :: Rotation -> Rotation+rrev r (False,True)   = conjugate (r (True,False)) +rrev r (True,False)   = conjugate (r (False,True))+rrev r xy             = conjugate (r xy)++-- | Rotations can be compared for equality.+-- They are equal if the define the same matrix.+instance Eq Rotation where+    f == g =    (f (False,False)  == g (False,False)) +             && (f (False,True)   == g (False,True)) +             && (f (True,False)   == g (True,False)) +             && (f (True,True)    == g (True,True))+    f /= g =    (f (False,False)  /= g (False,False)) +             || (f (False,True)   /= g (False,True)) +             || (f (True,False)   /= g (True,False)) +             || (f (True,True)    /= g (True,True))++-- | The non-recursive data type definition of a unitary operation+data UFunctor u = UReturn +                | Rot Qbit Rotation u+                | Swap Qbit Qbit u +                | Cond Qbit (Bool -> u) u +                | Ulet Bool (Qbit -> u) u++-- | In order to define an F-Algebra, 'UFunctor' must be a functor.+instance Functor UFunctor where+    fmap eval UReturn = UReturn+    fmap eval (Rot q r u) = Rot q r (eval u)  +    fmap eval (Swap q1 q2 u) = Swap q1 q2 (eval u)+    fmap eval (Cond q f u) = Cond q (eval . f) (eval u)+    fmap eval (Ulet b f u) = Ulet b (eval . f) (eval u)++-- | The fix point type construtor.+newtype Fix f = Fx (f (Fix f))++-- | We can define the inverse of Fx+unFix :: Fix f -> f (Fix f)+unFix (Fx x) = x++-- | We fix the non-recursice data-type in order to get our type 'U'+-- of unitary operations.+type U = Fix UFunctor++-- | The type of an F-Algebra.+type Algebra f a = f a -> a++-- | The type of the initial algebra for UFunctor+type UInitialAlgebra = Algebra UFunctor U++-- | We can now define the initial algebra for U+uInitialAlgebra :: UInitialAlgebra+uInitialAlgebra = Fx++-- | We can use a catamorphism to abstract evaluation over a given+-- algebra+cata :: Functor f => Algebra f a -> Fix f -> a+cata algebra = algebra . fmap (cata algebra) . unFix++-- | The type "U" forms a Monoid. +instance Monoid U where+    mempty = Fx UReturn+    mappend (Fx UReturn) u = u+    mappend (Fx (Rot x a u)) u' = Fx $ Rot x a (mappend u u')+    mappend (Fx (Swap x y u)) u' = Fx $ Swap x y (mappend u u')+    mappend (Fx (Cond x br u')) u'' = Fx $ Cond x br (mappend u' u'')+    mappend (Fx (Ulet b f u)) u' = Fx $ Ulet b f (mappend u u')++-- | Apply the given rotation to the given qubit+rot :: Qbit -> Rotation -> U+rot x r = Fx $ Rot x r mempty++-- | Swap the state of the two given qubits+swap :: Qbit -> Qbit -> U+swap x y = Fx $ Swap x y mempty++-- | Apply the conditional unitary, depending on the value of the given qubit+cond :: Qbit -> (Bool -> U) -> U+cond x br = Fx $ Cond x br mempty++-- | Introduce an Ancilla qubit in the given state, for use in the sub-unitary+ulet :: Bool -> (Qbit -> U) -> U+ulet b ux = Fx $ Ulet b ux mempty++-- | Apply a not rotation to the given qubit+unot :: Qbit -> U+unot x = rot x rnot++-- | Apply a hadamard rotation to the given qubit+uhad :: Qbit -> U+uhad x = rot x rhad++-- | Apply a phase rotation (of the given angle) to the given qubit+uphase :: Qbit -> RR -> U+uphase x r = rot x (rphase r) ++-- | Returns the inverse (or reverse) of the given unitary operation,+-- using an F-Algebra+urev :: U -> U+urev = cata urev_algebra+ where +   urev_algebra :: UFunctor U -> U+   urev_algebra UReturn = Fx UReturn+   urev_algebra (Rot x r u) = u `mappend` rot x (rrev r)+   urev_algebra (Swap x y u) = u `mappend` swap x y+   urev_algebra (Cond x br u) = u `mappend` cond x br+   urev_algebra (Ulet b xu u) = u `mappend` ulet b xu++-- | We can display a representation of a unitary, using an F-Algebra+instance Show U where+  show = cata showU_algebra+    where+      showU_algebra :: UFunctor String -> String+      showU_algebra UReturn = ""+      showU_algebra (Rot q r u) = +       "Rotate " ++ show q ++ ":" ++ show r ++ "\n" ++ u +      showU_algebra (Swap q1 q2 u) = +       "Swap " ++ show q1 ++ " and " ++ show q2 ++ "\n" ++ u+      showU_algebra (Cond q br u) = +       "Cond " ++ show q ++ " \\b -> if b then (\n"+       ++ unlines (map (' ':) (lines $ br True))+       ++ ") else (\n"+       ++ unlines (map (' ':) (lines $ br False))+       ++ ")\n" ++ u+      showU_algebra (Ulet b xu u) =+       let fv = find_fv xu in+       "Ulet " ++ show fv ++ " = " ++ (if b then "1" else "0") ++ " in (\n" +       ++ unlines (map (' ':) (lines $ xu fv))+       ++ ")\n" ++ u+      -- this is currently a dummy function+      find_fv :: (Qbit -> String) -> Qbit+      find_fv _ = -1      ++-- | The non-recursive data type definition of a QIO computation+data QIOFunctor a q = QReturn a +                    | MkQbit Bool (Qbit -> q) +                    | ApplyU U q+                    | Meas Qbit (Bool -> q)++-- | In order to define an F-Algebra, 'UF' must be a functor.+instance Functor (QIOFunctor a) where+    fmap eval (QReturn a) = QReturn a+    fmap eval (MkQbit b f) = MkQbit b (eval . f)+    fmap eval (ApplyU u q) = ApplyU u (eval q)+    fmap eval (Meas q f) = Meas q (eval . f)++-- | We fix the non-recursice data-type in order to get our type 'U'+-- of unitary operations.+type QIOprim a = Fix (QIOFunctor a)++-- | The type of the initial algebra for UFunctor+type QIOInitialAlgebra a = Algebra (QIOFunctor a) (QIOprim a)++-- | We can now define the initial algebra for U+qioInitialAlgebra :: QIOInitialAlgebra a+qioInitialAlgebra = Fx+  +-- | The "QIO" type forms a Monad, by wrapping 'QIOprim'+data QIO a = Apply (Fix (QIOFunctor a))++-- | We can remove the wrapper.+primQIO :: QIO a -> QIOprim a+primQIO (Apply q) = q++instance Functor QIO where+    fmap = liftM+ +instance Applicative QIO  where+    pure  = Apply . Fx . QReturn+    (<*>) = ap++-- | The wrapper type 'ApplyFix' forms a Monad+instance Monad QIO where+    return = pure+    (Apply (Fx (QReturn a))) >>= f = f a+    (Apply (Fx (MkQbit b g))) >>= f = Apply . Fx $ +      MkQbit b  (\q -> primQIO $ (Apply (g q)) >>= f)+    (Apply (Fx (ApplyU u q))) >>= f = Apply . Fx $ +      ApplyU u $ primQIO (Apply q >>= f)+    (Apply (Fx (Meas x g))) >>= f = Apply . Fx $ +      Meas x (\b -> primQIO $ (Apply (g b)) >>= f)++-- | Initialise a qubit in the given state (adding it to the overall quantum state)+mkQbit :: Bool -> QIO Qbit+mkQbit b = Apply . Fx $ MkQbit b (\q -> primQIO (return q))++-- | Apply the given unitary operation to the current quantum state+applyU :: U -> QIO ()+applyU u =  Apply . Fx $ ApplyU u $ primQIO (return ())++-- | Measure the given qubit, and return the measurement outcome (note that this+-- operation may affect the overall quantum state, as a measurement is destructive)+measQbit :: Qbit -> QIO Bool+measQbit x = Apply . Fx $ Meas x (\b -> primQIO (return b))++-- | We can show a QIO computation, using an F-Algebra+instance (Show a) => Show (QIO a) where+  show = (cata showQIO_algebra) . primQIO+   where+    showQIO_algebra :: (Show a) => Algebra (QIOFunctor a) String       +    showQIO_algebra (QReturn a) = +       "Return: " ++ show a ++ "\n"+    showQIO_algebra (MkQbit b f) = +       "Init" ++ (if b then "1" else "0") ++ "\n"+       ++ f 0+    showQIO_algebra (ApplyU u qio) = +       "Apply Unitary: (\n"+       ++ unlines (map (' ':) (lines $ show u))+       ++ ")\n" ++ qio+    showQIO_algebra (Meas q f) = +       "Measure " ++ show q ++ " \\b -> if b then (\n"+       ++ unlines (map (' ':) (lines $ f True))+       ++ ") else (\n"+       ++ unlines (map (' ':) (lines $ f False))+       ++ ")\n"++-- | We can count the number of each primitive operation using an F-Algebra+count :: QIO a -> (Int,Int,Int)+count = (cata count_algebra) . primQIO+ where+  count_algebra :: Algebra (QIOFunctor a) (Int,Int,Int)+  count_algebra (QReturn _) = (0,0,0)+  count_algebra (MkQbit b f) = let (mk,ap,ms) = f 0 in+                               (mk+1,ap,ms)+  count_algebra (ApplyU _ (mk,ap,ms)) = (mk,ap+1,ms)+  count_algebra (Meas q f) = let (mk,ap,ms) = f False in+                             (mk,ap,ms+1)+     +toffoli :: Qbit -> Qbit -> Qbit -> U+toffoli q1 q2 q3 = +  cond q1 (\b1 -> if b1 then (+   cond q2 (\b2 -> if b2 then (unot q3) +    else mempty)) else mempty)++and :: Bool -> Bool -> QIO Bool+and a b = do+ q1 <- mkQbit a+ q2 <- mkQbit b+ q3 <- mkQbit False+ applyU (toffoli q1 q2 q3)+ measQbit q3

@@ -24,20 +24,16 @@  -- | The overall \"phase-estimation\" structure of Shor's algorithm. shorU :: QInt -> QInt -> Int -> Int -> U-shorU k i1 x n = hadamardsI k-		 `mappend` -		 modExp n x k i1-                 `mappend`-		 qftI k+shorU k i1 x n = hadamardsI k `mappend` modExp n x k i1 `mappend` qftI k  -- | A quantum computation the implementes shor's algorithm, returning the period -- of the function. shor :: Int -> Int -> QIO Int-shor x n = do i0 <- mkQ 0-              i1 <- mkQ 1-	      applyU (shorU i0 i1 x n)-	      p <- measQ i0-	      return p+shor x n = do +  i0 <- mkQ 0+  i1 <- mkQ 1+  applyU (shorU i0 i1 x n)+  measQ i0  -- | A classical (inefficient) implementation of the period finding subroutine  period :: Int -> Int -> Int@@ -49,40 +45,41 @@          | otherwise = do x <- rand_coprime n                           a <- shor x n                           let xa = x^(half a) -		              in if   odd a || xa == (n-1) `mod` n || a == 0-		                 then factor n-		                 else  return (gcd (xa+1) n,gcd (xa-1) n)+                            in if odd a || xa == (n-1) `mod` n || a == 0+                               then factor n+                               else return (gcd (xa+1) n,gcd (xa-1) n) --this function can only be run too, for similar reasons to the rand_co' --function below  -- | This function simulates the running of a QIO computation, whilst using -- System.Time functions to time how long the simulation took. runTime :: QIO a -> IO a-runTime a = do start <- getClockTime-               result <- run a-               stop <- getClockTime-               putStr ("The total time taken was " ++ -                      (timeDiffToString (diffClockTimes stop start) ++ "\n"))-               return result+runTime a = do +  start <- getClockTime+  result <- run a+  stop <- getClockTime+  putStr ("The total time taken was " ++ (timeDiffToString (diffClockTimes stop start) ++ "\n"))+  return result  -- | Times the running of various subroutines within the factorisation algorithm. factorV' :: Int -> IO (Int,Int) factorV' n | even n = return (2,2)-           | otherwise = do start <- getClockTime-                            putStr ("Started at " ++ (show start) ++ "\n")-                            x <- run (rand_coprime n)-                            putStr ("Calling \"shor " ++ show x ++ " " ++ show n ++ "\"\n")-                            a <- run (shor x n)-                            stop <- getClockTime-                            putStr ("Shor took " ++ (timeDiffToString (diffClockTimes stop start)) ++ "\n")-                            putStr ("period a = " ++ show a)-                            let xa = x^(half a) -		               in do putStr (", giving xa = " ++ show xa ++ "\n")-                                     if odd a || xa == (n-1) `mod` n || (gcd (xa+1) n,gcd (xa-1) n) == (1,n) || (gcd (xa+1) n,gcd (xa-1) n) == (n,1) || (gcd (xa+1) n,gcd (xa-1) n) == (1,1)-		                      then do putStr "Recalling factorV\n"-                                              factorV' n-		                      else do putStr "Result: " -                                              return (gcd (xa+1) n,gcd (xa-1) n)+           | otherwise = do +  start <- getClockTime+  putStr ("Started at " ++ (show start) ++ "\n")+  x <- run (rand_coprime n)+  putStr ("Calling \"shor " ++ show x ++ " " ++ show n ++ "\"\n")+  a <- run (shor x n)+  stop <- getClockTime+  putStr ("Shor took " ++ (timeDiffToString (diffClockTimes stop start)) ++ "\n")+  putStr ("period a = " ++ show a)+  let xa = x^(half a) +    in do putStr (", giving xa = " ++ show xa ++ "\n")+          if odd a || xa == (n-1) `mod` n || (gcd (xa+1) n,gcd (xa-1) n) == (1,n) || (gcd (xa+1) n,gcd (xa-1) n) == (n,1) || (gcd (xa+1) n,gcd (xa-1) n) == (1,1)+           then do putStr "Recalling factorV\n"+                   factorV' n+           else do putStr "Result: " +                   return (gcd (xa+1) n,gcd (xa-1) n)  -- | Calls the 'factorV'', and times the overall factorisation. factorV :: Int -> IO ()@@ -109,8 +106,9 @@ -- computation that returns a random number between 2 and \n\, that is then -- returned if it is co-prime to \n\. rand_co' :: Int -> QIO Int-rand_co' n = do x <- randomQIO (2,n)-		if gcd x n == 1 then return x else rand_co' n+rand_co' n = do +  x <- randomQIO (2,n)+  if gcd x n == 1 then return x else rand_co' n --simulating this (with the sim function) gives rise to infinite paths in --the computation, e.g. each path where gcd x n /= 1. However, this function --can still be run (with the run function) always returning a single value.

@@ -1,9 +1,11 @@- -- | This module defines a Vector as a list of pairs.  -- In the context of QIO, a Vector is the type used to represent a probability -- distribution. module QIO.Vec where +import Control.Applicative (Applicative(..))+import Control.Monad       (liftM, ap)+ -- | A Vector over types 'x' and 'a' is a wrapper around list of  -- pairs of 'a' and 'x'. newtype Vec x a = Vec {unVec :: [(a,x)]} deriving Show@@ -26,8 +28,21 @@ (<++>) :: (Vec x a) -> (Vec x a) -> Vec x a (Vec as) <++> (Vec bs) = (Vec (as ++ bs)) +instance Num n => Functor (Vec n) where+    fmap = liftM+ +instance Num n => Applicative (Vec n)  where+    pure  a = Vec [(a,1)]+    (<*>) = ap+ -- | Vectors, over Numeric types, can be defined as a Monad. instance Num n => Monad (Vec n) where-    return a = Vec [(a,1)]+    return = pure     (Vec ms) >>= f = Vec [(b,i*j) | (a,i) <- ms, (b,j) <- unVec (f a)]    ++ ++++

@@ -8,6 +8,8 @@  import QIO.QioSyn import QIO.Heap+import Control.Applicative (Applicative(..))+import Control.Monad       (liftM, ap)  -- | Any type that fulfills this type class is a Vector over types with equality class VecEq v where@@ -16,7 +18,7 @@     -- | Two Vectors can be combined     (<+>) :: (Eq a, Num x) => v x a -> v x a -> v x a     -- | A Vector can be multiplied by a scalar-    (<*>) :: (Num x, Eq x) => x -> v x a -> v x a+    (<.>) :: (Num x, Eq x) => x -> v x a -> v x a     -- | The amplitude of a given element can be accessed     (<@>) :: (Eq a, Num x) => a -> v x a -> x     -- | The vector can be created from a list of pairs@@ -64,7 +66,7 @@ instance VecEq VecEqL where       vzero = vEqZero       (<+>) = vEqPlus-      (<*>) = vEqTimes+      (<.>) = vEqTimes       (<@>) = vEqAt       fromList as = VecEqL as       toList (VecEqL as) = as@@ -80,8 +82,8 @@     eqReturn a = fromList [(a,1)]     eqBind va f = case toList va of                    ([]) -> vzero-                   ((a,x):[]) -> x <*> f a-                   ((a,x):vas) -> (x <*> f a) <+> ((fromList vas) `eqBind` f)+                   ((a,x):[]) -> x <.> f a+                   ((a,x):vas) -> (x <.> f a) <+> ((fromList vas) `eqBind` f)  -- | We can define a datatype that holds EqMonad operations, so that it can -- be defined as a Monad. @@ -90,9 +92,16 @@    Return :: EqMonad m => a -> AsMonad m a    Bind   :: EqMonad m => AsMonad m a -> (a -> AsMonad m b) -> AsMonad m b  +instance EqMonad m => Functor (AsMonad m) where+    fmap = liftM+ +instance EqMonad m => Applicative (AsMonad m)  where+    pure  = Return+    (<*>) = ap+ -- | We can define an AsMonad over an EqMonad, as a Monad instance EqMonad m => Monad (AsMonad m) where-   return = Return+   return = pure    (>>=) = Bind  -- | Given Equality, we can unembed the EqMonad operations from an AsMonad