packages feed

QIO (empty) → 1.0

raw patch · 15 files changed

+1186/−0 lines, 15 filesdep +basedep +containersdep +haskell98setup-changed

Dependencies added: base, containers, haskell98, mtl, old-time, random

Files

+ LICENSE view
@@ -0,0 +1,10 @@+Copyright (c) 2010, 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:++    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.+    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.+
+ QIO.cabal view
@@ -0,0 +1,18 @@+Name:		QIO+Version:	1.0+Cabal-Version:	>= 1.2+License:	BSD3+License-File:	LICENSE+Author:		Alexander S. Green+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+Build-Type:	Simple++Library+  Build-Depends:		+    base >= 4 && < 5, containers, haskell98, 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/Heap.hs view
@@ -0,0 +1,27 @@+{-# OPTIONS_GHC -fglasgow-exts #-}++module QIO.Heap where++import qualified Data.Map as Map+import Data.Maybe as Maybe+import QIO.QioSyn++class Eq h => Heap h where+      initial :: h+      update :: h -> Qbit -> Bool -> h+      (?) :: h -> Qbit -> Maybe Bool+      forget :: h -> Qbit -> h+      hswap :: h -> Qbit -> Qbit -> h+      hswap h x y = update (update h y (fromJust (h ? x))) x (fromJust (h ? y)) +type HeapMap = Map.Map Qbit Bool++instance Heap HeapMap where+       initial = Map.empty+       update h q b = Map.insert q b h+       h ? q = Map.lookup q h+       forget h q = Map.delete q h+++    ++
+ QIO/QArith.hs view
@@ -0,0 +1,184 @@+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}++module QIO.QArith where++import Data.Monoid as Monoid+import QIO.QioSyn+import QIO.Qdata+import QIO.QioClass+import QIO.Qio+import QIO.QExamples++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++ifElseQ :: Qbit -> U -> U -> U+ifElseQ qa t f = cond qa (\ qa -> if qa then t else f)++ifQ :: Qbit -> U -> U+ifQ qa t = ifElseQ qa t mempty++cnot :: Qbit -> Qbit -> U+cnot qa qb = ifQ qa (unot qb)++addBit :: Qbit -> Qbit -> Qbit -> U+addBit qc qa qb = +    cnot qa qb `mappend`+    cnot qc qb++carry :: Qbit -> Qbit -> Qbit -> Qbit -> U+carry qci qa qb qcsi = +    cond qci (\ ci ->+              cond qa (\ a ->+                       cond qb (\ b ->+                                if ci && a || ci && b || a && b+                                then unot qcsi+                                else mempty)))++addBits :: [Qbit] -> [Qbit] -> Qbit -> U+addBits qas qbs qc' = +    letU False (addBits' qas qbs)+    where addBits' [] [] qc = ifQ qc (unot qc')+          addBits' (qa:qas) (qb:qbs) qc =+              letU False (\ qc' -> carry qc qa qb qc' `mappend`+                                   addBits' qas qbs qc'`mappend`+                                   urev (carry qc qa qb qc')) `mappend`+              addBit qc qa qb++addBits' :: [Qbit] -> [Qbit] -> [Qbit] -> Qbit -> U+addBits' [] [] [] qc = mempty+addBits' (qa:qas) (qb:qbs) (qc':qcs') qc =+    (carry qc qa qb qc' `mappend`+     addBits' qas qbs qcs' qc'`mappend`+     urev (carry qc qa qb qc')) `mappend`+    addBit qc qa qb++adder :: QInt -> QInt -> Qbit -> U+adder (QInt qas) (QInt qbs) qc = addBits qas qbs qc ++tadder :: (Int,(Int,Bool)) -> QIO (Int,(Int,Bool))+tadder xyc = do q @ (qx,(qy,qc)) <- mkQ xyc+                applyU (adder qx qy qc)+                xyc <- measQ q+                return xyc++tRadder :: (Int,(Int,Bool)) -> QIO (Int,(Int,Bool))+tRadder xyc = do q @ (qx,(qy,qc)) <- mkQ xyc+                 applyU (urev (adder qx qy qc))+                 xyc <- measQ q+                 return xyc++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++adderMod :: Int -> QInt -> QInt -> U+adderMod n qa qb =+    letU n (\ qn ->+       letU False (\ qz ->+          letU False (\ qc -> +             adder qa qb qc+             `mappend` -- b = a+b, c=False+             urev (adder qn qb qc)+             `mappend` -- b = a+b-N, c = a+b < N+             cond qc (\ c -> if c then unot qz else mempty)+             `mappend` -- z = c = a+b < N+             cond qz (\ z -> if z then adder qn qb qc else mempty)+             `mappend` -- b = a+b mod N, c = False, z = a+b < N+             urev (adder qa qb qc)+             `mappend` -- if a+b < N then a=a,b=b,c=False +                       -- else a=a,b=a+b mod N - b,c=True+                       -- z = not c+             cond qc (\ c -> if c then mempty else unot qz)+             `mappend` -- z = False+             adder qa qb qc))) -- b = a+b mod N, c=False, z=False++tadderMod :: Int -> (Int,Int) -> QIO (Int,Int)+tadderMod n ab = do q @ (qa,qb) <- mkQ ab+                    applyU (adderMod n qa qb)+                    ab <- measQ q+                    return ab++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)))+		               +-- output is a*x mod n+tmultMod :: Int -> Int -> Int -> QIO (Int,Int)+tmultMod n a x = do y <- mkQ 0+                    x' <- mkQ x+                    applyU(multMod n a x' y)+                    qy <- measQ y+                    qx <- measQ x'+                    return (qx,qy)++condMultMod :: Qbit -> Int -> Int -> QInt -> QInt -> U+condMultMod q n a x y = ifQ q (multMod n a x y)++------------------------------------------------------------------------------++inverseMod :: Int -> Int -> Int+inverseMod n a = inverseMod'' n a (inverseMod' n a)++inverseMod' :: Int -> Int -> [Int]+inverseMod' n a = [x | x <- [1..n], ((x*a) `mod` n) == 1]+++inverseMod'' :: Int -> Int -> [Int] -> Int+inverseMod'' n a [] = error ("inverseMod: no inverse of "++(show a)++" mod "++(show n)++ " found")+inverseMod'' _ _ xs = head xs++-------------------------------------------------------------------------------++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++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++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))++--a^x mod N++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++++		       ++++++++                     +++++
+ QIO/QExamples.hs view
@@ -0,0 +1,138 @@+module QIO.QExamples where++import Data.Monoid as Monoid+import QIO.QioSyn+import QIO.Qdata+import QIO.QioClass+import QIO.Qio+++q0 :: QIO Qbit+q0 = mkQ False++q1 :: QIO Qbit+q1  =  mkQ True+          +qPlus :: QIO Qbit+qPlus  =  do  qa <- q0+              applyU (uhad qa)+              return qa++qMinus :: QIO Qbit+qMinus  =  do  qa <- q1+               applyU (uhad qa)+               return qa++randBit :: QIO Bool+randBit  =  do  qa <- qPlus+                x <- measQbit qa+                return x++share :: Qbit -> QIO Qbit+share  qa  =  do  qb <- q0+                  applyU (cond qa (\a -> if a then (unot qb)+                                              else (mempty)  )  )+                  return qb++bell :: QIO (Qbit, Qbit)+bell = do  qa <- qPlus+           qb <- share qa+           return (qa,qb)++test_bell :: QIO (Bool,Bool)+test_bell  =  do  qb <- bell+                  b <- measQ qb+                  return b++hadTwice :: Bool -> QIO Bool+hadTwice x = do q <- mkQ x+		applyU (uhad q `mappend` uhad q)+		b <- measQ q+		return b++hadTwice' :: Bool -> QIO Bool+hadTwice' x = do q <- mkQ x+	 	 applyU (uhad q)+                 applyU (uhad q)+		 b <- measQ q+		 return b++----------------------------------------------+---- Teleportation ---------------------------+----------------------------------------------++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++uZZ :: Qbit -> U+uZZ qb = (uphase qb pi)++bobsU :: (Bool,Bool) -> Qbit -> U+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))++bob :: Qbit -> (Bool,Bool) -> QIO Qbit+bob eq cd  =  do  applyU (bobsU cd eq)+                  return eq++teleportation :: Qbit -> QIO Qbit+teleportation iq  =  do  (eq1,eq2) <- bell+                         cd <- alice iq eq1+                         tq <- bob eq2 cd+                         return tq++test_teleport :: QIO (Bool,Bool)+test_teleport = do (q1,q2) <- bell+                   tq2 <- teleportation q2+		   result <- measQ (q1,tq2)+		   return result++teleport_true' :: QIO Qbit+teleport_true' = do q <- q1+                    tq <- teleportation q+                    return tq++teleport_true :: QIO Bool+teleport_true = do q <- teleport_true'+                   result <- measQ q+                   return result++teleport_random' :: QIO Qbit+teleport_random' = do q <- qPlus+                      tq <- teleportation q+                      return tq++teleport_random :: QIO Bool+teleport_random = do q <- teleport_random'+                     result <- measQ q+                     return result++-----------------------------------------------+----- Deutsch's Algorithm ---------------------+-----------------------------------------------++u :: (Bool -> Bool) -> Qbit -> Qbit -> U+u f x y = cond x (\ b -> if f b then unot y else mempty)++deutsch :: (Bool -> Bool) -> QIO Bool+deutsch f = do+              x <- qPlus+              y <- qMinus+              applyU (u f x y)+              applyU (uhad x)+              measQ x++-----------------------------------------------++problem :: QIO Bool+problem = do q <- qPlus+             x <- measQ q+             if x then return x else problem+-- can be run returning True, but cannot be simulated!
+ QIO/QIORandom.hs view
@@ -0,0 +1,114 @@+{-# 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++
+ QIO/Qdata.hs view
@@ -0,0 +1,71 @@+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}++module QIO.Qdata where++import Data.Monoid as Monoid+import QIO.QioSyn++class Qdata a qa | a -> qa, qa -> a where+    mkQ :: a -> QIO qa+    measQ :: qa -> QIO a+    letU :: a -> (qa -> U) -> U+    condQ :: qa -> (a -> U) -> U++instance Qdata Bool Qbit where+    mkQ = mkQbit+    measQ = measQbit+    letU b xu = ulet b xu+    condQ q br = cond q br+                         +instance (Qdata a qa,Qdata b qb) => Qdata (a,b) (qa,qb) where++    mkQ (a,b) = do qa <- mkQ a+                   qb <- mkQ b+                   return (qa,qb)++    measQ (qa,qb) = do a <- measQ qa+                       b <- measQ qb+                       return (a,b)++    letU (a,b) xyu = letU a (\ x -> letU b (\ y -> xyu (x,y)))++    condQ (qa,qb) br = condQ qa (\x -> condQ qb (\y -> br (x,y))) ++instance Qdata a qa => Qdata [a] [qa] where+    mkQ n = sequence (map mkQ n)+    measQ qs = sequence (map measQ qs)+    letU as xsu = letU' as []+        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]))++condQRec :: Qdata a qa => [qa] -> [(a -> U)] -> U+condQRec [] [] = mempty+condQRec (q:qs) (u:us) = (condQ q u) `mappend` condQRec qs us+++qIntSize :: Int+qIntSize = 4++newtype QInt = QInt [Qbit] deriving Show++int2bits :: Int -> [Bool]+int2bits n = int2bits' n qIntSize+             where int2bits' 0 0 = []+                   int2bits' _ 0 = error "int2bits: too large"+                   int2bits' n l = ((n `mod` 2) /= 0) : int2bits' (n `div` 2) (l-1)++bits2int :: [Bool] -> Int+bits2int [] = 0+bits2int (b:bs) = (2*bits2int bs)+(if b then 1 else 0)++instance Qdata Int QInt where+    mkQ n = do qn <- mkQ (int2bits n)+               return (QInt qn)+    measQ (QInt qbs) = +        do bs <- measQ qbs+           return (bits2int bs)+    letU n xu = letU (int2bits n) (\ bs -> xu (QInt bs))+    condQ (QInt qi) qiu = condQ qi (\ x -> qiu (bits2int x))
+ QIO/Qft.hs view
@@ -0,0 +1,57 @@+module QIO.Qft where++import Data.Monoid as Monoid+import QIO.QioSyn+import QIO.Qio+import QIO.Qdata++qft :: [Qbit] -> U+qft qs = condQ qs (\bs -> qftAcu qs bs [])+++qftAcu :: [Qbit] -> [Bool] -> [Bool] -> U+qftAcu [] [] _ = mempty+qftAcu (q:qs) (b:bs) cs = qftBase cs q `mappend` qftAcu qs bs (b:cs)++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)++--need to change this into a conQRec???++++-- e.g. qft [Qbit 0]+-- = condQ [Qbit 0] (\(b:bs) -> uhad 0 `mappend` mempty)+-- but gives  cond 0 (\x -> if x then uhad 0 else uhad 0) which is forbidden++testCond :: [Qbit] -> U+testCond [] = mempty+testCond (q:qs) = condQ (q:qs) (\bs -> uhad q)++testCondOk :: [Qbit] -> U+testCondOk [] = mempty+testCondOk (q:qs) = condQ (qs) (\bs -> uhad q)++rotK :: Int -> Qbit -> U+rotK k q = uphase q (1.0/(2.0^k))++tryQft :: Int -> QIO Int+tryQft n = do QInt qs <- mkQ n+	      applyU(qft qs)+	      x <- measQ (QInt qs)+	      return x++tC :: (Qbit,Qbit) -> U+tC qxy = condQ qxy (\xy -> tC' qxy xy)++tC' :: (Qbit,Qbit) -> (Bool,Bool) -> U+tC' (qx,qy) (x,y) = if x then unot qy else mempty++testTC :: QIO (Bool,Bool)+testTC = do (qx,qy) <- mkQ (False,False)+	    applyU (uhad qx)+            applyU (tC (qx,qy))+            measQ (qx,qy)
+ QIO/Qio.hs view
@@ -0,0 +1,133 @@+module QIO.Qio where++import List+import qualified System.Random as Random+import Data.Monoid as Monoid+import Data.Maybe as Maybe+import Control.Monad.State+import QIO.QioSyn+import QIO.Vec+import QIO.VecEq+import QIO.Heap++type Pure = VecEqL CC HeapMap++updateP :: Pure -> Qbit -> Bool -> Pure+updateP p x b = VecEqL (map (\ (h,pa) -> (update h x b,pa)) (unVecEqL p))++newtype Unitary =  U {unU :: Int -> HeapMap -> Pure }++instance Monoid Unitary where+    mempty = U (\ fv h -> unEmbed $ return h)+    mappend (U f) (U g) = U (\ fv h -> unEmbed $ do h' <- Embed $ f fv h+                                                    h'' <- Embed $ g fv  h'+                                                    return h''+                                                 )++uRot :: Qbit -> Rotation -> Unitary+uRot q r = if (unitaryRot r) then (uMatrix q (r (False,False),+                                              r (False,True),+                                              r (True,False),+                                              r (True,True)))+                             else error "Non unitary Rotation!"++unitaryRot :: Rotation -> Bool+unitaryRot r = True+-- update to check that the rotation is unitary...++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)))))++uSwap :: Qbit -> Qbit -> Unitary+uSwap x y = U (\ fv h -> unEmbed $ return (hswap h x y ))++uCond :: Qbit -> (Bool -> Unitary) -> Unitary+--uCond x us = U (\ fv h -> updateP (unU (us (h ? x)) fv (forget h x)) x (h ? x))+uCond x us = U (\ fv h -> unU (us (fromJust(h ? x))) fv h )+--whether or not to forget? (if not then no runtime error for conditionals)++uLet :: Bool -> (Qbit -> Unitary) -> Unitary+uLet b ux = U (\fv h -> unU (ux (Qbit fv)) (fv + 1) (update h (Qbit fv) b))+--doesn't enforce unitary+-- need Unitary -> [Qbit] ???++runU :: U -> Unitary+runU UReturn = mempty+runU (Rot x a u) = uRot x a `mappend` runU u+runU (Swap x y u) = uSwap x y `mappend` runU u+runU (Cond x us u) = uCond x (runU.us) `mappend` runU u+runU (Ulet b xu u) = uLet b (runU.xu) `mappend` runU u++data StateQ = StateQ { free :: Int, pure :: Pure }++initialStateQ :: StateQ+initialStateQ = StateQ 0 (unEmbed $ return initial)++pa :: Pure -> RR+pa (VecEqL as) = foldr (\ (_,k) p -> p + amp k) 0 as++data Split = Split { p :: RR, ifTrue,ifFalse :: Pure }++split :: Pure -> Qbit -> Split+split (VecEqL as) x =+    let pas = pa (VecEqL as)+        (ift',iff') = partition (\ (h,_) -> (fromJust(h ? x))) as+        ift = VecEqL ift'+        iff = VecEqL iff'+        p_ift = if pas==0 then 0 else (pa ift)/pas+    in Split p_ift ift iff++class Monad m => PMonad m where+    merge :: RR -> m a -> m a -> m a++instance PMonad IO where+    merge pr ift iff = do pp <- Random.randomRIO (0,1.0)+                          if pr > pp then ift else iff++data Prob a = Prob {unProb :: Vec RR a}++instance Show a => Show (Prob a) where+    show (Prob (Vec ps)) = show (filter (\ (a,p) -> p>0) ps)++instance Monad Prob where+    return = Prob . return+    (Prob ps) >>= f = Prob (ps >>= unProb . f)++instance PMonad Prob where+    merge pr (Prob ift) (Prob iff) = Prob ((pr <**> ift) <++> ((1-pr) <**> iff))+++evalWith :: PMonad m => QIO a -> State StateQ (m a)+evalWith (QReturn a) = return (return a)+evalWith (MkQbit b g) = do (StateQ f p) <- get +                           put (StateQ (f+1) (updateP p (Qbit f) b))+                           evalWith (g (Qbit f))+evalWith (ApplyU u q) = do (StateQ f p) <- get+                           put (StateQ f (unEmbed $ do x <- Embed $ p+                                                       x' <-Embed $ uu f x+                                                       return x'+                                          )+                                )+                           evalWith q  +                               where U uu = runU u+evalWith (Meas x g) = do (StateQ f p) <- get+                         (let Split pr ift iff = split p x+                          in if pr < 0 || pr > 1 then error "pr < 0 or >1" +                             else do put (StateQ f ift)+                                     pift <- evalWith (g True)+                                     put (StateQ f iff)+                                     piff <- evalWith (g False)+                                     return (merge pr pift piff))++eval :: PMonad m => QIO a -> m a+eval p = evalState (evalWith p) initialStateQ++run :: QIO a -> IO a+run = eval++sim :: QIO a -> Prob a+sim = eval
+ QIO/QioClass.hs view
@@ -0,0 +1,61 @@+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}++module QIO.QioClass where++import Data.Maybe as Maybe+import Data.Monoid as Monoid+import Control.Monad.State+import QIO.QioSyn+import QIO.Heap+import Complex++newtype UnitaryC =  U {unU :: Int -> HeapMap -> HeapMap}++instance Monoid UnitaryC where+    mempty = U (\ fv bs -> bs)+    mappend (U f) (U g) = U (\ fv h -> g fv (f fv h))++uRotC :: Qbit -> Rotation -> UnitaryC+uRotC x f | f==rnot    = U (\ _ h -> update h x (not (fromJust (h ? x))))+          | f==rid     = mempty+          | otherwise  = error "not classical"++uSwapC :: Qbit -> Qbit -> UnitaryC +uSwapC x y = U (\ _ h -> hswap h x y )++uCondC :: Qbit -> (Bool -> UnitaryC) -> UnitaryC+uCondC x br = U (\ fv h -> update (unU (br (fromJust (h ? x))) fv (forget h x)) x (fromJust (h ? x)))++uLetC :: Bool -> (Qbit -> UnitaryC) -> UnitaryC+uLetC b ux = U (\ fv h -> unU (ux (Qbit fv)) (fv+1) (update h (Qbit fv) b))++runUC :: U -> UnitaryC+runUC UReturn = mempty+runUC (Rot x r u) = uRotC x r `mappend` runUC u+runUC (Swap x y u) = uSwapC x y `mappend` runUC u+runUC (Cond x us u) = uCondC x (runUC.us) `mappend` runUC u+runUC (Ulet b xu u) = uLetC b (runUC.xu) `mappend` runUC u++data StateC = StateC {fv :: Int, heap :: HeapMap}++initialStateC :: StateC +initialStateC = StateC 0 initial++runQStateC :: QIO a -> State StateC a+runQStateC (QReturn a) = return a+runQStateC (MkQbit b xq) = do (StateC fv h) <- get +                              put (StateC (fv+1) (update h (Qbit fv) b)) +                              runQStateC (xq (Qbit fv))+runQStateC (ApplyU u q) = do (StateC fv h) <- get+                             put (StateC fv (unU (runUC u) fv h))+                             runQStateC q +runQStateC (Meas x qs) = do (StateC _ h) <- get+                            runQStateC (qs (fromJust (h ? x)))++runC :: QIO a -> a+runC q = evalState (runQStateC q) initialStateC++          +++
+ QIO/QioSyn.hs view
@@ -0,0 +1,146 @@+{-# OPTIONS_GHC -fglasgow-exts #-}++module QIO.QioSyn where++import Data.Monoid as Monoid+import Complex+++-- complex numbers+type CC = Complex RR++amp :: CC -> RR+amp k = (magnitude k)*(magnitude k)+++-- real numbers+type RR = Float+++-- Qubits are references+newtype Qbit = Qbit Int deriving (Num, Enum, Eq, Ord)+++-- QIO and U as traces+type Rotation = ((Bool,Bool) -> CC)++data U = UReturn | Rot Qbit Rotation U+       | Swap Qbit Qbit U | Cond Qbit (Bool -> U) U | Ulet Bool (Qbit -> U) U++data QIO a = QReturn a | MkQbit Bool (Qbit -> QIO a) | ApplyU U (QIO a) +           | Meas Qbit (Bool -> QIO a)+++-- U functions+instance Monoid U where+    mempty = UReturn+    mappend UReturn u = u+    mappend (Rot x a u) u' = Rot x a (mappend u u')+    mappend (Swap x y u) u' = Swap x y (mappend u u')+    mappend (Cond x br u') u'' = Cond x br (mappend u' u'')+    mappend (Ulet b f u) u' = Ulet b f (mappend u u')++rot :: Qbit -> Rotation -> U+rot x r = Rot x r UReturn++swap :: Qbit -> Qbit -> U+swap x y = Swap x y UReturn++cond :: Qbit -> (Bool -> U) -> U+cond x br = Cond x br UReturn++ulet :: Bool -> (Qbit -> U) -> U+ulet b ux = Ulet b ux UReturn++urev :: U -> U+urev UReturn = UReturn+urev (Rot x r u) = urev u `mappend` rot x (rrev r)+urev (Swap x y u) = urev u `mappend` swap x y+urev (Cond x br u) = urev u `mappend` cond x (urev.br)+urev (Ulet b xu u) = urev u `mappend` ulet b (urev.xu)++unot :: Qbit -> U+unot x = rot x rnot++uhad :: Qbit -> U+uhad x = rot x rhad++uphase :: Qbit -> RR -> U+uphase x r = rot x (rphase r) +++--- QIO functions+instance Monad QIO where+    return = QReturn+    (QReturn a) >>= f = f a+    (MkQbit b g) >>= f = MkQbit b (\ x -> g x >>= f)+    (ApplyU u q) >>= f = ApplyU u (q >>= f)+    (Meas x g) >>= f = Meas x (\ b -> g b >>= f)++mkQbit :: Bool -> QIO Qbit+mkQbit b = MkQbit b return++applyU :: U -> QIO ()+applyU u =  ApplyU u (return ())++measQbit :: Qbit -> QIO Bool+measQbit x = Meas x return+++-- rotations+rid :: Rotation+rid (x,y) = if x==y then 1 else 0++rnot :: Rotation+rnot (x,y) = if x==y then 0 else 1++rhad :: Rotation+rhad (x,y) = if x && y then -h else h where h = (1/sqrt 2)++rphase :: RR -> Rotation+rphase _ (False,False)  = 1+rphase r (True,True)    = exp(0:+r)+rphase _ (_,_)          = 0++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)++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))+++-- show functions (for Qbit, Rotation and U)+instance Show Qbit where+	show (Qbit q) = "(Qbit:" ++ show q ++ ")"++instance Show Rotation where+    show f = "(" ++ (show (f (False,False))) ++ "," ++ (show (f (False,True))) ++ "," ++ (show (f (True,False))) ++ "," ++ (show (f (True,True))) ++ ")"++instance Show U where+	show u = show' u 0 (-1)++show' :: U -> Int -> Int -> String++show' (UReturn) x fv = ""++show' (Rot q a u) x fv = spaces x ++ "Rotate " ++ show q ++ " by " ++ show a ++ ".\n" ++ show' u x fv++show' (Swap q1 q2 u) x fv = spaces x ++ "Swap " ++ show q1 ++ " and " ++ show q2 ++ ".\n" ++ show' u x fv++show' (Cond q f u) x fv = spaces x ++ "Cond (if " ++ show q ++ " then \n" ++ spaces (x+1) ++ "(\n" ++ show' (f True) (x+1) fv ++ spaces (x+1) ++ ")\n" ++ spaces x ++ "else \n" ++ spaces (x+1) ++ "(\n" ++ show' (f False) (x+1) fv ++ spaces (x+1) ++ ")\n" ++ show' u x fv++show' (Ulet b f u) x fv = spaces x ++ "Ulet " ++ show b ++ " (\\" ++ show (Qbit fv) ++ "->\n " ++ show' (f (Qbit fv)) x (fv-1) ++ ")\n" ++ show' u x fv++spaces :: Int -> String+spaces 0 = ""+spaces (n+1) = "    " ++ spaces n   +
+ QIO/Shor.hs view
@@ -0,0 +1,115 @@+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}++module QIO.Shor where++import Data.Monoid as Monoid+import QIO.QIORandom+import QIO.QioSyn+import QIO.Qdata+import QIO.Qio+import QIO.QExamples+import QIO.QArith+import QIO.Qft+import System.Time++qftI :: QInt -> U+qftI (QInt i) = urev (qft i)++hadamardsI :: QInt -> U+hadamardsI (QInt xs) = hadamards xs++shorU :: QInt -> QInt -> Int -> Int -> U+shorU k i1 x n = hadamardsI k+		 `mappend` +		 modExp n x k i1+                 `mappend`+		 qftI k++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++period :: Int -> Int -> Int+period m q = r where (_,r) = reduce (m,q)+++factor :: Int -> QIO (Int,Int)+factor n | even n = return (2,2)+         | 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)+--this function can only be run too, for similar reasons to the rand_co'+--function below++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++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)++factorV :: Int -> IO ()+factorV n = do start <- getClockTime+               (a,b) <- factorV' n+               stop <- getClockTime+               putStr ( "Factors of "+                       ++ (show n) +                       ++ " include "+                       ++ (show a)+                       ++ " and "+                       ++ (show b)+                       ++ ".\n The total time taken was "+                       ++ (timeDiffToString (diffClockTimes stop start) ++ "\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+--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.++rand_coprime :: Int -> QIO Int+rand_coprime n = do x <- randomQIO (0,(length cps)-1)+                    return (cps!!x)+                     where cps = [x | x <- [0..n], gcd x n == 1]++++half :: Int -> Int+half x = floor (fromIntegral x/2.0)++reduce :: (Int,Int) -> (Int,Int)+reduce (x,y) = if g == 1 then (x,y) else (floor ((fromIntegral x)/(fromIntegral g)),floor ((fromIntegral y)/(fromIntegral g)))+               where g = gcd x y++++++
+ QIO/Vec.hs view
@@ -0,0 +1,22 @@+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}++module QIO.Vec where++newtype Vec x a = Vec {unVec :: [(a,x)]} deriving Show++empty :: Vec x a+empty = Vec []+ +(<@@>) :: (Num x,Eq a) => Vec x a -> a -> x+(Vec ms) <@@> a = foldr (\(b,k) m -> if a == b then m + k else m) 0 ms++(<**>) :: Num x => x -> (Vec x a) -> Vec x a+l <**> (Vec as) = (Vec (map (\ (a,k) -> (a,l*k)) as))++(<++>) :: (Vec x a) -> (Vec x a) -> Vec x a+(Vec as) <++> (Vec bs) = (Vec (as ++ bs))++instance Num n => Monad (Vec n) where+    return a = Vec [(a,1)]+    (Vec ms) >>= f = Vec [(b,i*j) | (a,i) <- ms, (b,j) <- unVec (f a)]+
+ QIO/VecEq.hs view
@@ -0,0 +1,87 @@+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}++module QIO.VecEq where++import QIO.QioSyn+import QIO.Heap++class VecEq v where+    vzero :: v x a+    (<+>) :: (Eq a, Num x) => v x a -> v x a -> v x a+    (<*>) :: (Num x) => x -> v x a -> v x a+    (<@>) :: (Eq a, Num x) => a -> v x a -> x+    fromList  :: [(a,x)] -> v x a+    toList    :: v x a -> [(a,x)] ++newtype VecEqL x a = VecEqL {unVecEqL :: [(a,x)]} deriving Show++vEqZero :: VecEqL x a+vEqZero = VecEqL []+          ++vEqPlus :: (Eq a, Num x) => VecEqL x a -> VecEqL x a -> VecEqL x a+(VecEqL as) `vEqPlus` vbs = foldr add vbs as+++vEqTimes :: (Num x) => x -> VecEqL x a -> VecEqL x a+l `vEqTimes` (VecEqL bs) | l==0 = VecEqL []+                         | otherwise = VecEqL (map (\ (b,k) -> (b,l*k)) bs)+          ++vEqAt :: (Eq a, Num x) => a -> VecEqL x a -> x+a `vEqAt` (VecEqL []) = 0+a `vEqAt` (VecEqL ((a',b):abs)) | a == a' = b+                                | otherwise = a `vEqAt` (VecEqL abs)+          ++add :: (Eq a,Num x) => (a,x) -> VecEqL x a -> VecEqL x a+add (a,x) (VecEqL axs) = VecEqL (addV' axs)+    where addV' [] = [(a,x)]+          addV' ((by @ (b,y)):bys) | a == b = (b,x+y):bys+                                   | otherwise = by:(addV' bys)++instance VecEq VecEqL where+      vzero = vEqZero+      (<+>) = vEqPlus+      (<*>) = vEqTimes+      (<@>) = vEqAt+      fromList as = VecEqL as+      toList (VecEqL as) = as++class EqMonad m where+    eqReturn :: Eq a => a -> m a+    eqBind   :: (Eq a, Eq b) => m a -> (a -> m b) -> m b ++instance (VecEq v, Num x) => EqMonad (v x) where+    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)+++data AsMonad m a where+   Embed  :: (EqMonad m, Eq a) => m a -> AsMonad m a+   Return :: EqMonad m => a -> AsMonad m a+   Bind   :: EqMonad m => AsMonad m a -> (a -> AsMonad m b) -> AsMonad m b+ +instance EqMonad m => Monad (AsMonad m) where+   return = Return+   (>>=) = Bind++unEmbed :: Eq a => AsMonad m a -> m a+unEmbed (Embed m) = m+unEmbed (Return a) = eqReturn a+unEmbed (Bind (Embed m) f) = m `eqBind` (unEmbed.f)+unEmbed (Bind (Return a) f) = unEmbed (f a)+unEmbed (Bind (Bind m f) g) = unEmbed (Bind m (\x -> Bind (f x) g))++++++++    ++
+ Setup.hs view
@@ -0,0 +1,3 @@+import Distribution.Simple+main = defaultMain+