-- | This module defines the QIO computation that represents Shor's factorisation -- algorithm. It makes use of the Arithmetic library, and the Quantum Fourier -- Transform. 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 -- | The inverse Quantum Fourier Transform is defined by reversing the QFT qftI :: QInt -> U qftI (QInt i) = urev (qft i) -- | The Hadamard transform can be mapped over the qubits in a Quantum Integer. hadamardsI :: QInt -> U hadamardsI (QInt xs) = hadamards xs -- | 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 -- | 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) measQ i0 -- | A classical (inefficient) implementation of the period finding subroutine period :: Int -> Int -> Int period m q = r where (_,r) = reduce (m,q) -- | A wrapper for Shor's algorithm, that returns prime factors of \n\. 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 -- | 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 -- | 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) -- | Calls the 'factorV'', and times the overall factorisation. 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")) -- | This function defines a quantum computation that returns a random index, that -- is used to pick from a list of integers that are co-prime to \n\. 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] -- | A different implementation of "rand_coprime", that defines a quantum -- 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 --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. -- | Integer division by 2. half :: Int -> Int half x = floor (fromIntegral x/2.0) -- | Reduces a pair of integers, by dividing them by thier gcd, -- until their gcd is 1. 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