diff --git a/ChangeLog b/ChangeLog --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,17 @@ ChangeLog +v0.2 2014/03/12 + (2014/03/12) NJR, PS1 - Added the new gridsynth algorithm from + N. J. Ross and P. Selinger, "Optimal ancilla-free Clifford+T + approximation of z-rotations", arXiv:1403.2975. Added a module + GridProblem for solving one- and two-dimensional grid equations, + and a module Diophantine for solving a Diophantine equation. + Removed the now obsolete Newsynth algorithm, and replaced it by a + backward compatibile interface to the new algorithm. New modules + GridSynth, QuadraticEquation, and StepComp. Additions and minor + improvements to EuclideanDomain, LaTeX, Matrix, and Ring. New + executable gridsynth. + v0.1.1.0 2014/02/05 (2014/02/04) PS1 - new functions euclid_divides and euclid_associates. (2014/02/04) PS1 - changed 'lobit' and improved its asymptotic diff --git a/Quantum/Synthesis/Diophantine.hs b/Quantum/Synthesis/Diophantine.hs new file mode 100644 --- /dev/null +++ b/Quantum/Synthesis/Diophantine.hs @@ -0,0 +1,432 @@ +-- | This module provides some number-theoretic functions, +-- particularly functions for solving the Diophantine equation +-- +-- \[center /t/[sup †]/t/ = ξ,] +-- +-- where ξ ∈ ℤ[√2] and /t/ ∈ ℤ[ω], or ξ ∈ [bold D][√2] and /t/ ∈ [bold D][ω]. +-- +-- In general, solving this equation can be hard, as it depends on the +-- ability to factor the integer /n/ = ξ[sup •]ξ into primes. We +-- formulate the solution as a step computation (see +-- "Quantum.Synthesis.StepComp"), so that the caller can dynamically +-- determine how much time to spend on solving the equation, or can +-- attempt to solve several such equations in parallel. +-- +-- In many cases, even a partial factorization of /n/ is sufficient to +-- determine that no solution exists. This implementation is written +-- to take advantage of such cases. + +module Quantum.Synthesis.Diophantine ( + -- * Diophantine solvers + diophantine, + diophantine_dyadic, + diophantine_associate, + + -- * Factoring + find_factor, + relatively_prime_factors, + + -- * Computations in ℤ[sub /n/] + power_mod, + root_of_negative_one, + root_mod, + ) where + +import Quantum.Synthesis.StepComp +import Quantum.Synthesis.EuclideanDomain +import Quantum.Synthesis.Ring + +import System.Random +import Control.Exception + +-- ---------------------------------------------------------------------- +-- * Diophantine solvers + +-- | Given ξ ∈ ℤ[√2], find /t/ ∈ ℤ[ω] such that /t/[sup †]/t/ = ξ, if +-- such /t/ exists, or return 'Nothing' otherwise. +diophantine :: (RandomGen g) => g -> ZRootTwo -> StepComp (Maybe ZOmega) +diophantine g xi + | xi == 0 = return (Just 0) + | xi < 0 = return Nothing + | adj2 xi < 0 = return Nothing + | otherwise = do + t <- diophantine_associate g xi + case t of + Nothing -> return Nothing + Just t -> do + let xi_associate = zroottwo_of_zomega (adj t * t) + let u = euclid_div xi xi_associate + case zroottwo_root u of + Nothing -> return Nothing + Just v -> return (Just (fromZRootTwo v * t)) + +-- | Given an element ξ ∈ [bold D][√2], find /t/ ∈ [bold D][ω] such +-- that /t/[sup †]/t/ = ξ, if such /t/ exists, or return 'Nothing' +-- otherwise. + +-- Implementation note: In the reduction from [bold D][√2] to ℤ[√2], +-- we can multiply by a power of λ√2 instead of √2. This has the same +-- effect of reducing the denominator exponent to 0 (note that λ is a +-- unit of the ring ℤ[√2]), but has the additional advantage that λ√2 +-- (unlike √2) is doubly positive, thereby preserving solvability of +-- the Diophantine equation. +-- +-- Similarly, in translating the solution back from ℤ[ω] to +-- [bold D][ω], we use the fact that 1\/(λ√2) = /u/[sup †]/u/, where +-- /u/ = (ω - /i/)\/√2. Also note that /u/ = δ⁻¹, where δ = 1 + ω. +diophantine_dyadic :: (RandomGen g) => g -> DRootTwo -> StepComp (Maybe DOmega) +diophantine_dyadic g xi = do + let k = denomexp xi + let (k',k'') = k `divMod` 2 + let xi' = to_whole ((lambda * roottwo)^k'' * 2^k' * xi) + t' <- diophantine g xi' + case t' of + Nothing -> return Nothing + Just t' -> return (Just (u^k'' * roothalf^k' * from_whole t')) + where + u = roothalf * (omega - i) + lambda = 1 + roottwo + +-- | Given ξ ∈ ℤ[√2], find /t/ ∈ ℤ[ω] such that /t/[sup †]/t/ ~ ξ, if +-- such /t/ exists, or 'Nothing' otherwise. Unlike 'diophantine', the +-- equation is only solved up to associates, i.e., up to a unit of the +-- ring. +diophantine_associate :: (RandomGen g) => g -> ZRootTwo -> StepComp (Maybe ZOmega) +diophantine_associate g xi + | xi == 0 = return (Just 0) + | otherwise = do + let d = euclid_gcd xi (adj2 xi) + let xi' = euclid_div xi d + res <- parallel_maybe (dioph_zroottwo_selfassociate g1 d) (dioph_zroottwo_assoc g2 xi') + case res of + Nothing -> return Nothing + Just (t1, t2) -> return (Just (t1*t2)) + where + (g1, g2) = split g + +-- ---------------------------------------------------------------------- +-- * Factoring + +-- | Given a positive composite integer /n/, find a non-trivial factor +-- of /n/ using a simple Pollard-rho method. The expected runtime is +-- O(√/p/), where /p/ is the size of the smallest prime factor. If +-- /n/ is not composite (i.e., if /n/ is prime or 1), this function +-- diverges. +find_factor :: (RandomGen g) => g -> Integer -> StepComp Integer +find_factor g n + | even n && n > 2 = return 2 + | otherwise = tick >> aux 2 (f 2) + where + (a, g2) = randomR (1, n-1) g + f x = (x^2 + a) `mod` n + aux x y + | d == 1 = tick >> aux (f x) (f (f y)) + | d == n = find_factor g2 n + | otherwise = return d + where + d = gcd (x-y) n + +-- | Given a factorization /n/ = /ab/ of some element of a Euclidean domain, find a factorization of /n/ into relatively prime factors, +-- +-- \[center /n/ = /u/ /c/[sub 1][sup /k/[sub 1]] ⋅ … ⋅ /c/[sub /m/][sup /k/[sub /m/]],] +-- +-- where /m/ ≥ 2, /u/ is a unit, and /c/[sub 1], …, /c/[sub /m/] are +-- pairwise relatively prime. +-- +-- While this is not quite a prime factorization of /n/, it can be a +-- useful intermediate step for computations that proceed by recursion +-- on relatively prime factors (such as Euler's φ-function, the +-- solution of Diophantine equations, etc.). +relatively_prime_factors :: (EuclideanDomain a) => a -> a -> (a, [(a, Integer)]) +relatively_prime_factors a b = aux 1 [a,b] [] where + aux u [] fs = (u, fs) + aux u (h:t) fs + | is_unit h = aux (h*u) t fs + aux u (h:t) fs = aux (u'*u) (hs ++ t) fs' where + (u', hs, fs') = aux2 h fs + + aux2 h [] = (1, [], [(h,1)]) + aux2 h ((f,k) : fs) + | euclid_associates h f = (u', [], (f,k+1) : fs) + | is_unit d = (u, hs, (f,k) : fs') + | otherwise = (1, [h `euclid_div` d, d] ++ replicate (fromInteger k) (f `euclid_div` d) ++ replicate (fromInteger k) d, fs) + where + d = euclid_gcd h f + (u, hs, fs') = aux2 h fs + u' = h `euclid_div` f + +-- ---------------------------------------------------------------------- +-- * Computations in ℤ[sub /n/] + +-- | Modular exponentiation, using the method of repeated squaring. +-- 'power_mod' /a/ /k/ /n/ computes /a/[sup /k/] (mod /n/). +power_mod :: Integer -> Integer -> Integer -> Integer +power_mod a k n + | k == 0 = 1 + | k == 1 = a `mod` n + | even k = (b*b) `mod` n + | otherwise = (b*b*a) `mod` n + where + b = power_mod a (k `div` 2) n + +-- | Compute a root of −1 in ℤ[sub /n/], where /n/ > 0. If /n/ is a +-- positive prime satisfying /n/ ≡ 1 (mod 4), this succeeds within an +-- expected number of 2 ticks. Otherwise, it probably diverges. +-- +-- As a special case, if this function notices that /n/ is not prime, +-- then it diverges without doing any additional work. +root_of_negative_one :: (RandomGen g) => g -> Integer -> StepComp Integer +root_of_negative_one g n = do + tick + let (b, g') = randomR (1, n-1) g + let h = power_mod b ((n-1) `div` 4) n + let r = (h*h) `mod` n + if r == n-1 then + return h + else + if r /= 1 then + diverge + else + root_of_negative_one g' n + +-- | Compute a root of /a/ in ℤ[sub /n/], where /n/ > 0. If /n/ is an +-- odd prime and /a/ is a non-zero square in ℤ[sub /n/], then this +-- succeeds in an expected number of 2 ticks. Otherwise, it probably +-- diverges. +root_mod :: (RandomGen g) => g -> Integer -> Integer -> StepComp Integer +root_mod g n a + | a `mod` n == -1 -- handle this special case more efficiently + = root_of_negative_one g n + | otherwise = tick >> res + where + (b, g') = randomR (0, n-1) g + (r,s) = (2*b `mod` n, b^2-a `mod` n) + (c,d) = pow (1,0) ((n-1) `div` 2) + res = case inv_mod n c of + Just c' + | (t1^2 - a) `mod` n == 0 -> return t1 + where + t = (1-d) * c' + t1 = (t+b) `mod` n + _ -> root_mod g' n a + + -- | 'mul' performs a multiplication in the ring + -- ℤ[sub /n/][t]\/(/t/²+/rt/+/s/). The elements /at/+/b/ are + -- represented as pairs (/a/,/b/). + mul :: (Integer, Integer) -> (Integer, Integer) -> (Integer, Integer) + mul (a,b) (c,d) = (a'',b'') where + (x,y,z) = (a*c, a*d+b*c, b*d) -- multiply polynomials + (a',b') = (y - x*r,z - x*s) -- reduce modulo /t/²+/rt/+/s/ + (a'',b'') = (a' `mod` n, b' `mod` n) -- reduce modulo /n/ + + -- | 'pow' takes a power in the ring + -- ℤ[sub /n/][t]\/(/t/²+/rt/+/s/. The elements /at/+/b/ are + -- represented as pairs (/a/,/b/). + pow :: (Integer, Integer) -> Integer -> (Integer, Integer) + pow x m + | m <= 0 = (0,1) + | odd m = x `mul` (x `pow` (m-1)) + | otherwise = y `mul` y where y = x `pow` (m `div` 2) + +-- ---------------------------------------------------------------------- +-- * Implementation details + +-- $ Our implementation of the top-level Diophantine equation solvers +-- proceeds through a series of special cases. The following functions +-- handle the special cases, and are not of independent interest. + +-- ---------------------------------------------------------------------- +-- ** Case: ξ is an integer + +-- | Given an integer /n/ ∈ ℤ, attempt to find /t/ ∈ ℤ[ω] such that +-- /t/[sup †]/t/ ~ /n/, or return 'Nothing' if no such /t/ exists. +-- +-- This function is optimized for the case when /n/ is prime, and +-- succeeds in an expected number of 2 ticks in this case. If /n/ is +-- not prime, this function probably diverges. +dioph_int_assoc_prime :: (RandomGen g) => g -> Integer -> StepComp (Maybe ZOmega) +dioph_int_assoc_prime g n + | n < 0 = dioph_int_assoc_prime g (-n) + | n == 0 = return (Just 0) + | n == 2 = return (Just roottwo) + | n_mod_4 == 1 = do + h <- root_of_negative_one g n + let t = euclid_gcd (fromInteger h+i) (fromInteger n) :: ZOmega + assert (adj t * t == fromInteger n) $ return (Just t) + | n_mod_8 == 3 = do + h <- root_mod g n (-2) + let t = euclid_gcd (fromInteger h+i*roottwo) (fromInteger n) :: ZOmega + assert (adj t * t == fromInteger n) $ return (Just t) + | n_mod_8 == 7 = do + h <- root_mod g n 2 + -- if n is prime, then 2 is a square. Conversely, if 2 is a + -- square, even if n is not prime, it implies that the Diophantine + -- equation has no solution. Because in this case, 2 is a square + -- for every prime divisor of n, so each such divisor must be + -- congruent to 1 or 7 (mod 8), so there must be at least one + -- prime divisor that occurs as an odd power and is congruent to 7 + -- mod n. + return Nothing + where + n_mod_4 = n `mod` 4 + n_mod_8 = n `mod` 8 + +-- | Given an integer /n/ ∈ ℤ, find /t/ ∈ ℤ[ω] such that /t/[sup †]/t/ +-- ~ /n/, if such /t/ exists, or return 'Nothing' if no such /t/ +-- exists. +-- +-- This function alternately calls 'dioph_int_assoc_prime' and +-- attempts to factor /n/. Therefore, it will eventually succeed; +-- however, the runtime depends on how hard it is to factor ξ. +dioph_int_assoc :: (RandomGen g) => g -> Integer -> StepComp (Maybe ZOmega) +dioph_int_assoc g n + | n < 0 = dioph_int_assoc g (-n) + | n == 0 = return (Just 0) + | n == 1 = return (Just 1) + | otherwise = interleave prime_solver factor_solver where + interleave p f = do + p <- subtask 4 p + case p of + Done res -> return res + _ -> do + f <- subtask 1000 f + case f of + Done (a, k) -> do + let b = n `div` a + let (u, facs) = relatively_prime_factors a b + forward (k `div` 2) $ dioph_int_assoc_powers g3 facs + _ -> interleave p f + + (g1, g') = split g + (g2, g3) = split g' + prime_solver = dioph_int_assoc_prime g1 n + factor_solver = with_counter $ speedup 30 $ find_factor g2 n + +-- | Given a factorization /n/ = /q/[sub 1][sup /k/[sub 1]]⋅…⋅/q/[sub +-- /m/][sup /k/[sub /m/]] of an integer /n/, where /q/[sub 1], …, +-- /q/[sub /m/] are pairwise relatively prime, find /t/ ∈ ℤ[ω] such +-- that /t/[sup †]/t/ ~ /n/, if such /t/ exists, or return 'Nothing' +-- if no such /t/ exists. +dioph_int_assoc_powers :: (RandomGen g) => g -> [(Integer, Integer)] -> StepComp (Maybe ZOmega) +dioph_int_assoc_powers g facs = do + res <- parallel_list_maybe [dioph_int_assoc_power g (n,k) | (n,k) <- facs] + case res of + Nothing -> return Nothing + Just sols -> return (Just (product sols)) + +-- | Given a pair of integers (/n/, /k/), find /t/ ∈ ℤ[ω] such that +-- /t/[sup †]/t/ ~ /n/[sup /k/], if such /t/ exists, or return +-- 'Nothing' if no such /t/ exists. +dioph_int_assoc_power :: (RandomGen g) => g -> (Integer, Integer) -> StepComp (Maybe ZOmega) +dioph_int_assoc_power g (n,k) + | even k = return (Just (fromInteger (n^(k `div` 2)))) + | otherwise = do + t <- dioph_int_assoc g n + case t of + Nothing -> return Nothing + Just t -> return (Just (t^k)) + +-- ---------------------------------------------------------------------- +-- ** Case: ξ ~ ξ[sup •] + +-- | Given ξ ∈ ℤ[√2] such that ξ ~ ξ[sup •], find /t/ ∈ ℤ[ω] such that +-- /t/[sup †]/t/ ~ ξ, if such /t/ exists, or return 'Nothing' if no +-- such /t/ exists. +dioph_zroottwo_selfassociate :: (RandomGen g) => g -> ZRootTwo -> StepComp (Maybe ZOmega) +dioph_zroottwo_selfassociate g xi + | xi == 0 = return (Just 0) + | otherwise = do + res <- dioph_int_assoc g n + case res of + Nothing -> return Nothing + Just t -> if euclid_divides roottwo r then + return (Just ((1+omega) * t)) + else + return (Just t) + where + RootTwo a b = xi + n = gcd a b + r = euclid_div xi (fromInteger n) + +-- ---------------------------------------------------------------------- +-- ** Case: gcd(ξ, ξ[sup •]) = 1 + +-- | Given ξ ∈ ℤ[√2] such that gcd(ξ, ξ[sup •]) = 1, attempt to find +-- /t/ ∈ ℤ[ω] such that /t/[sup †]/t/ ~ ξ, or return 'Nothing' if no +-- such /t/ exists. +-- +-- This function is optimized for the case when ξ is a prime in the +-- ring ℤ[√2]. In this case, it succeeds quickly, in an expected +-- number of 2 ticks. If ξ is not prime, this function probably +-- diverges. +dioph_zroottwo_assoc_prime :: (RandomGen g) => g -> ZRootTwo -> StepComp (Maybe ZOmega) +dioph_zroottwo_assoc_prime g xi + | xi == 0 = return (Just 0) + | n_mod_8 == 1 = do + h <- root_of_negative_one g n + let t = euclid_gcd (fromInteger h+i) (fromZRootTwo xi) :: ZOmega + assert ((adj t * t) `euclid_associates` fromZRootTwo xi) $ return (Just t) + | n_mod_8 == 7 = return Nothing + | otherwise = diverge + where + n_mod_8 = n `mod` 8 + n = abs (norm xi) + +-- | Given ξ ∈ ℤ[√2] such that gcd(ξ, ξ[sup •]) = 1, find /t/ ∈ ℤ[ω] +-- such that /t/[sup †]/t/ ~ ξ, if such /t/ exists, or return +-- 'Nothing' if no such /t/ exists. +-- +-- This function alternately calls 'dioph_int_assoc_prime' and +-- attempts to factor ξ. Therefore, it will eventually succeed. +-- However, the runtime depends on how hard it is to factor ξ. +dioph_zroottwo_assoc :: (RandomGen g) => g -> ZRootTwo -> StepComp (Maybe ZOmega) +dioph_zroottwo_assoc g xi + | xi == 0 = return (Just 0) + | otherwise = interleave prime_solver factor_solver where + interleave p f = do + p <- subtask 4 p + case p of + Done res -> return res + _ -> do + f <- subtask 1000 f + case f of + Done (a, k) -> do + let alpha = euclid_gcd xi (fromInteger a) + let beta = xi `euclid_div` alpha + let (u, facs) = relatively_prime_factors alpha beta + forward (k `div` 2) $ dioph_zroottwo_assoc_powers g3 facs + _ -> interleave p f + + (g1, g') = split g + (g2, g3) = split g' + prime_solver = dioph_zroottwo_assoc_prime g1 xi + factor_solver = with_counter $ speedup 30 $ find_factor g2 n + + n = abs (norm xi) + +-- | Given a factorization ξ = /q/[sub 1][sup /k/[sub 1]]⋅…⋅/q/[sub +-- /m/][sup /k/[sub /m/]] of some ξ ∈ ℤ[√2], where /q/[sub 1], …, +-- /q/[sub /m/] are pairwise relatively prime, find /t/ ∈ ℤ[ω] such +-- that /t/[sup †]/t/ ~ ξ, if such /t/ exists, or return 'Nothing' +-- if it can be proven not to exist. +dioph_zroottwo_assoc_powers :: (RandomGen g) => g -> [(ZRootTwo, Integer)] -> StepComp (Maybe ZOmega) +dioph_zroottwo_assoc_powers g facs = do + res <- parallel_list_maybe [dioph_zroottwo_assoc_power g (q,k) | (q,k) <- facs] + case res of + Nothing -> return Nothing + Just sols -> return (Just (product sols)) + +-- | Given a pair (ξ, /k/), with ξ ∈ ℤ[√2] and /k/ ≥ 0, find /t/ ∈ +-- ℤ[ω] such that /t/[sup †]/t/ ~ ξ[sup /k/], if such /t/ exists, or +-- return 'Nothing' if no such /t/ exists. +dioph_zroottwo_assoc_power :: (RandomGen g) => g -> (ZRootTwo, Integer) -> StepComp (Maybe ZOmega) +dioph_zroottwo_assoc_power g (xi,k) + | even k = return (Just (fromZRootTwo (xi^(k `div` 2)))) + | otherwise = do + t <- dioph_zroottwo_assoc g xi + case t of + Nothing -> return Nothing + Just t -> return (Just (t^k)) + + + diff --git a/Quantum/Synthesis/EuclideanDomain.hs b/Quantum/Synthesis/EuclideanDomain.hs --- a/Quantum/Synthesis/EuclideanDomain.hs +++ b/Quantum/Synthesis/EuclideanDomain.hs @@ -150,6 +150,20 @@ euclid_associates :: (EuclideanDomain a) => a -> a -> Bool euclid_associates a b = (a `euclid_divides` b) && (b `euclid_divides` a) +-- | Given elements /x/ and /y/ of a Euclidean domain, find the +-- largest /k/ such that /x/ can be written as /y/[sup /k/]/z/. +-- Return the pair (/k/, /z/). +-- +-- If /x/=0 or /y/ is a unit, return (/0/, /x/). +euclid_extract_power :: EuclideanDomain a => a -> a -> (Integer, a) +euclid_extract_power x y + | x == 0 = (0, x) + | is_unit y = (0, x) + | y `euclid_divides` x = (k+1, z) + | otherwise = (0, x) + where + (k, z) = euclid_extract_power (x `euclid_div` y) y + -- ---------------------------------------------------------------------- -- * Auxiliary functions diff --git a/Quantum/Synthesis/GridProblems.hs b/Quantum/Synthesis/GridProblems.hs new file mode 100644 --- /dev/null +++ b/Quantum/Synthesis/GridProblems.hs @@ -0,0 +1,931 @@ +-- | This module provides functions for solving one- and +-- two-dimensional grid problems. + +module Quantum.Synthesis.GridProblems where + +import Quantum.Synthesis.Ring +import Quantum.Synthesis.Matrix +import Quantum.Synthesis.QuadraticEquation + +import System.Random + +-- ---------------------------------------------------------------------- +-- * 1-dimensional grid problems + +-- $ The /1-dimensional grid problem/ is the following: given closed +-- intervals /A/ and /B/ of the real numbers, find all α ∈ ℤ[√2] such +-- that α ∈ /A/ and α[sup •] ∈ /B/. +-- +-- Let Δx be the size of /A/, and Δy the size of /B/. It is a theorem +-- that the 1-dimensional grid problem has at least one solution if +-- ΔxΔy ≥ (1 + √2)², and at most one solution if ΔxΔy < 1. +-- Asymptotically, the expected number of solutions is ΔxΔy/\√8. +-- +-- The following functions provide solutions to a number of variations +-- of the grid problem. All functions are formulated so that the +-- intervals can be specified exactly (using a type such as +-- 'QRootTwo'), or approximately (using a type such as 'Double' or +-- 'FixedPrec' /e/). + +-- ---------------------------------------------------------------------- +-- ** General solutions + +-- | Given two intervals /A/ = [/x/₀, /x/₁] and /B/ = [/y/₀, /y/₁] of +-- real numbers, output all solutions α ∈ ℤ[√2] of the 1-dimensional +-- grid problem for /A/ and /B/. The list is produced lazily, and is +-- sorted in order of increasing α. +gridpoints :: (RootTwoRing r, Fractional r, Floor r, Ord r) => (r, r) -> (r, r) -> [ZRootTwo] +gridpoints (x0, x1) (y0, y1) + | dy <= 0 && dx > 0 = + map adj2 $ gridpoints (y0, y1) (x0, x1) + | dy >= lambda && even n = + map (lambda_inv_n *) $ gridpoints (lambda_n*x0, lambda_n*x1) (lambda_bul_n*y0, lambda_bul_n*y1) + | dy >= lambda && odd n = + map (lambda_inv_n *) $ gridpoints (lambda_n*x0, lambda_n*x1) (lambda_bul_n*y1, lambda_bul_n*y0) + | dy > 0 && dy < 1 && even n = + map (lambda_m *) $ gridpoints (lambda_inv_m*x0, lambda_inv_m*x1) (lambda_bul_inv_m*y0, lambda_bul_inv_m*y1) + | dy > 0 && dy < 1 && odd n = + map (lambda_m *) $ gridpoints (lambda_inv_m*x0, lambda_inv_m*x1) (lambda_bul_inv_m*y1, lambda_bul_inv_m*y0) + | otherwise = + [ RootTwo a b | a <- [amin..amax], b <- [bmin a..bmax a], test a b ] + where + dx = x1 - x0 + dy = y1 - y0 + (n, _) = floorlog lambda dy + m = -n + + lambda_m = lambda^m + lambda_n = lambda^n + lambda_bul_n = (-lambda_inv)^n + lambda_inv_m = lambda_inv^m + lambda_bul_inv_m = (-lambda)^m + lambda_inv_n = lambda_inv^n + + within x (x0, x1) = x0 <= x && x <= x1 + amin = ceiling_of ((x0 + y0) / 2) + amax = floor_of ((x1 + y1) / 2) + bmin a = ceiling_of ((fromInteger a - y1) / roottwo) + bmax a = floor_of ((fromInteger a - y0) / roottwo) + test a b = fromZRootTwo x `within` (x0, x1) && fromZRootTwo (adj2 x) `within` (y0, y1) + where x = RootTwo a b + +-- | Like 'gridpoints', but only produce solutions /a/ + /b/√2 where +-- /a/ has the same parity as the given integer. +gridpoints_parity :: (RootTwoRing r, Fractional r, Floor r, Ord r) => Integer -> (r, r) -> (r, r) -> [ZRootTwo] +gridpoints_parity e (x0,x1) (y0,y1) = do + z' <- gridpoints (x0', x1') (-y1', -y0') + return (roottwo * z' + fromInteger e2) + where + x0' = (x0 - e') / roottwo + x1' = (x1 - e') / roottwo + y0' = (y0 - e') / roottwo + y1' = (y1 - e') / roottwo + e' = fromInteger e2 + e2 = e `mod` 2 + +-- ---------------------------------------------------------------------- +-- ** Randomized solutions + +-- | Given two intervals /A/ = [/x/₀, /x/₁] and /B/ = [/y/₀, /y/₁] of +-- real numbers, and a source of randomness, output a random solution +-- α ∈ ℤ[√2] of the 1-dimensional grid problem for /A/ and /B/. +-- +-- Note: the randomness is not uniform. To ensure that the set of +-- solutions is non-empty, we must have ΔxΔy ≥ (1 + √2)², where Δx = +-- /x/₁ − /x/₀ ≥ 0 and Δy = /y/₁ − /y/₀ ≥ 0. If there are no solutions +-- at all, the function returns 'Nothing'. +gridpoint_random :: (RootTwoRing r, Fractional r, Floor r, Ord r, RandomGen g) => (r, r) -> (r, r) -> g -> Maybe ZRootTwo +gridpoint_random (x0, x1) (y0, y1) g = z + where + dx = max 0 (x1 - x0) + dy = max 0 (y1 - y0) + area = dx * dy + n = floor_of (area + 1) + (i,_) = randomR (0, n-1) g + r = fromInteger i / fromInteger n + pts = gridpoints (x0 + r * dx, x1) (y0, y1) ++ gridpoints (x0, x1) (y0, y1) + z = case pts of + h:t -> Just h + [] -> Nothing + +-- | Like 'gridpoint_random', but only produce solutions /a/ + /b/√2 +-- where /a/ has the same parity as the given integer. +gridpoint_random_parity :: (RootTwoRing r, Fractional r, Floor r, Ord r, RandomGen g) => Integer -> (r, r) -> (r, r) -> g -> Maybe ZRootTwo +gridpoint_random_parity e (x0,x1) (y0,y1) g = do + z' <- gridpoint_random (x0', x1') (-y1', -y0') g + return (roottwo * z' + fromInteger e2) + where + x0' = (x0 - e') / roottwo + x1' = (x1 - e') / roottwo + y0' = (y0 - e') / roottwo + y1' = (y1 - e') / roottwo + e' = fromInteger e2 + e2 = e `mod` 2 + +-- ---------------------------------------------------------------------- +-- ** Scaled solutions + +-- $ The scaled version of the 1-dimensional grid problem is the +-- following: given closed intervals /A/ and /B/ of the real numbers, +-- and /k/ ≥ 0, find all α ∈ ℤ[√2] \/ √2[sup /k/] such that α ∈ /A/ +-- and α[sup •] ∈ /B/. + +-- | Given intervals /A/ = [/x/₀, /x/₁] and /B/ = [/y/₀, /y/₁], and an +-- integer /k/ ≥ 0, output all solutions α ∈ ℤ[√2] \/ √2[sup /k/] of +-- the scaled 1-dimensional grid problem for /A/, /B/, and /k/. The +-- list is produced lazily, and is sorted in order of increasing /α/. +gridpoints_scaled :: (RootTwoRing r, Fractional r, Floor r, Ord r) => (r, r) -> (r, r) -> Integer -> [DRootTwo] +gridpoints_scaled (x0, x1) (y0, y1) k = do + w <- gridpoints (x0', x1') (y0', y1') + return (scale * fromZRootTwo w) + where + scale = roothalf^k + scale_inv = roottwo^k + (x0', x1') = (scale_inv * x0, scale_inv * x1) + (y0', y1') + | even k = (scale_inv * y0, scale_inv * y1) + | otherwise = (-scale_inv * y1, -scale_inv * y0) + +-- | Like 'gridpoints_scaled', but assume /k/ ≥ 1, take an additional +-- parameter β ∈ ℤ[√2] \/ √2[sup /k/], and return only those α such +-- that β − α ∈ ℤ[√2] \/ √2[sup /k-1/]. +gridpoints_scaled_parity :: (RootHalfRing r, Fractional r, Floor r, Ord r) => DRootTwo -> (r, r) -> (r, r) -> Integer -> [DRootTwo] +gridpoints_scaled_parity beta (x0, x1) (y0, y1) k + | denomexp beta <= k-1 = gridpoints_scaled (x0, x1) (y0, y1) (k-1) + | otherwise = do + z' <- gridpoints_scaled (x0+offs', x1+offs') (y0+offs_bul', y1+offs_bul') (k-1) + return (z' - offs) + where + offs = roothalf^k + offs_bul = adj2 offs + offs' = fromDRootTwo offs + offs_bul' = fromDRootTwo offs_bul + +-- ---------------------------------------------------------------------- +-- * 2-dimensional grid problems + +-- $ The /2-dimensional grid problem/ is the following: given bounded +-- convex subsets /A/ and /B/ of ℂ with non-empty interior, find all +-- /u/ ∈ ℤ[ω] such that /u/ ∈ /A/ and /u/[sup •] ∈ /B/. + +-- ---------------------------------------------------------------------- +-- ** Representation of convex sets + +-- $ Since convex sets /A/ and /B/ are inputs of the 2-dimensional +-- grid problem, we need a way to specify convex subsets of ℂ. Our +-- specification of a convex sets consists of three parts: +-- +-- * a /bounding ellipse/ for the convex set; +-- +-- * a /characteristic function/, which tests whether any given point +-- is an element of the convex set; and +-- +-- * a /line intersector/, which estimates the intersection of any +-- given straight line and the convex set. + +-- | A point in the plane. +type Point r = (r,r) + +-- | Convert a point with coordinates in 'DRootTwo' to a point with +-- coordinates in any 'RootHalfRing'. +point_fromDRootTwo :: (RootHalfRing r) => Point DRootTwo -> Point r +point_fromDRootTwo (x, y) = (fromDRootTwo x, fromDRootTwo y) + +-- | An operator is a real 2×2-matrix. +type Operator a = Matrix Two Two a + +-- | An /ellipse/ is given by an operator /D/ and a center /p/; the +-- ellipse in this case is +-- +-- /A/ = { /v/ | (/v/-/p/)[sup †] /D/ (/v/-/p/) ≤ 1}. +data Ellipse r = Ellipse (Operator r) (Point r) + deriving (Show) + +-- | The /characteristic function/ of a set /A/ inputs a point /p/, +-- and outputs 'True' if /p/ ∈ /A/ and 'False' otherwise. +-- +-- The point /p/ is given of an exact type, so characteristic +-- functions have the opportunity to use infinite precision. +type CharFun = Point DRootTwo -> Bool + +-- | A /line intersector/ knows about some compact convex set +-- /A/. Given a straight line /L/, it computes an approximation of the +-- intersection of /L/ and /A/. +-- +-- More specifically, /L/ is given as a parametric equation /p/(/t/) = +-- /v/ + /tw/, where /v/ and /w/ ≠ 0 are vectors. Given /v/ and /w/, the +-- line intersector returns /t/₀ and /t/₁ such that /p/(/t/) ∈ /A/ implies +-- /t/ ∈ [/t/₀, /t/₁]. +-- +-- Line intersectors should overestimate (\"fatten\") the convex set +-- slightly, to guard against possible round-off errors. +type LineIntersector r = (Point DRootTwo -> Point DRootTwo -> (r, r)) + +-- | A compact convex set is given by a bounding ellipse, a +-- characteristic function, and a line intersector. +data ConvexSet r = ConvexSet (Ellipse r) CharFun (LineIntersector r) + +instance (Show r) => Show (ConvexSet r) where + show (ConvexSet ell tst int) = "ConvexSet (" ++ show ell ++ ", ..., ...)" + +-- ---------------------------------------------------------------------- +-- ** Specific convex sets + +-- | The closed unit disk. +unitdisk :: (Fractional r, Ord r, RootHalfRing r, Quadratic r) => ConvexSet r +unitdisk = ConvexSet ell tst int where + ell = Ellipse 1 (0,0) + + int p v + | q == Nothing = (1, 0) + | otherwise = (t0, t1) + where + a = iprod v v + b = 2 * iprod v p + c = iprod p p - 1 + q = quadratic (fromDRootTwo a) (fromDRootTwo b) (fromDRootTwo c) + Just (t0, t1) = q + + tst (x,y) = x^2 + y^2 <= 1 + +-- | A closed rectangle with the given dimensions. +rectangle :: (Fractional r, Ord r, RootHalfRing r) => (r,r) -> (r,r) -> ConvexSet r +rectangle (x0,x1) (y0,y1) = ConvexSet ell tst int where + w = x1-x0 + h = y1-y0 + center = ((x0+x1) / 2, (y0+y1) / 2) + mat = toOperator ((2/w^2,0), (0,2/h^2)) + ell = Ellipse mat center + tst (x, y) = (fromDRootTwo x `within` (x0, x1)) && (fromDRootTwo y `within` (y0, y1)) + int p v = int_internal (point_fromDRootTwo p) (point_fromDRootTwo v) + int_internal p v + | vx == 0 && px `within` (x0, x1) = (min t0y t1y, max t0y t1y) + | vx == 0 = (1, 0) + | vy == 0 && py `within` (y0, y1) = (min t0x t1x, max t0x t1x) + | vy == 0 = (1, 0) + | otherwise = (t0, t1) + where + (px, py) = p + (vx, vy) = v + t0x = (x0 - px) / vx + t1x = (x1 - px) / vx + t0y = (y0 - py) / vy + t1y = (y1 - py) / vy + t0 = max (min t0x t1x) (min t0y t1y) + t1 = min (max t0x t1x) (max t0y t1y) + +-- ---------------------------------------------------------------------- +-- ** General solutions + +-- | Given bounded convex sets /A/ and /B/, enumerate all solutions +-- /u/ ∈ ℤ[ω] of the 2-dimensional grid problem for /A/ and /B/. +gridpoints2 :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> [DOmega] +gridpoints2 setA setB = gridpoints2_scaled setA setB 0 + +-- ---------------------------------------------------------------------- +-- ** Scaled solutions + +-- $ The scaled version of the 2-dimensional grid problem is the +-- following: given bounded convex subsets /A/ and /B/ of ℂ with +-- non-empty interior, and /k/ ≥ 0, find all /u/ ∈ ℤ[ω] \/ √2[sup /k/] +-- such that /u/ ∈ /A/ and /u/[sup •] ∈ /B/. + +-- | Given bounded convex sets /A/ and /B/, return a function that can +-- input a /k/ and enumerate all solutions of the two-dimensional +-- scaled grid problem for /A/, /B/, and /k/. +-- +-- Note: a large amount of precomputation is done on the sets /A/ and +-- /B/, so it is beneficial to call this function only once for a +-- given pair of sets, and then possibly call the result many times +-- for different /k/. In other words, for optimal performance, the +-- function should be used like this: +-- +-- > let solver = gridpoints2_scaled setA setB +-- > let solutions0 = solver 0 +-- > let solutions1 = solver 1 +-- > ... +-- +-- Note: the gridpoints are computed in some deterministic (but +-- unspecified) order. They are not randomized. +gridpoints2_scaled :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> Integer -> [DOmega] +gridpoints2_scaled setA setB = solutions_fun + where + ConvexSet ellA tstA intA = setA + ConvexSet ellB tstB intB = setB + Ellipse matA ctrA = ellA + Ellipse matB ctrB = ellB + + -- Find the grid operator + opG = to_upright (matA, matB) + opG_inv = special_inverse opG + + -- Change the coordinate system + setA' = convex_transform opG_inv setA + setB' = convex_transform (adj2 opG_inv) setB + bboxA' = boundingbox setA' + bboxB' = boundingbox setB' + ConvexSet ellA' tstA' intA' = setA' + ConvexSet ellB' tstB' intB' = setB' + ((x0A, x1A), (y0A, y1A)) = bboxA' + ((x0B, x1B), (y0B, y1B)) = bboxB' + + solutions_fun k = do + + -- Enumerate the solutions in the y-coordinate + beta' <- gridpoints_scaled (fatten_interval (y0A, y1A)) (fatten_interval (y0B, y1B)) (k+1) + let beta'_bul = adj2 beta' + + let xs = gridpoints_scaled (x0A, x1A) (x0B, x1B) (k+1) + x0 <- take 1 xs + let x0_bul = adj2 x0 + let dx = roothalf^k + let dx_bul = adj2 dx + + -- Intersect that y-coordinate with the convex sets + let (t0A, t1A) = intA' (x0, beta') (dx, 0) + let (t0B, t1B) = intB' (x0_bul, beta'_bul) (dx_bul, 0) + + -- offsets for slightly fattening the intervals, in a way that + -- does not add more than a small constant number of candidates + -- alpha' on both sides of the interval. + let dtA = min 1 (10 / (2^k * (x1B - x0B))) + let dtB = min 1 (10 / (2^k * (x1A - x0A))) + + -- Enumerate the solutions in the x-coordinate (ensuring correct + -- parity to make sure it's a grid point) + -- + -- For parity, we need: + -- 1/√2^k | alpha' - beta' + -- <=> 1/√2^k | alpha'_offs * dx + x0 - beta' + -- <=> 1 | alpha'_offs + (x0 - beta') * √2^k + alpha'_offs <- gridpoints_scaled_parity ((beta'-x0)*roottwo^k) (t0A-dtA, t1A+dtA) (t0B-dtB, t1B+dtB) 1 + let alpha' = alpha'_offs * dx + x0 + + -- Convert back to the original coordinate system + let (alpha,beta) = point_transform opG (alpha',beta') + + case tstA (alpha,beta) && tstB (adj2 alpha, adj2 beta) of + True -> do + let z = fromDRootTwo alpha + i * fromDRootTwo beta :: DOmega + return z + False -> do + [] + +-- | Given bounded convex sets /A/ and /B/, enumerate all solutions of +-- the two-dimensional scaled grid problem for all /k/ ≥ 0. Each +-- solution is only enumerated once, and the solutions are enumerated +-- in order of increasing /k/. +gridpoints2_increasing :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> [DOmega] +gridpoints2_increasing setA setB = solutions + where + solutions_fun = gridpoints2_scaled setA setB + solutions = solutions_fun 0 ++ additional_solutions 1 + additional_solutions k = exact_solutions k ++ additional_solutions (k+1) + exact_solutions k = [ z | z <- solutions_fun k, denomexp z == k ] + +-- ---------------------------------------------------------------------- +-- * Implementation details + +-- $ Our solution of the 2-dimensional grid problem follows the paper +-- +-- * N. J. Ross and P. Selinger, \"Optimal ancilla-free Clifford+/T/ +-- approximation of /z/-rotations\". . + +-- ---------------------------------------------------------------------- +-- ** Positive operators and ellipses + +-- | Construct a 2×2-matrix, by rows. +toOperator :: ((a, a), (a, a)) -> Operator a +toOperator ((a, b), (c, d)) = matrix2x2 (a,b) (c,d) + +-- | Extract the entries of a 2×2-matrix, by rows. +fromOperator :: Operator a -> ((a, a), (a, a)) +fromOperator m = from_matrix2x2 m + +-- | Convert an operator with entries in 'DRootTwo' to an operator with +-- entries in any 'RootHalfRing'. +op_fromDRootTwo :: (RootHalfRing r) => Operator DRootTwo -> Operator r +op_fromDRootTwo m = matrix_map fromDRootTwo m + +-- | The (/b/,/z/)-representation of a positive operator with determinant 1 is +-- +-- \[image bz.png] +-- +-- where /b/, /z/ ∈ ℝ and /e/ > 0 with /e/² = /b/² + 1. Create such an +-- operator from parameters /b/ and /z/. +operator_from_bz :: (RootTwoRing a, Floating a) => a -> a -> Operator a +operator_from_bz b z = toOperator ((a, b), (b, d)) where + lambda_z = lambda**z + a = e/lambda_z + d = e*lambda_z + e = sqrt(1 + b^2) + +-- | Conversely, given a positive definite real operator of +-- determinant 1, return the parameters (/b/, /z/). This is the +-- inverse of 'operator_from_bz'. For efficiency reasons, the +-- parameter /z/, which is a logarithm, is modeled as a 'Double'. +operator_to_bz :: (Fractional a, Real a, RootTwoRing a) => Operator a -> (a, Double) +operator_to_bz m = (b, z) where + ((a, b), (c, d)) = fromOperator m + lambda_z_squared = d / a + z = 0.5 * logBase_double lambda lambda_z_squared + +-- | A version of 'operator_to_bz' that returns (/b/, λ[sup 2/z/]) +-- instead of (/b/, /z/). +-- This is a critical optimization, as this function is called often, +-- and logarithms are relatively expensive to compute. +operator_to_bl2z :: (Floating a, Real a) => Operator a -> (a, a) +operator_to_bl2z m = (b, l2z) where + ((a, b), (c, d)) = fromOperator m + l2z = d / a + +-- | The determinant of a 2×2-matrix. +det :: (Ring a) => Operator a -> a +det m = a*d - b*c where + ((a, b), (c, d)) = fromOperator m + +-- | Compute the skew of a positive operator of determinant 1. We +-- define the /skew/ of a positive definite real operator /D/ to be +-- +-- \[image skew.png] +operator_skew :: (Ring a) => Operator a -> a +operator_skew m = b*c where + ((a, b), (c, d)) = fromOperator m + +-- | Compute the uprightness of a positive operator /D/. +-- +-- The /uprightness/ of /D/ is the ratio of the area of the ellipse +-- /E/ = {/v/ | /v/[sup †]/D//v/ ≤ 1} to the area of its bounding box +-- /R/. It is given by +-- +-- \[image area.png] +-- +-- \[image ellipse-rectangle.png] +uprightness :: (Floating a) => Operator a -> a +uprightness m = pi/4 * sqrt(det m / (a*d)) + where + ((a, b), (c, d)) = fromOperator m + +-- ---------------------------------------------------------------------- +-- ** States + +-- | A state is a pair (/D/, Δ) of real positive definite matrices of +-- determinant 1. It encodes a pair of ellipses. +type OperatorPair a = (Operator a, Operator a) + +-- | The /skew/ of a state is the sum of the skews of the two +-- operators. +skew :: (Ring a) => OperatorPair a -> a +skew (m1,m2) = operator_skew m1 + operator_skew m2 + +-- | The /bias/ of a state is ζ - /z/. +bias :: (Fractional a, Real a, RootTwoRing a) => OperatorPair a -> Double +bias (matA,matB) = (zeta - z) where + (b, z) = operator_to_bz matA + (beta, zeta) = operator_to_bz matB + +-- ---------------------------------------------------------------------- +-- ** Grid operators + +-- $ Consider the set ℤ[ω] ⊆ ℂ. In identifying ℂ with ℝ², we can +-- alternatively identify ℤ[ω] with the set of all vectors (/x/, +-- /y/)[sup †] ∈ ℝ² of the form +-- +-- * /x/ = /a/ + /a/'\/√2, +-- +-- * /y/ = /b/ + /b/'\/√2, +-- +-- such that /a/, /a/', /b/, /b/' ∈ ℤ and /a/' ≡ /b/' (mod 2). +-- +-- A /grid operator/ is a linear operator /G/ : ℝ² → ℝ² such that /G/ +-- maps ℤ[ω] to itself. We can characterize the grid operators as +-- the operators of the form +-- +-- \[image gridop.png] +-- +-- such that: +-- +-- * /a/ + /b/ + /c/ + /d/ ≡ 0 (mod 2) and +-- +-- * /a/' ≡ /b/' ≡ /c/' ≡ /d/' (mod 2). +-- +-- A /special grid operator/ is a grid operator of determinant ±1. All +-- special grid operators are invertible, and the inverse is again a +-- special grid operator. +-- +-- Since all coordinates of ℤ[ω] (as a subset of ℝ²), and all entries +-- of grid operators, can be represented as elements of the ring [bold +-- D][√2], the automorphism /x/ ↦ /x/[sup •], which maps /a/ + /b/√2 +-- to /a/ - /b/√2 (for rational /a/ and /b/), is well-defined for +-- them. +-- +-- In this section, we define some particular special grid operators +-- that are used in the Step Lemma. + +-- | The special grid operator /R/: a clockwise rotation by 45°. +-- +-- \[image gridop-R.png] +opR :: (RootHalfRing r) => Operator r +opR = roothalf * toOperator ((1, -1), (1, 1)) + +-- | The special grid operator /A/: a clockwise shearing with offset +-- 2, parallel to the /x/-axis. +-- +-- \[image gridop-A.png] +opA :: (Ring r) => Operator r +opA = matrix2x2 (1, -2) (0, 1) + +-- | The special grid operator /A/⁻¹: a counterclockwise shearing with offset +-- 2, parallel to the /x/-axis. +-- +-- \[image gridop-Ai.png] +opA_inv :: (Ring r) => Operator r +opA_inv = matrix2x2 (1, 2) (0, 1) + +-- | The operator /A/[sup /k/]. +opA_power :: (RootTwoRing r) => Integer -> Operator r +opA_power k + | k >= 0 = opA^k + | otherwise = opA_inv^(-k) + +-- | The special grid operator /B/: a clockwise shearing with offset +-- √2, parallel to the /x/-axis. +-- +-- \[image gridop-B.png] +opB :: (RootTwoRing r) => Operator r +opB = matrix2x2 (1, roottwo) (0, 1) + +-- | The special grid operator /B/⁻¹: a counterclockwise shearing with offset +-- √2, parallel to the /x/-axis. +-- +-- \[image gridop-Bi.png] +opB_inv :: (RootTwoRing r) => Operator r +opB_inv = matrix2x2 (1, -roottwo) (0, 1) + +-- | The operator /B/[sup /k/]. +opB_power :: (RootTwoRing r) => Integer -> Operator r +opB_power k + | k >= 0 = opB^k + | otherwise = opB_inv^(-k) + +-- | The special grid operator /K/. +-- +-- \[image gridop-K.png] +opK :: (RootHalfRing r) => Operator r +opK = roothalf * matrix2x2 (-lambda_inv, -1) (lambda, 1) + +-- | The Pauli /X/ operator is a special grid operator. +-- +-- \[image gridop-X.png] +opX :: (Ring r) => Operator r +opX = matrix2x2 (0, 1) (1, 0) + +-- | The Pauli operator /Z/ is a special grid operator. +-- +-- \[image gridop-Z.png] +opZ :: (Ring r) => Operator r +opZ = matrix2x2 (1, 0) (0, -1) + +-- | The special grid operator /S/: a scaling by λ = 1+√2 in the +-- /x/-direction, and by λ⁻¹ = -1+√2 in the /y/-direction. +-- +-- \[image gridop-S.png] +-- +-- The operator /S/ is not used in the paper, but we use it here for +-- a more efficient implementation of large shifts. The point is that +-- /S/ is a grid operator, but shifts in increments of 4, whereas the +-- Shift Lemma uses non-grid operators but shifts in increments of 2. +opS :: (RootTwoRing r) => Operator r +opS = toOperator((lambda, 0), (0, lambda_inv)) + +-- | The special grid operator /S/⁻¹, the inverse of 'opS'. +-- +-- \[image gridop-Si.png] +opS_inv :: (RootTwoRing r) => Operator r +opS_inv = matrix2x2 (lambda_inv, 0) (0, lambda) + +-- | Return /S/[sup /k/]. +opS_power :: (RootTwoRing r) => Integer -> Operator r +opS_power k + | k >= 0 = opS^k + | otherwise = opS_inv^(-k) + +-- ---------------------------------------------------------------------- +-- ** Action of grid operators on states + +-- | Compute the right action of a grid operator /G/ on a state (/D/, +-- Δ). This is defined as: +-- +-- (/D/, Δ) ⋅ /G/ := (/G/[sup †]/D//G/, /G/[sup •T]Δ/G/[sup •]). +action :: (RealFrac r, RootHalfRing r, Adjoint r) => (Operator r, Operator r) -> Operator DRootTwo -> (Operator r, Operator r) +action (a,b) g = (g1 * a * g2, g3 * b * g4) where + g1 = adj g2 + g2 = op_fromDRootTwo g + g3 = adj g4 + g4 = op_fromDRootTwo (adj2 g) + +-- ---------------------------------------------------------------------- +-- ** Shifts + +-- $ A shift is not quite the application of a grid operator, because +-- the shifts σ and τ actually involve a square root of λ. However, +-- they can be used to define an operation on states. + +-- | Given an operator /D/, compute σ[sup /k/]/D/σ[sup /k/]. +shift_sigma :: (RootTwoRing a) => Integer -> Operator a -> Operator a +shift_sigma k m = matrix2x2 (lambdapower k * a, b) (c, lambdapower (-k) * d) where + ((a,b),(c,d)) = fromOperator m + +-- | Given an operator Δ, compute τ[sup /k/]Δτ[sup /k/]. +shift_tau :: (RootTwoRing a) => Integer -> Operator a -> Operator a +shift_tau k m = matrix2x2 (lambdapower (-k) * a, signpower k * b) (c * signpower k, lambdapower k * d) where + ((a,b),(c,d)) = fromOperator m + +-- | Compute the /k/-shift of a state (/D/,Δ). +shift_state :: (RootTwoRing a) => Integer -> OperatorPair a -> OperatorPair a +shift_state k (d,delta) = (shift_sigma k d, shift_tau k delta) + +-- ---------------------------------------------------------------------- +-- ** Skew reduction + +-- | An implementation of the /A/-Lemma. Given /z/ and ζ, compute the +-- integer /m/ such that the operator /A/[sup /m/] reduces the skew. +lemma_A :: (RealFrac r, RootTwoRing r, Floating r) => r -> r -> Integer +lemma_A z zeta = n where + n = max 1 (floor (lambda ** c / 2)) + c = min z zeta + +-- | An implementation of the /B/-Lemma. Given /z/ and ζ, compute the +-- integer /m/ such that the operator /B/[sup /m/] reduces the skew. +lemma_B :: (RealFrac r, RootTwoRing r, Floating r) => r -> r -> Integer +lemma_B z zeta = n where + n = max 1 (floor (lambda ** c / roottwo)) + c = min z zeta + +-- | A version of 'lemma_A' that inputs λ[sup 2/z/] instead of /z/ and +-- λ[sup 2ζ] instead of ζ. Compute the constant /m/ such that the +-- operator /A/[sup /m/] reduces the skew. +lemma_A_l2 :: (RealFrac r, RootTwoRing r, Floating r) => r -> r -> Integer +lemma_A_l2 l2z l2zeta = n where + n = max 1 (intsqrt (floor (l2c / 4))) + l2c = min l2z l2zeta + +-- | A version of 'lemma_B' that inputs λ[sup 2/z/] instead of /z/ and +-- λ[sup 2ζ] instead of ζ. Compute the constant /m/ such that the +-- operator /B/[sup /m/] reduces the skew. +lemma_B_l2 :: (RealFrac r, RootTwoRing r, Floating r) => r -> r -> Integer +lemma_B_l2 l2z l2zeta = n where + n = max 1 (intsqrt (floor (l2c / 2))) + l2c = min l2z l2zeta + +-- | An implementation of the Step Lemma. Input a state (/D/,Δ). If +-- the skew is > 15, produce a special grid operator whose action +-- reduces Skew(/D/,Δ) by at least 5%. If the skew is ≤ 15 and β ≥ 0 +-- and z + ζ ≥ 0, do nothing. Otherwise, produce a special grid +-- operator that ensures β ≥ 0 and z + ζ ≥ 0. +step_lemma :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r) => OperatorPair r -> Maybe (Operator DRootTwo) +step_lemma (matA,matB) + -- First ensure that β ≥ 0, by applying /Z/ if necessary. + | beta < 0 + = wlog_using opZ + + -- Then ensure that z + ζ ≥ 0, by applying /X/ if necessary. + -- Case: z + zeta < 0. + | l2z * l2zeta < 1 + = wlog_using opX + + -- If the bias is greater than 2, use the grid operator /S/. This is + -- more efficient than applying the Shift Lemma. + -- Todo: ensure numeric stability + -- Case: |z-ζ| > 2. + | l2z_minus_zeta > 33.971 || l2z_minus_zeta < 0.029437 + = wlog_using (opS_power (round (logLambda l2z_minus_zeta / 8))) + + -- If the skew is below threshold, stop. + | skew (matA,matB) <= 15 + = Nothing + + -- If the bias is greater than 1, apply a shift. + -- Todo: ensure numeric stability + -- Case: |z-ζ| > 1. + | l2z_minus_zeta > 5.8285 || l2z_minus_zeta < 0.17157 + = with_shift (round (logLambda l2z_minus_zeta / 4)) + + -- Cases 1.1 and 2.1: z ∈ [-0.8, 0.8] and ζ ∈ [-0.8, 0.8]. + -- Region R. + | l2z `within` (0.24410, 4.0968) && l2zeta `within` (0.24410, 4.0968) + = Just opR + + -- Case 1.2: b ≥ 0 and z ≤ 0.3 and ζ ≥ 0.8. + -- Region K. + | b >= 0 && l2z <= 1.6969 + = Just opK + + -- Case 1.4: b ≥ 0 and z ≥ 0.8 and ζ ≤ 0.3. + -- Region K•. + | b >= 0 && l2zeta <= 1.6969 + = Just (adj2 opK) + + -- Case 1.6: b ≥ 0 and z ≥ 0.3 and zeta ≥ 0.3. + -- Region A^m. + | b >= 0 + = Just (opA_power (lemma_A_l2 l2z l2zeta)) + + -- Case 2.2: b ≤ 0 and z ≥ -0.2 and zeta ≥ -0.2. + -- Region B^m. + | otherwise + = Just (opB_power (lemma_B_l2 l2z l2zeta)) + + where + (b, l2z) = operator_to_bl2z matA + (beta, l2zeta) = operator_to_bl2z matB + + logLambda a = logBase_double lambda a + l2z_minus_zeta = l2z / l2zeta -- λ[sup 2(/z/-ζ)] + + wlog_using op = + let (matA',matB') = action (matA, matB) op + maybe_op2 = step_lemma (matA',matB') + in + case maybe_op2 of + Nothing -> Just op + Just op2 -> Just (op * op2) + + with_shift k = + let (matA', matB') = shift_state k (matA, matB) + maybe_op2 = step_lemma (matA', matB') + in + case maybe_op2 of + Nothing -> Nothing + Just op2 -> Just (shift_sigma k op2) + +-- | Repeatedly apply the Step Lemma to the given state, until the +-- skew is 15 or less. +reduction :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r) => OperatorPair r -> Operator DRootTwo +reduction st = + case step_lemma st of + Nothing -> 1 + Just opG -> opG * opG' + where + opG' = reduction (action st opG) + +-- | Given a pair of ellipses, return a grid operator /G/ such that +-- the uprightness of each ellipse is greater than 1\/6. This is +-- essentially the same as 'reduction', except we do not assume that +-- the input operators have determinant 1. +to_upright :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r) => OperatorPair r -> Operator DRootTwo +to_upright (a,b) = opG + where + a' = a `scalardiv` (sqrt (det a)) + b' = b `scalardiv` (sqrt (det b)) + opG = reduction (a',b') + +-- ---------------------------------------------------------------------- +-- ** Action of special grid operators on convex sets + +-- | Apply a linear transformation /G/ to a point /p/. +point_transform :: (Ring r) => Operator r -> Point r -> Point r +point_transform opG (x,y) = (x',y') where + ((a,b), (c,d)) = fromOperator opG + x' = a * x + b * y + y' = c * x + d * y + +-- | Apply a special linear transformation /G/ to an ellipse /A/. This +-- results in the new ellipse /G/(/A/) = { /G/(/z/) | /z/ ∈ /A/ }. +ellipse_transform :: (Ring r, Adjoint r) => Operator r -> Ellipse r -> Ellipse r +ellipse_transform opG (Ellipse matA ctrA) = (Ellipse matA' ctrA') where + matA' = adj opG_inv * matA * opG_inv + ctrA' = point_transform opG ctrA + opG_inv = special_inverse opG + +-- | Apply a special grid operator /G/ to a characteristic function. +charfun_transform :: Operator DRootTwo -> CharFun -> CharFun +charfun_transform opG f = f' where + f' p = f (point_transform opG_inv p) + opG_inv = special_inverse opG + +-- | Apply a special linear transformation /G/ to a line +-- intersector. If the input line intersector was for a convex set +-- /A/, then the output line intersector is for the set /G/(/A/) +-- = { /G/(/z/) | /z/ ∈ /A/ }. +lineintersector_transform :: (Ring r) => Operator DRootTwo -> LineIntersector r -> LineIntersector r +lineintersector_transform opG intA = intA' + where + opG_inv = special_inverse opG + intA' v' w' = intA v w + where + v = point_transform opG_inv v' + w = point_transform opG_inv w' + +-- | Apply a special linear transformation /G/ to a convex set +-- /A/. This results in the new convex set /G/(/A/) = { /G/(/z/) | /z/ +-- ∈ /A/ }. +convex_transform :: (Ring r, Adjoint r, RootHalfRing r) => Operator DRootTwo -> ConvexSet r -> ConvexSet r +convex_transform opG (ConvexSet ellA tstA intA) = (ConvexSet ellA' tstA' intA') where + ellA' = ellipse_transform (op_fromDRootTwo opG) ellA + intA' = lineintersector_transform (op_fromDRootTwo opG) intA + tstA' = charfun_transform opG tstA + +-- ---------------------------------------------------------------------- +-- ** Bounding boxes + +-- | Calculate the bounding box for an ellipse. +boundingbox_ellipse :: (Floating r) => Ellipse r -> ((r, r), (r, r)) +boundingbox_ellipse (Ellipse matA ctrA) = ((x-w, x+w), (y-h, y+h)) where + (x,y) = ctrA + ((a, b), (c, d)) = fromOperator matA + w = sqrt d / sqrt_det + h = sqrt a / sqrt_det + sqrt_det = sqrt (det matA) + +-- | Calculate a bounding box for a convex set. Returns ((/x/₀, /x/₁), +-- (/y/₀, /y/₁)). +boundingbox :: (Floating r) => ConvexSet r -> ((r, r), (r, r)) +boundingbox (ConvexSet ell tst int) = boundingbox_ellipse ell + +-- ---------------------------------------------------------------------- +-- * Auxiliary functions + +-- | We write /x/ \`@within@\` (/a/,/b/) for /a/ ≤ /x/ ≤ /b/, or +-- equivalently, /x/ ∈ [/a/, /b/]. +within :: (Ord a) => a -> (a,a) -> Bool +within x (a,b) = a <= x && x <= b + +-- | Given an interval, return a slightly bigger one. +fatten_interval :: (Fractional a) => (a,a) -> (a,a) +fatten_interval (x,y) = (x - epsilon, y + epsilon) where + epsilon = 0.0001 * (y-x) + +-- | The constant λ = 1 + √2. +lambda :: (RootTwoRing r) => r +lambda = 1 + roottwo + +-- | The constant λ⁻¹ = √2 - 1. +lambda_inv :: (RootTwoRing r) => r +lambda_inv = roottwo - 1 + +-- | Return λ[sup /k/], where /k/ ∈ ℤ. This works in any 'RootTwoRing'. +-- +-- Note that we can't use '^', because it requires /k/ ≥ 0, nor '**', +-- because it requires the 'Floating' class. +lambdapower :: (RootTwoRing r) => Integer -> r +lambdapower k + | k >= 0 = lambda^k + | otherwise = lambda_inv^(-k) + +-- | Return (-1)[sup /k/], where /k/ ∈ ℤ. +signpower :: (Num r) => Integer -> r +signpower k + | even k = 1 + | otherwise = -1 + +-- | Given positive numbers /b/ and /x/, return (/n/, /r/) such that +-- +-- * /x/ = /r/ /b/[sup /n/] and +-- +-- * 1 ≤ /r/ < /b/. +-- +-- In other words, let /n/ = ⌊log[sub /b/] /x/⌋ and +-- /r/ = /x/ /b/[sup −/n/]. This can be more efficient than 'floor' +-- ('logBase' /b/ /x/) depending on the type; moreover, it also works +-- for exact types such as 'Rational' and 'QRootTwo'. +floorlog :: (Fractional b, Ord b) => b -> b -> (Integer, b) +floorlog b x + | x <= 0 = error "floorlog: argument not positive" + | 1 <= x && x < b = (0, x) + | 1 <= x*b && x < 1 = (-1, b*x) + | r < b = (2*n, r) + | otherwise = (2*n+1, r/b) + where + (n, r) = floorlog (b^2) x + +-- | A version of the natural logarithm that returns a 'Double'. The +-- logarithm of just about any value can fit into a 'Double'; so if +-- not a lot of precision is required in the mantissa, this function +-- is often faster than 'log'. +logBase_double :: (Fractional a, Real a) => a -> a -> Double +logBase_double b x + | b > 1 = y + | b <= 0 = 0/0 -- NaN + | b == 1 = 1/0 -- Infinity + | otherwise = - logBase_double (1/b) x + where + (n, r) = floorlog b x + y = fromInteger n + logBase (to_double b) (to_double r) + to_double = fromRational . toRational + +-- | The inner product of two points. +iprod :: (Num r) => Point r -> Point r -> r +iprod (x,y) (a,b) = x*a + y*b + +-- | Subtract two points. +point_sub :: (Num r) => Point r -> Point r -> Point r +point_sub (x,y) (a,b) = (x-a, y-b) + +-- | Calculute the inverse of an operator of determinant 1. Note: this +-- does not work correctly for operators whose determinant is not 1. +special_inverse :: (Ring r) => Operator r -> Operator r +special_inverse opG = opG_inv where + ((a,b), (c,d)) = fromOperator opG + opG_inv = det opG `scalarmult` toOperator ((d,-b), (-c,a)) + diff --git a/Quantum/Synthesis/GridSynth.hs b/Quantum/Synthesis/GridSynth.hs new file mode 100644 --- /dev/null +++ b/Quantum/Synthesis/GridSynth.hs @@ -0,0 +1,185 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +-- | This module implements the approximate single-qubit synthesis +-- algorithm of +-- +-- * N. J. Ross and P. Selinger, \"Optimal ancilla-free Clifford+/T/ +-- approximation of /z/-rotations\". . +-- +-- The algorithm is near-optimal in the following sense: it produces +-- an operator whose expected /T/-count exceeds the /T/-count of the +-- second-to-optimal solution to the approximate synthesis problem by +-- at most /O/(log(log(1/ε))). + +module Quantum.Synthesis.GridSynth where + +import Quantum.Synthesis.Ring +import Quantum.Synthesis.Ring.FixedPrec +import Quantum.Synthesis.Matrix +import Quantum.Synthesis.CliffordT +import Quantum.Synthesis.SymReal +import Quantum.Synthesis.GridProblems +import Quantum.Synthesis.Diophantine +import Quantum.Synthesis.StepComp +import Quantum.Synthesis.QuadraticEquation + +import System.Random +import Data.Number.FixedPrec + +-- ---------------------------------------------------------------------- +-- * Approximate synthesis + +-- ---------------------------------------------------------------------- +-- ** User-friendly functions + +-- | Output a unitary operator in the Clifford+/T/ group that +-- approximates /R/[sub /z/](θ) = [exp −/i/θ/Z/\/2] to within ε in the +-- operator norm. This operator can then be converted to a list of +-- gates with 'to_gates'. +-- +-- The parameters are: +-- +-- * a source of randomness /g/; +-- +-- * the angle θ; +-- +-- * the precision /b/ ≥ 0 in bits, such that ε = 2[sup -/b/]; +-- +-- * an integer that determines the amount of \"effort\" to put into +-- factoring. A larger number means more time spent on factoring. +-- A good default for this is 25. +-- +-- Note: the argument /theta/ is given as a symbolic real number. It +-- will automatically be expanded to as many digits as are necessary +-- for the internal calculation. In this way, the caller can specify, +-- e.g., an angle of pi\/128 @::@ 'SymReal', without having to worry +-- about how many digits of π to specify. +gridsynth :: (RandomGen g) => g -> Double -> SymReal -> Int -> U2 DOmega +gridsynth g prec theta effort = m where + (m, _, _) = gridsynth_stats g prec theta effort + +-- | A version of 'gridsynth' that returns a list of gates instead of a +-- matrix. +-- +-- Note: the list of gates will be returned in right-to-left order, +-- i.e., as in the mathematical notation for matrix multiplication. +-- This is the opposite of the quantum circuit notation. +gridsynth_gates :: (RandomGen g) => g -> Double -> SymReal -> Int -> [Gate] +gridsynth_gates g prec theta effort = synthesis_u2 (gridsynth g prec theta effort) + +-- | A version of 'gridsynth' that also returns some statistics: +-- log[sub 0.5] of the actual approximation error (or 'Nothing' if the +-- error is 0), and a data structure with information on the +-- candidates tried. +gridsynth_stats :: (RandomGen g) => g -> Double -> SymReal -> Int -> (U2 DOmega, Maybe Double, [(DOmega, DStatus)]) +gridsynth_stats g prec theta effort = dynamic_fixedprec2 digits f prec theta where + digits = ceiling (15 + 2 * prec * logBase 10 2) + f prec theta = gridsynth_internal g prec theta effort + +-- | Information about the status of an attempt to solve a Diophantine +-- equation. 'Success' means the Diophantine equation was solved; +-- 'Fail' means that it was proved that there was no solution; +-- 'Timeout' means that the question was not decided within the +-- allotted time. +data DStatus = Success | Fail | Timeout + deriving (Eq, Show) + +-- ---------------------------------------------------------------------- +-- * Implementation details + +-- ---------------------------------------------------------------------- +-- ** The ε-region + +-- | The ε-/region/ for given ε and θ is a convex subset of the closed +-- unit disk, given by [nobr /u/ ⋅ /z/ ≥ 1 - ε²\/2], where [nobr /z/ = +-- [exp −/i/θ\/2]], and “⋅” denotes the dot product of ℝ² (identified +-- with ℂ). +-- +-- \[center [image Re.png]] + +epsilon_region :: (Floating r, Ord r, RootHalfRing r, Quadratic r) => r -> r -> ConvexSet r +epsilon_region epsilon theta = ConvexSet ell tst int where + + -- A bounding ellipse for the ε-region. + ell = Ellipse mat ctr + ctr = (d*zx, d*zy) + mat = bmat * mmat * special_inverse bmat + mmat = toOperator ((ev1, 0), (0, ev2)) + bmat = toOperator ((zx, -zy), (zy, zx)) + ev1 = 4 * (1 / epsilon)^4 + ev2 = (1 / epsilon)^2 + + -- A line intersector for the ε-region. + int p v + | q == Nothing = (1, 0) + | vz == 0 && rhs <= 0 = (t0, t1) + | vz == 0 && otherwise = (1, 0) + | vz > 0 = (max t0 t2, t1) + | otherwise = (t0, min t1 t2) + where + a = iprod v v + b = 2 * iprod v p + c = iprod p p - 1 + q = quadratic (fromDRootTwo a) (fromDRootTwo b) (fromDRootTwo c) + Just (t0, t1) = q + + -- solve (p + tv) * z >= d + -- equivalently, t * vz >= d - pz + vz = iprod (point_fromDRootTwo v) z + rhs = d - iprod (point_fromDRootTwo p) z + t2 = rhs / vz + + -- The characteristic function of the ε-region. + tst (x,y) = x^2 + y^2 <= 1 && zx * fromDRootTwo x + zy * fromDRootTwo y >= d + + zx = cos (-theta/2) + zy = sin (-theta/2) + d = 1 - epsilon^2/2 + z = (zx, zy) + +-- ---------------------------------------------------------------------- +-- ** Main algorithm implementation + +-- | The internal implementation of the ellipse-based approximate +-- synthesis algorithm. The parameters are a source of randomness /g/, +-- the angle θ, the precision /b/ ≥ 0 in bits, and an amount of +-- \"effort\" to put into factoring. +-- +-- The outputs are a unitary operator in the Clifford+/T/ group that +-- approximates /R/[sub /z/](θ) to within ε in the operator norm; +-- log[sub 0.5] of the actual error, or 'Nothing' if the error is 0; +-- and the number of candidates tried. +-- +-- Note: the parameter θ must be of a real number type that has enough +-- precision to perform intermediate calculations; this typically +-- requires precision O(ε[sup 2]). A more user-friendly function that +-- selects the required precision automatically is 'gridsynth'. +gridsynth_internal :: forall r g.(RootHalfRing r, Ord r, Floating r, Adjoint r, Floor r, RealFrac r, Quadratic r, RandomGen g) => g -> r -> r -> Int -> (U2 DOmega, Maybe Double, [(DOmega, DStatus)]) +gridsynth_internal g prec theta effort = (uU, log_err, candidate_info) where + epsilon = 2 ** (-prec) + region = epsilon_region epsilon theta + candidates = gridpoints2_increasing region unitdisk + (uU, log_err, candidate_info) = first_solvable [] g candidates + + first_solvable candidate_info g [] = error "gridsynth: internal error: finite list of candidates?" + first_solvable candidate_info g (u : us) = case answer_t of + Just (Just t) -> let (uU, log_err) = with_successful_candidate u t in (uU, log_err, ((u, Success) : candidate_info)) + Just Nothing -> first_solvable ((u, Fail) : candidate_info) g2 us + Nothing -> first_solvable ((u, Timeout) : candidate_info) g2 us + where + (g1, g2) = split g + xi = real (1 - adj u * u) + answer_t = run_bounded effort $ diophantine_dyadic g1 xi + + with_successful_candidate u t = (uU, log_err) where + uU | denomexp (u + t) < denomexp (u + omega * t) + = matrix2x2 (u, -(adj t)) (t, adj u) + | otherwise + = matrix2x2 (u, -(adj (omega*t))) (omega*t, adj u) + log_err + | err <= 0 = Nothing + | otherwise = Just (logBase_double 0.5 err) + err = sqrt (real (hs_sqnorm (uU_fixed - zrot_fixed)) / 2) + uU_fixed = matrix_map fromDOmega uU + zrot_fixed = zrot (theta :: r) + diff --git a/Quantum/Synthesis/LaTeX.hs b/Quantum/Synthesis/LaTeX.hs --- a/Quantum/Synthesis/LaTeX.hs +++ b/Quantum/Synthesis/LaTeX.hs @@ -53,7 +53,7 @@ showlatex = show instance ShowLaTeX ZOmega where - showlatex (Omega a b c d) = format_signed_list list2 where + showlatex_p prec (Omega a b c d) = showParen (prec > 6) $ showString $ format_signed_list list2 where list = map signedunit [(a,"\\omega^3"),(b,"\\omega^2"),(c,"\\omega"),(d,"")] list2 = filter (\(s,a) -> s /= 0) list signedunit (a, u) @@ -75,7 +75,7 @@ cont ((_,a):t) = "-" ++ a ++ cont t instance (ShowLaTeX a, Nat n) => ShowLaTeX (Matrix n m a) where - showlatex (Matrix a) = "\\zmatrix{" ++ replicate m 'c' ++ "}{" ++ entries ++ "}" where + showlatex (Matrix a) = "\\begin{pmatrix}" ++ entries ++ "\\end{pmatrix}" where m = length (list_of_vector a) entries = concat $ list_of_vector $ vector_map showcolumn (vector_transpose a) showcolumn :: ShowLaTeX a => Vector m a -> String @@ -121,21 +121,24 @@ instance ShowLaTeX Double where showlatex x = printf "%0.10f" x +instance ShowLaTeX DOmega where + showlatex_p = showlatex_denomexp_p + -- This is an overlapping instance instance Nat n => ShowLaTeX (Matrix n m DOmega) where - showlatex = showlatex_denomexp + showlatex_p = showlatex_denomexp_p -- This is an overlapping instance instance Nat n => ShowLaTeX (Matrix n m DRComplex) where - showlatex = showlatex_denomexp + showlatex_p = showlatex_denomexp_p -- | Generic showlatex-like method that factors out a common -- denominator exponent. -showlatex_denomexp :: (WholePart a b, ShowLaTeX b, DenomExp a) => a -> String -showlatex_denomexp a - | k == 0 = showlatex b - | k == 1 = "\\frac{1}{\\sqrt{2}}" ++ showlatex b - | otherwise = "\\frac{1}{\\sqrt{2}^{" ++ show k ++ "}}" ++ showlatex b +showlatex_denomexp_p :: (WholePart a b, ShowLaTeX b, DenomExp a) => Int -> a -> ShowS +showlatex_denomexp_p d a + | k == 0 = showlatex_p d b + | k == 1 = showParen (d > 7) $ showString "\\frac{1}{\\sqrt{2}}" . showlatex_p 7 b + | otherwise = showParen (d > 7) $ showString ("\\frac{1}{\\sqrt{2}^{" ++ show k ++ "}}") . showlatex_p 7 b where (b, k) = denomexp_decompose a instance ShowLaTeX [Gate] where diff --git a/Quantum/Synthesis/Matrix.hs b/Quantum/Synthesis/Matrix.hs --- a/Quantum/Synthesis/Matrix.hs +++ b/Quantum/Synthesis/Matrix.hs @@ -353,6 +353,12 @@ infixl 7 `scalarmult` +-- | Division of an /m/×/n/-matrix by a scalar. +scalardiv :: (Fractional a) => Matrix m n a -> a -> Matrix m n a +scalardiv m x = matrix_map (/ x) m + +infixl 7 `scalardiv` + -- | Multiplication of /m/×/n/-matrices. We use a special symbol -- because /m/×/n/-matrices do not form a ring; only /n/×/n/-matrices -- form a ring (in which case the normal symbol \"'*'\" also works). diff --git a/Quantum/Synthesis/Newsynth.hs b/Quantum/Synthesis/Newsynth.hs --- a/Quantum/Synthesis/Newsynth.hs +++ b/Quantum/Synthesis/Newsynth.hs @@ -1,420 +1,28 @@ -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE ScopedTypeVariables #-} - --- | This module implements an efficient single-qubit Clifford+/T/ --- approximation algorithm. The algorithm is described here: +-- | This module provides backward compatibility with older versions +-- of the newsynth package. Formerly, it contained an implementation of +-- the single-qubit Clifford+/T/ approximation algorithm of -- -- * Peter Selinger. Efficient Clifford+/T/ approximation of -- single-qubit operators. . +-- +-- Since the new algorithm in "Quantum.Synthesis.GridSynth" is better +-- in all cases, we now simply provide a compatible interface to that +-- algorithm. +-- +-- New software should not use this module, and it may eventually be +-- removed. module Quantum.Synthesis.Newsynth where import Quantum.Synthesis.Ring -import Quantum.Synthesis.Ring.FixedPrec +import Quantum.Synthesis.SymReal import Quantum.Synthesis.Matrix import Quantum.Synthesis.CliffordT -import Quantum.Synthesis.EuclideanDomain -import Quantum.Synthesis.SymReal +import Quantum.Synthesis.GridSynth import System.Random -import Data.Number.FixedPrec --- ---------------------------------------------------------------------- --- * Miscellaneous functions - --- | A useful operation for the 'Maybe' monad, used to ensure that --- some condition holds (i.e., return 'Nothing' if the condition is --- false). To be used like this: --- --- > do --- > x <- something --- > y <- something_else --- > ensure (x > y) --- > ... -ensure :: Bool -> Maybe () -ensure True = Just () -ensure False = Nothing - --- | Return the head of a list, if non-empty, or else 'Nothing'. -maybe_head :: [a] -> Maybe a -maybe_head [] = Nothing -maybe_head (h:t) = Just h - --- | Exponentiation via repeated squaring, parameterized by a --- multiplication function and a unit. Given an associative --- multiplication function @*@ with unit @e@, the function 'power' --- @(*)@ /e/ /a/ /n/ efficiently computes /a/[sup /n/] = /a/ @*@ (/a/ --- @*@ (… @*@ (/a/ @*@ /e/)…)). -power :: (a -> a -> a) -> a -> a -> Integer -> a -power mul unit = aux where - aux x n - | n <= 0 = unit - | n == 1 = x - | odd n = x `mul` (x `aux` (n-1)) - | otherwise = y `mul` y where y = x `aux` (n `div` 2) - --- | Given positive numbers /b/ and /x/, return (/n/, /r/) such that --- --- * /x/ = /r/ /b/[sup /n/] and --- --- * 1 ≤ /r/ < /b/. --- --- In other words, let /n/ = ⌊log[sub /b/] /x/⌋ and --- /r/ = /x/ /b/[sup −/n/]. This can be more efficient than 'floor' --- ('logBase' /b/ /x/) depending on the type; moreover, it also works --- for exact types such as 'Rational' and 'QRootTwo'. -floorlog :: (Fractional b, Ord b) => b -> b -> (Integer, b) -floorlog b x - | x <= 0 = error "floorlog: argument not positive" - | 1 <= x && x < b = (0, x) - | 1 <= x*b && x < 1 = (-1, b*x) - | r < b = (2*n, r) - | otherwise = (2*n+1, r/b) - where - (n, r) = floorlog (b^2) x - --- ---------------------------------------------------------------------- --- * Randomized algorithms - --- | A combinator for turning a probabilistic function that succeeds --- with some small probability into a probabilistic function that --- always succeeds, by trying again and again. -keeptrying :: (RandomGen g) => (g -> Maybe a) -> (g -> a) -keeptrying f g = case f g1 of - Just res -> res - Nothing -> keeptrying f g2 - where - (g1, g2) = split g - --- | Like 'keeptrying', but also returns a count of the number of attempts. -keeptrying_count :: (RandomGen g) => (g -> Maybe a) -> (g -> (a, Integer)) -keeptrying_count f g = aux g 1 where - aux g n = case f g1 of - Just res -> (res, n) - Nothing -> aux g2 n1 - where - (g1, g2) = split g - !n1 = n + 1 - --- | A combinator for turning a probabilistic function that succeeds --- with some small probability into a probabilistic function that --- succeeds with a higher probability, by repeating it /n/ times. -try_for :: (RandomGen g) => Integer -> (g -> Maybe a) -> (g -> Maybe a) -try_for n f g - | n <= 0 = Nothing - | otherwise = case f g1 of - Just res -> Just res - Nothing -> try_for (n-1) f g2 - where - (g1, g2) = split g - --- ---------------------------------------------------------------------- --- * Square roots in ℤ[√2] - --- | Return a square root of an element of ℤ[√2], if such a square --- root exists, or else 'Nothing'. -zroottwo_root :: ZRootTwo -> Maybe ZRootTwo -zroottwo_root z@(RootTwo a b) = res where - d = a^2 - 2*b^2 - r = intsqrt d - x1 = intsqrt ((a + r) `div` 2) - x2 = intsqrt ((a - r) `div` 2) - y1 = intsqrt ((a - r) `div` 4) - y2 = intsqrt ((a + r) `div` 4) - w1 = RootTwo x1 y1 - w2 = RootTwo x2 y2 - w3 = RootTwo x1 (-y1) - w4 = RootTwo x2 (-y2) - res - | w1*w1 == z = Just w1 - | w2*w2 == z = Just w2 - | w3*w3 == z = Just w3 - | w4*w4 == z = Just w4 - | otherwise = Nothing - --- ---------------------------------------------------------------------- --- * Roots of −1 in ℤ[sub /p/] - --- | Input an integer /p/, and maybe output a root of −1 modulo /p/. --- This succeeds with probability at least 1\/2 if /p/ is a positive --- prime ≡ 1 (mod 4); otherwise, the success probability is --- unspecified (and may be 0). -root_minus_one_step :: (RandomGen g) => Integer -> g -> Maybe Integer -root_minus_one_step p g = do - let (b, _) = randomR (1, p-1) g - let h = power mul 1 b ((p-1) `div` 4) - ensure $ h `mul` h == p-1 -- succeeds with probability 1/2 - return h - where - mul :: Integer -> Integer -> Integer - mul a b = (a*b) `mod` p - --- | Input a positive prime /p/ ≡ 1 (mod 4), and output a root of −1. -root_minus_one :: (RandomGen g) => Integer -> g -> Integer -root_minus_one p = keeptrying (root_minus_one_step p) - --- ---------------------------------------------------------------------- --- * Solving a Diophantine equation - --- | Input ξ ∈ ℤ[√2], and maybe output some /t/ ∈ ℤ[ω] such that --- /t/[sup †]/t/ = ξ. If ξ ≥ 0, ξ[sup •] ≥ 0 and /p/ = ξ[sup •]ξ is a --- prime ≡ 1 (mod 4) in ℤ, then this succeeds with probability at least --- 1\/2. Otherwise, the success probability is unspecified and may be --- 0. -dioph_step :: (RandomGen g) => ZRootTwo -> g -> Maybe ZOmega -dioph_step xi g = do - h <- root_minus_one_step (norm xi) g - let s = euclid_gcd (fromInteger h+i) (fromZRootTwo xi) :: ZOmega - ss = zroottwo_of_zomega (adj s * s) - u = euclid_div xi ss - v <- zroottwo_root u - let t = fromZRootTwo v * s - ensure $ adj t * t == fromZRootTwo xi -- check the answer, just in case - return t - --- | Input ξ ∈ ℤ[√2] such that ξ ≥ 0, ξ[sup •] ≥ 0, and /p/ = --- ξ[sup •]ξ is a prime ≡ 1 (mod 4) in ℤ. Output /t/ ∈ ℤ[ω] such that --- /t/[sup †]/t/ = ξ. If the hypotheses are not satisfied, this will --- likely loop forever. -dioph :: (RandomGen g) => ZRootTwo -> g -> ZOmega -dioph xi = keeptrying (dioph_step xi) - --- ---------------------------------------------------------------------- --- * Approximations in ℤ[√2] - --- | Input two intervals [/x/₀, /x/₁] ⊆ ℝ and [/y/₀, /y/₁] ⊆ ℝ. Output --- a list of all points /z/ = /a/ + √2/b/ ∈ ℤ[√2] such that /z/ ∈ --- [/x/₀, /x/₁] and /z/[sup •] ∈ [/y/₀, /y/₁]. The list will be --- produced lazily, and will be sorted in order of increasing /z/. --- --- It is a theorem that there will be at least one solution if ΔxΔy ≥ (1 --- + √2)², and at most one solution if ΔxΔy < 1, where Δx = /x/₁ − /x/₀ ≥ 0 --- and Δy = /y/₁ − /y/₀ ≥ 0. Asymptotically, the expected number of --- solutions is ΔxΔy/\√8. --- --- This function is formulated so that the intervals can be specified --- exactly (using a type such as 'QRootTwo'), or approximately (using a --- type such as 'Double' or 'FixedPrec' /e/). -gridpoints :: (RootTwoRing r, Fractional r, Floor r, Ord r) => (r, r) -> (r, r) -> [ZRootTwo] -gridpoints (x0, x1) (y0, y1) - | dy <= 0 && dx > 0 = - map adj2 $ gridpoints (y0, y1) (x0, x1) - | dy >= lambda && even n = - map (lambdainv_n *) $ gridpoints (lambda_n*x0, lambda_n*x1) (lambda'_n*y0, lambda'_n*y1) - | dy >= lambda && odd n = - map (lambdainv_n *) $ gridpoints (lambda_n*x0, lambda_n*x1) (lambda'_n*y1, lambda'_n*y0) - | dy > 0 && dy < 1 && even n = - map (lambda_m *) $ gridpoints (lambdainv_m*x0, lambdainv_m*x1) (lambdainv'_m*y0, lambdainv'_m*y1) - | dy > 0 && dy < 1 && odd n = - map (lambda_m *) $ gridpoints (lambdainv_m*x0, lambdainv_m*x1) (lambdainv'_m*y1, lambdainv'_m*y0) - | otherwise = - [ RootTwo a b | a <- [amin..amax], b <- [bmin a..bmax a], test a b ] - where - dx = x1 - x0 - dy = y1 - y0 - (n, _) = floorlog lambda dy - m = -n - - lambda :: (RootTwoRing r) => r - lambda = 1 + roottwo - lambda' :: (RootTwoRing r) => r - lambda' = 1 - roottwo - lambdainv :: (RootTwoRing r) => r - lambdainv = -1 + roottwo - lambdainv' :: (RootTwoRing r) => r - lambdainv' = -1 - roottwo - lambda_m = lambda^m - lambda_n = lambda^n - lambda'_n = lambda'^n - lambdainv_m = lambdainv^m - lambdainv'_m = lambdainv'^m - lambdainv_n = lambdainv^n - - within x (x0, x1) = x0 <= x && x <= x1 - amin = ceiling_of ((x0 + y0) / 2) - amax = floor_of ((x1 + y1) / 2) - bmin a = ceiling_of ((fromInteger a - y1) / roottwo) - bmax a = floor_of ((fromInteger a - y0) / roottwo) - test a b = fromZRootTwo x `within` (x0, x1) && fromZRootTwo (adj2 x) `within` (y0, y1) - where x = RootTwo a b - --- | Input two intervals [/x/₀, /x/₁] ⊆ ℝ and [/y/₀, /y/₁] ⊆ ℝ and a --- source of randomness. Output a random element /z/ = /a/ + √2/b/ --- ∈ ℤ[√2] such that /z/ ∈ [/x/₀, /x/₁] and /z/[sup •] ∈ [/y/₀, --- /y/₁]. --- --- Note: the randomness will not be uniform. To ensure that the set of --- solutions is non-empty, we must have ΔxΔy ≥ (1 + √2)², where Δx = --- /x/₁ − /x/₀ ≥ 0 and Δy = /y/₁ − /y/₀ ≥ 0. If there are no solutions --- at all, the function will return 'Nothing'. --- --- This function is formulated so that the intervals can be specified --- exactly (using a type such as 'QRootTwo'), or approximately (using a --- type such as 'Double' or 'FixedPrec' /e/). -gridpoint_random :: (RootTwoRing r, Fractional r, Floor r, Ord r, RandomGen g) => (r, r) -> (r, r) -> g -> Maybe ZRootTwo -gridpoint_random (x0, x1) (y0, y1) g = z - where - dx = max 0 (x1 - x0) - dy = max 0 (y1 - y0) - area = dx * dy - n = floor_of (area + 1) - (i,_) = randomR (0, n-1) g - r = fromInteger i / fromInteger n - pts = gridpoints (x0 + r * dx, x1) (y0, y1) ++ gridpoints (x0, x1) (y0, y1) - z = maybe_head pts - --- | Input an integer /e/, two intervals [/x/₀, /x/₁] ⊆ ℝ and [/y/₀, --- /y/₁] ⊆ ℝ, and a source of randomness. Output random /z/ = /a/ + --- √2/b/ ∈ ℤ[√2] such that /a/ + √2/b/ ∈ [/x/₀, /x/₁], /a/ - √2/b/ ∈ --- [/y/₀, /y/₁], and /a/-/e/ is even. --- --- Note: the randomness will not be uniform. To ensure that the set of --- solutions is non-empty, we must have ΔxΔy ≥ 2(√2 + 1)², where Δx = --- /x/₁ − /x/₀ ≥ 0 and Δy = /y/₁ − /y/₀ ≥ 0. If there are no solutions --- at all, the function will return 'Nothing'. --- --- This function is formulated so that the intervals can be specified --- exactly (using a type such as 'QRootTwo'), or approximately (using a --- type such as 'Double' or 'FixedPrec' /e/). -gridpoint_random_parity :: (RootTwoRing r, Fractional r, Floor r, Ord r, RandomGen g) => Integer -> (r, r) -> (r, r) -> g -> Maybe ZRootTwo -gridpoint_random_parity e (x0,x1) (y0,y1) g = do - z' <- gridpoint_random (x0', x1') (-y1', -y0') g - return (roottwo * z' + fromInteger e2) - where - x0' = (x0 - e') / roottwo - x1' = (x1 - e') / roottwo - y0' = (y0 - e') / roottwo - y1' = (y1 - e') / roottwo - e' = fromInteger e2 - e2 = e `mod` 2 - --- ---------------------------------------------------------------------- --- * Approximate synthesis - --- ---------------------------------------------------------------------- --- ** The main algorithm - --- | The internal implementation of the approximate synthesis --- algorithm. The parameters are: --- --- * an angle θ, to implement a /R/[sub /z/](θ) = [exp −/i/θ/Z/\/2] --- gate; --- --- * a precision /p/ ≥ 0 in bits, such that ε = 2[sup -/p/]; --- --- * a source of randomness /g/. --- --- With some probability, output a unitary operator in the --- Clifford+/T/ group that approximates /R/[sub /z/](θ) to within ε in --- the operator norm. This operator can then be converted to a list of --- gates with 'to_gates'. Also output log[sub 0.1] of the actual --- error, or 'Nothing' if the error is 0. --- --- This implementation does not use seeding. --- --- As a special case, if the /R/[sub /z/](θ) is a Clifford operator --- (to within the given ε), always return this operator directly. --- --- Note: the parameter θ must be of a real number type that has enough --- precision to perform intermediate calculations; this typically --- requires precision O(ε[sup 2]). A more user-friendly function that --- does this automatically is 'newsynth'. -newsynth_step :: forall r g.(RealFrac r, Floating r, RootHalfRing r, Floor r, Adjoint r, RandomGen g) => r -> r -> g -> Maybe (U2 DOmega, Maybe Double) -newsynth_step prec theta = payload where - -- We are careful to do all computations that depend only on epsilon - -- and theta (but not g) outside of aux, to avoid re-computing them - -- with each attempt. - - -- Calculate ε. - epsilon = 2 ** (-prec) - - -- Convert prec to a Double - dprec = fromRational (toRational prec) - - -- Determine k. - const = 3 + 2 * logBase 2 (1 + sqrt 2) :: Double - k = ceiling (const + 2 * dprec) - scale = roottwo^k - - -- Normalize θ to be in [-π/4, π/4]. - n = round(theta / (pi/2)) - theta1 = theta - fromInteger n * pi/2 - - -- Describe the ε-region. - z @ (x,y) = (cos (theta1 / 2), -sin (theta1 / 2)) - e2 = 1 - epsilon^2/2 - e4 = 1 - epsilon^2/4 - z1 @ (x1,y1) = (e4 * x, e4 * y) - e' = epsilon / roottwo - f = e' * sqrt((1+e'/2)*(1-e'/2)) -- == sqrt(1-e4^2) - w @ (wx,wy) = (-f * y, f * x) - y_min = y1 - wy - y_max = y1 + wy - y'_min = y_min * scale - y'_max = y_max * scale - dx = (e4 - e2) * x - - find_uU_step = - -- As a special case, if (1,0) is in the ε-region, return the - -- identity operator. - if x >= e2 then \g -> Just 1 else aux - - -- The rest of the computation depends on the random seed g. - payload g = do - uU1 <- find_uU_step g - let uU = correct uU1 n - let err = calc_error uU theta - return (uU, err) - - aux g = do - -- Find a random grid point in the ε-region. - let (g0,g1) = split g - beta <- gridpoint_random (y'_min, y'_max) (-roothalf * scale, roothalf * scale) g0 - let - beta' = fromZRootTwo beta / scale - tmp = (beta' - e2 * y) / wy - x0 = e2 * x + tmp * wx - x1 = x0 + dx - x0' = x0 * scale - x1' = x1 * scale - (g2,g3) = split g1 - RootTwo c _ = beta - alpha <- gridpoint_random_parity (c+1) (x0', x1') (-roothalf * scale, roothalf * scale) g2 - - -- Calculate u, ξ, and solve Diophantine equation to calculate t. - let - u = (fromZRootTwo alpha) + i * (fromZRootTwo beta) :: ZOmega - xi = zroottwo_of_zomega (2^k - u * adj u) - t <- dioph_step xi g3 - - -- If Diophantine equation solved successfully, calculate matrix U. - let - u' = fromZOmega u * roothalf^k :: DOmega - t' = fromZOmega t * roothalf^k :: DOmega - uU1 = matrix2x2 (u', -(adj t')) - (t', (adj u')) - - return uU1 - - -- Correct for when θ wasn't in [-π/4, π/4]. - correct uU1 n = uU1 * rR^(n `mod` 8) where - rR = matrix2x2 (omega^7, 0) - (0, omega) - - -- Calculate the actual error. Since this is done lazily, this - -- incurs no overhead in case the error is not actually used. - calc_error uU theta = log_err where - uU_fixed = matrix_map fromDOmega uU :: U2 (Cplx r) - zrot_fixed = zrot theta :: U2 (Cplx r) - err = sqrt (real (hs_sqnorm (uU_fixed - zrot_fixed)) / 2) - log_err - | err <= 0 = Nothing - | otherwise = Just (log_double err / log 0.1) - --- ---------------------------------------------------------------------- --- ** User-friendly functions - --- | A user-friendly interface to the approximate synthesis +-- | Backward compatible interface to the approximate synthesis -- algorithm. The parameters are: -- -- * an angle θ, to implement a /R/[sub /z/](θ) = [exp −/i/θ/Z/\/2] @@ -429,26 +37,25 @@ -- operator can then be converted to a list of gates with -- 'to_gates'. -- --- This implementation does not use seeding. --- -- Note: the argument /theta/ is given as a symbolic real number. It -- will automatically be expanded to as many digits as are necessary -- for the internal calculation. In this way, the caller can specify, -- e.g., an angle of 'pi'\/128 @::@ 'SymReal', without having to worry -- about how many digits of π to specify. newsynth :: (RandomGen g) => Double -> SymReal -> g -> U2 DOmega -newsynth prec theta g = m where - (m, _, _) = newsynth_stats prec theta g +newsynth prec theta g = gridsynth g prec theta 25 -- | A version of 'newsynth' that also returns some statistics: -- log[sub 0.1] of the actual approximation error (or 'Nothing' if the -- error is 0), and the number of candidates tried. newsynth_stats :: (RandomGen g) => Double -> SymReal -> g -> (U2 DOmega, Maybe Double, Integer) -newsynth_stats prec theta g = dynamic_fixedprec2 digits f prec theta where - digits = ceiling (10 + 2 * prec * logBase 10 2) - f prec theta = (m, err, ct) where - ((m, err), ct) = keeptrying_count (newsynth_step prec theta) g - +newsynth_stats prec theta g = (op, err_d, n) where + (op, err_b, cinfo) = gridsynth_stats g prec theta 25 + err_d = case err_b of + Nothing -> Nothing + Just b -> Just (b * logBase 10 2) + n = fromIntegral (length cinfo) + -- | A version of 'newsynth' that returns a list of gates instead of a -- matrix. The inputs are the same as for 'newsynth'. -- @@ -456,4 +63,4 @@ -- i.e., as in the mathematical notation for matrix multiplication. -- This is the opposite of the quantum circuit notation. newsynth_gates :: (RandomGen g) => Double -> SymReal -> g -> [Gate] -newsynth_gates prec theta g = synthesis_u2 (newsynth prec theta g) +newsynth_gates prec theta g = gridsynth_gates g prec theta 25 diff --git a/Quantum/Synthesis/QuadraticEquation.hs b/Quantum/Synthesis/QuadraticEquation.hs new file mode 100644 --- /dev/null +++ b/Quantum/Synthesis/QuadraticEquation.hs @@ -0,0 +1,86 @@ +-- | This module provides a type class 'Quadratic', for solving +-- quadratic equations. + +module Quantum.Synthesis.QuadraticEquation ( + Quadratic (..) + ) where + +import Data.Number.FixedPrec +import Quantum.Synthesis.Ring + +-- | This type class provides a primitive method for solving quadratic +-- equations. For many floating-point or fixed-precision +-- representations of real numbers, using the usual \"quadratic +-- formula\" results in a significant loss of precision. Instances of +-- the 'Quadratic' class should provide an efficient high-precision +-- method when possible. +class Quadratic a where + -- | 'qroottwo_quadratic' /a/ /b/ /c/: solve the quadratic equation + -- /ax/² + /bx/ + /c/ = 0. Return the pair of solutions (/x/₁, /x/₂) + -- with /x/₁ ≤ /x/₂, or 'Nothing' if no solution exists. Note that + -- the coefficients /a/, /b/, and /c/ are taken to be of an exact + -- type; therefore instances have the opportunity to work with + -- infinite precision. + quadratic :: QRootTwo -> QRootTwo -> QRootTwo -> Maybe (a, a) + +-- ---------------------------------------------------------------------- +-- FixedPrec instance + +-- | Given /b/, /c/ ∈ ℚ[√2], consider the quadratic function /f/(/t/) +-- = /t/² + /b//t/ + /c/. +-- +-- * If /f/(/t/) = 0 has no real solutions, return 'Nothing'. +-- +-- * If /f/(/t/) = 0 has real solutions /t/₀ ≤ /t/₁, return /t/'₀, +-- /t/'₁ ∈ ℤ such that /t/'₀ ≤ /t/₀, /t/₁ ≤ /t/'₁, and |/t/'₀ - /t/₀|, +-- |/t/'₁ - /t/₁| ≤ 1. +int_quadratic :: QRootTwo -> QRootTwo -> Maybe (Integer, Integer) +int_quadratic b c + | radix < 0 = Nothing + | otherwise = Just (t0, t1) + where + radix = b^2/4 - c + tm = -b / 2 + rootradix' = intsqrt (floor_of radix) + t1' = floor_of tm + rootradix' + t1 + | is_solution1 (t1'+2) = t1'+2 + | is_solution1 (t1'+1) = t1'+1 + | otherwise = t1' + t0' = ceiling_of tm - rootradix' + t0 + | is_solution0 (t0'-2) = t0'-2 + | is_solution0 (t0'-1) = t0'-1 + | otherwise = t0' + is_solution1 x = f x' >= 0 && (f (x'-1) < 0 || x'-1 < tm) where + x' = fromInteger x + is_solution0 x = f x' >= 0 && (f (x'+1) < 0 || x'-1 > tm) where + x' = fromInteger x + f x = x^2 + b*x + c + +-- | Given /a/, /b/, /c/ ∈ ℚ[√2] with /a/ > 0, consider the quadratic +-- function /f/(/t/) = /a//t/² + /b//t/ + /c/. +-- +-- * If /f/(/t/) = 0 has no real solutions, return 'Nothing'. +-- +-- * If /f/(/t/) = 0 has real solutions /t/₀ ≤ /t/₁, return (/t/'₀, +-- /t/'₁) such that /t/'₀ ≤ /t/₀, /t/₁ ≤ /t/'₁, and |/t/'₀ - /t/₀|, +-- |/t/'₁ - /t/₁| ≤ 10[sup -/d/], where /d/ is the precision of the +-- fixed-point real number type. +qroottwo_quadratic_fixedprec :: (Precision e) => QRootTwo -> QRootTwo -> QRootTwo -> Maybe (FixedPrec e, FixedPrec e) +qroottwo_quadratic_fixedprec a b c + | False = Just (r, r) + | otherwise = do + (x0, x1) <- int_quadratic b' c' + return (fromInteger x0 / prec, fromInteger x1 / prec) + where + r = 0 + d = getprec r + prec = 10^d + prec' = 10^d + b' = prec' * b/a + c' = prec'^2 * c/a + q = int_quadratic b' c' + +instance (Precision e) => Quadratic (FixedPrec e) where + quadratic = qroottwo_quadratic_fixedprec diff --git a/Quantum/Synthesis/Ring.hs b/Quantum/Synthesis/Ring.hs --- a/Quantum/Synthesis/Ring.hs +++ b/Quantum/Synthesis/Ring.hs @@ -453,6 +453,27 @@ fromZRootTwo :: (RootTwoRing a) => ZRootTwo -> a fromZRootTwo (RootTwo x y) = fromInteger x + roottwo * fromInteger y +-- | Return a square root of an element of ℤ[√2], if such a square +-- root exists, or else 'Nothing'. +zroottwo_root :: ZRootTwo -> Maybe ZRootTwo +zroottwo_root z@(RootTwo a b) = res where + d = a^2 - 2*b^2 + r = intsqrt d + x1 = intsqrt ((a + r) `div` 2) + x2 = intsqrt ((a - r) `div` 2) + y1 = intsqrt ((a - r) `div` 4) + y2 = intsqrt ((a + r) `div` 4) + w1 = RootTwo x1 y1 + w2 = RootTwo x2 y2 + w3 = RootTwo x1 (-y1) + w4 = RootTwo x2 (-y2) + res + | w1*w1 == z = Just w1 + | w2*w2 == z = Just w2 + | w3*w3 == z = Just w3 + | w4*w4 == z = Just w4 + | otherwise = Nothing + -- ---------------------------------------------------------------------- -- ** The ring [bold D][√2] @@ -1042,16 +1063,35 @@ where k = lobit n +-- | Return 1 + the position of the leftmost \"1\" bit of a +-- non-negative 'Integer'. Do this in time O(/n/ log /n/), where /n/ +-- is the size of the integer (in digits). +hibit :: Integer -> Int +hibit 0 = 0 +hibit n = aux 1 where + aux k + | n >= 2^k = aux (2*k) + | otherwise = aux2 k (k `div` 2) -- 2^(k/2) <= n < 2^k + aux2 upper lower + | upper - lower < 2 = upper + | n >= 2^middle = aux2 upper middle + | otherwise = aux2 middle lower + where + middle = (upper + lower) `div` 2 + -- | For /n/ ≥ 0, return the floor of the square root of /n/. This is -- done using integer arithmetic, so there are no rounding errors. intsqrt :: (Integral n) => n -> n intsqrt n | n <= 0 = 0 - | otherwise = iterate 1 + | otherwise = iterate a where iterate m | m_sq <= n && m_sq + 2*m + 1 > n = m | otherwise = iterate ((m + n `div` m) `div` 2) where m_sq = m*m + a = 2^(b `div` 2) + b = hibit (fromIntegral n) + diff --git a/Quantum/Synthesis/StepComp.hs b/Quantum/Synthesis/StepComp.hs new file mode 100644 --- /dev/null +++ b/Quantum/Synthesis/StepComp.hs @@ -0,0 +1,157 @@ +-- | This module provides /step computations/. These are computations +-- that can be run, stopped, resumed, parallelized, and/or bounded in +-- runtime. + +module Quantum.Synthesis.StepComp where + +-- ---------------------------------------------------------------------- +-- * A monad for step computations + +-- | A step computation can be run for a specified number of steps, +-- stopped, continued, and interleaved. Such a computation produces +-- \"ticks\" at user-defined intervals, which must be consumed by the +-- environment for the computation to continue. +data StepComp a = + Done a -- ^ Terminate with a result. + | Tick (StepComp a) -- ^ Produce a \"tick\", then resume the + -- computation. + +instance Monad StepComp where + return a = Done a + Done a >>= g = g a + Tick f >>= g = Tick (f >>= g) + +instance Show a => Show (StepComp a) where + show (Done a) = "Done(" ++ show a ++ ")" + show (Tick c) = "Incomplete" + +-- ---------------------------------------------------------------------- +-- * Basic operations + +-- | Issue a single tick. +tick :: StepComp () +tick = Tick (Done ()) + +-- | Run the step computation for one step. +untick :: StepComp a -> StepComp a +untick (Done a) = (Done a) +untick (Tick c) = c + +-- | Fast-forward a computation by /n/ steps. This is essentially +-- equivalent to doing /n/ 'untick' operations. +forward :: Int -> StepComp a -> StepComp a +forward 0 c = c +forward n (Done a) = Done a +forward n c = forward (n-1) (untick c) + +-- | Check whether a step computation is completed. +is_done :: StepComp a -> Bool +is_done (Done a) = True +is_done (Tick c) = False + +-- | Retrieve the result of a completed step computation (or 'Nothing' +-- if it is incomplete). +get_result :: StepComp a -> Maybe a +get_result (Done a) = Just a +get_result (Tick c) = Nothing + +-- | Run a subsidiary computation for up to /n/ steps, translated into +-- an equal number of steps of the parent computation. +subtask :: Int -> StepComp a -> StepComp (StepComp a) +subtask n c | n <= 0 = Done c +subtask n (Done a) = Done (Done a) +subtask n (Tick c) = Tick (subtask (n-1) c) + +-- | Run a subtask, speeding it up by a factor of /n/ ≥ 1. Every 1 tick of +-- the calling task corresponds to up to /n/ ticks of the subtask. +speedup :: Int -> StepComp a -> StepComp a +speedup n (Done a) = Done a +speedup n (Tick c) = do + tick + speedup n (forward (n-1) c) + +-- | Run two step computations in parallel, until one branch +-- terminates. Tick allocation is associative: each tick of the +-- parent function translates into one tick for each subcomputation. +-- Therefore, when running, e.g., three subcomputations in parallel, +-- they will each receive an approximately equal number of ticks. +parallel :: StepComp a -> StepComp b -> StepComp (Either (a, StepComp b) (StepComp a, b)) +parallel (Done a) c = Done (Left (a, c)) +parallel c (Done b) = Done (Right (c, b)) +parallel (Tick c) (Tick c') = Tick (parallel c c') + +-- | Wrap a step computation to return the number of steps, in +-- addition to the result. +with_counter :: StepComp a -> StepComp (a, Int) +with_counter c = aux 0 c where + aux n (Done a) = return (a, n) + aux n (Tick c) = do + n `seq` tick + aux (n+1) c + +-- ---------------------------------------------------------------------- +-- ** Run functions + +-- | Run a step computation until it finishes. +run :: StepComp a -> a +run (Done a) = a +run (Tick c) = run c + +-- | Run a step computation until it finishes, and also return the +-- number of steps it took. +run_with_steps :: StepComp a -> (a, Int) +run_with_steps = run . with_counter + +-- | Run a step computation for at most /n/ steps. +run_bounded :: Int -> StepComp a -> Maybe a +run_bounded n = get_result . forward n + +-- ---------------------------------------------------------------------- +-- * Other operations + +-- | Do nothing, forever. +diverge :: StepComp a +diverge = tick >> diverge + +-- | Run two step computations in parallel. The first one to complete +-- becomes the result of the computation. +parallel_first :: StepComp a -> StepComp a -> StepComp a +parallel_first c1 c2 = do + r <- parallel c1 c2 + case r of + Left (a, _) -> return a + Right (_, a) -> return a + +-- | Run two step computations in parallel. If either computation +-- returns 'Nothing', return 'Nothing'. Otherwise, return the pair of +-- results. +parallel_maybe :: StepComp (Maybe a) -> StepComp (Maybe b) -> StepComp (Maybe (a,b)) +parallel_maybe c1 c2 = do + res <- parallel c1 c2 + case res of + Left (Nothing, c2) -> return Nothing + Right (c1, Nothing) -> return Nothing + Left (Just a, c2) -> do + b <- c2 + case b of + Nothing -> return Nothing + Just b -> return (Just (a,b)) + Right (c1, Just b) -> do + a <- c1 + case a of + Nothing -> return Nothing + Just a -> return (Just (a,b)) + +-- | Run a list of step computations in parallel. If any computation +-- returns 'Nothing', return 'Nothing'. Otherwise, return the list of +-- results. +parallel_list_maybe :: [StepComp (Maybe a)] -> StepComp (Maybe [a]) +parallel_list_maybe [] = return (Just []) +parallel_list_maybe (h:t) = do + res <- parallel_maybe h c2 + return $ do + (h',t') <- res + return (h':t') + where + c2 = parallel_list_maybe t + diff --git a/images/Re.png b/images/Re.png new file mode 100644 Binary files /dev/null and b/images/Re.png differ diff --git a/images/area.png b/images/area.png new file mode 100644 Binary files /dev/null and b/images/area.png differ diff --git a/images/bz.png b/images/bz.png new file mode 100644 Binary files /dev/null and b/images/bz.png differ diff --git a/images/ellipse-rectangle.png b/images/ellipse-rectangle.png new file mode 100644 Binary files /dev/null and b/images/ellipse-rectangle.png differ diff --git a/images/gridop-A.png b/images/gridop-A.png new file mode 100644 Binary files /dev/null and b/images/gridop-A.png differ diff --git a/images/gridop-Ai.png b/images/gridop-Ai.png new file mode 100644 Binary files /dev/null and b/images/gridop-Ai.png differ diff --git a/images/gridop-B.png b/images/gridop-B.png new file mode 100644 Binary files /dev/null and b/images/gridop-B.png differ diff --git a/images/gridop-Bi.png b/images/gridop-Bi.png new file mode 100644 Binary files /dev/null and b/images/gridop-Bi.png differ diff --git a/images/gridop-K.png b/images/gridop-K.png new file mode 100644 Binary files /dev/null and b/images/gridop-K.png differ diff --git a/images/gridop-R.png b/images/gridop-R.png new file mode 100644 Binary files /dev/null and b/images/gridop-R.png differ diff --git a/images/gridop-S.png b/images/gridop-S.png new file mode 100644 Binary files /dev/null and b/images/gridop-S.png differ diff --git a/images/gridop-Si.png b/images/gridop-Si.png new file mode 100644 Binary files /dev/null and b/images/gridop-Si.png differ diff --git a/images/gridop-V.png b/images/gridop-V.png new file mode 100644 Binary files /dev/null and b/images/gridop-V.png differ diff --git a/images/gridop-X.png b/images/gridop-X.png new file mode 100644 Binary files /dev/null and b/images/gridop-X.png differ diff --git a/images/gridop-Z.png b/images/gridop-Z.png new file mode 100644 Binary files /dev/null and b/images/gridop-Z.png differ diff --git a/images/gridop.png b/images/gridop.png new file mode 100644 Binary files /dev/null and b/images/gridop.png differ diff --git a/images/skew.png b/images/skew.png new file mode 100644 Binary files /dev/null and b/images/skew.png differ diff --git a/newsynth.cabal b/newsynth.cabal --- a/newsynth.cabal +++ b/newsynth.cabal @@ -7,7 +7,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 0.1.1.0 +version: 0.2 -- A short (one-line) description of the package. synopsis: Exact and approximate synthesis of quantum circuits @@ -19,10 +19,10 @@ quantum circuits over the Clifford+T gate set. This includes, among other things: . - * "Quantum.Synthesis.Newsynth": an efficient single-qubit - approximate synthesis algorithm. From P. Selinger, \"Efficient - Clifford+T approximation of single-qubit operators\", - . + * "Quantum.Synthesis.GridSynth": an efficient single-qubit + approximate synthesis algorithm. From N. J. Ross and P. Selinger, + \"Optimal ancilla-free Clifford+/T/ approximation of + /z/-rotations\", . . * "Quantum.Synthesis.MultiQubitSynthesis": multi-qubit exact synthesis algorithms. From B. Giles and P. Selinger, \"Exact @@ -53,14 +53,14 @@ license-file: LICENSE -- The package author(s). -author: Peter Selinger +author: Neil J. Ross, Peter Selinger -- An email address to which users can send suggestions, bug reports, and -- patches. maintainer: selinger@mathstat.dal.ca -- A copyright notice. -copyright: Copyright (c) 2012-2014 Peter Selinger +copyright: Copyright (c) 2012-2014 Neil J. Ross and Peter Selinger -- A classification category for future use by the package catalogue -- Hackage. These categories have not yet been specified, but the @@ -79,21 +79,21 @@ library -- Modules exported by the library. - exposed-modules: Quantum.Synthesis.Newsynth, Quantum.Synthesis.Matrix, Quantum.Synthesis.LaTeX, Quantum.Synthesis.RotationDecomposition, Quantum.Synthesis.ArcTan2, Quantum.Synthesis.EulerAngles, Quantum.Synthesis.EuclideanDomain, Quantum.Synthesis.SymReal, Quantum.Synthesis.Ring, Quantum.Synthesis.Clifford, Quantum.Synthesis.MultiQubitSynthesis, Quantum.Synthesis.CliffordT, Quantum.Synthesis.Ring.FixedPrec, Quantum.Synthesis.Ring.SymReal + exposed-modules: Quantum.Synthesis.Newsynth, Quantum.Synthesis.Matrix, Quantum.Synthesis.LaTeX, Quantum.Synthesis.RotationDecomposition, Quantum.Synthesis.ArcTan2, Quantum.Synthesis.EulerAngles, Quantum.Synthesis.EuclideanDomain, Quantum.Synthesis.SymReal, Quantum.Synthesis.Ring, Quantum.Synthesis.Clifford, Quantum.Synthesis.MultiQubitSynthesis, Quantum.Synthesis.CliffordT, Quantum.Synthesis.Ring.FixedPrec, Quantum.Synthesis.Ring.SymReal, Quantum.Synthesis.GridProblems Quantum.Synthesis.GridSynth Quantum.Synthesis.Diophantine Quantum.Synthesis.StepComp Quantum.Synthesis.QuadraticEquation -- Modules included in this library but not exported. -- other-modules: -- Other library packages from which modules are imported. - build-depends: base ==4.6.*, random ==1.0.*, fixedprec ==0.2.*, superdoc ==0.1.* + build-depends: base ==4.6.*, random ==1.0.*, fixedprec >= 0.2.2 && < 0.3, superdoc ==0.1.*, containers ==0.5.* -- Additional options for GHC when the package is built with -- profiling enabled. ghc-prof-options: -prof -auto-all -executable newsynth +executable gridsynth -- .hs or .lhs file containing the Main module. - main-is: newsynth.hs + main-is: gridsynth.hs -- Root directories for the module hierarchy. hs-source-dirs: programs diff --git a/programs/gridsynth.hs b/programs/gridsynth.hs new file mode 100644 --- /dev/null +++ b/programs/gridsynth.hs @@ -0,0 +1,418 @@ +{-# LANGUAGE BangPatterns #-} + +-- | This module provides a command line interface to the +-- single-qubit approximate synthesis algorithm. + +module Main where + +import Quantum.Synthesis.SymReal +import Quantum.Synthesis.CliffordT +import Quantum.Synthesis.Ring +import Quantum.Synthesis.Matrix +import Quantum.Synthesis.LaTeX +import Quantum.Synthesis.GridSynth +import Quantum.Synthesis.GridProblems + +import CommandLine + +-- import other stuff +import Control.Monad +import Data.Time +import System.Console.GetOpt +import System.Environment +import System.Exit +import System.IO +import System.Random +import Text.Printf + +-- ---------------------------------------------------------------------- +-- * Option processing + +-- | A data type to hold values set by command line options. +data Options = Options { + opt_digits :: Maybe Double, -- ^ Requested precision in decimal digits (default: 10). + opt_theta :: Maybe SymReal, -- ^ The angle θ to approximate. + opt_effort :: Int, -- ^ The amount of \"effort\" to spend on factoring. + opt_hex :: Bool, -- ^ Output operator in hex coding? (default: ASCII). + opt_stats :: Bool, -- ^ Output statistics? + opt_latex :: Bool, -- ^ Use LaTeX format? + opt_table :: Bool, -- ^ Generate the table of results for the paper? + opt_count :: Maybe Int, -- ^ Repeat count for --table mode (default: 50). + opt_rseed :: Maybe StdGen -- ^ An optional random seed. +} deriving Show + +-- | The initial default options. +defaultOptions :: Options +defaultOptions = Options + { opt_digits = Nothing, + opt_theta = Nothing, + opt_effort = 25, + opt_hex = False, + opt_stats = False, + opt_latex = False, + opt_table = False, + opt_count = Nothing, + opt_rseed = Nothing + } + +-- | The list of command line options, in the format required by 'getOpt'. +options :: [OptDescr (Options -> IO Options)] +options = + [ Option ['h'] ["help"] (NoArg help) "print usage info and exit", + Option ['d'] ["digits"] (ReqArg digits "") "set precision in decimal digits (default: 10)", + Option ['b'] ["bits"] (ReqArg bits "") "set precision in bits", + Option ['e'] ["epsilon"] (ReqArg epsilon "") "set precision as epsilon (default: 1e-10)", + Option ['f'] ["effort"] (ReqArg effort "\"\"") "how hard to try to factor (default: 25)", + Option ['x'] ["hex"] (NoArg hex) "output hexadecimal coding (default: ASCII)", + Option ['s'] ["stats"] (NoArg stats) "output statistics", + Option ['l'] ["latex"] (NoArg latex) "use LaTeX output format", + Option ['t'] ["table"] (NoArg table) "generate the table of results for the article", + Option ['c'] ["count"] (ReqArg count "") "repeat count for --table mode (default: 50)", + Option ['r'] ["rseed"] (ReqArg rseed "\"\"") "set optional random seed (default: random)" + ] + where + help :: Options -> IO Options + help o = do + usage + exitSuccess + + digits :: String -> Options -> IO Options + digits string o = + case parse_double string of + Just n | n >= 0 -> return o { opt_digits = Just n } + Just n -> optfail ("Number of digits must not be negative -- " ++ string ++ "\n") + _ -> optfail ("Invalid digits -- " ++ string ++ "\n") + + bits :: String -> Options -> IO Options + bits string o = + case parse_double string of + Just n | n >= 0 -> return o { opt_digits = Just (n * logBase 10 2) } + Just n -> optfail ("Number of bits must not be negative -- " ++ string ++ "\n") + _ -> optfail ("Invalid bits -- " ++ string ++ "\n") + + epsilon :: String -> Options -> IO Options + epsilon string o = + case parse_double string of + Just eps | eps < 1 && eps > 0 -> return o { opt_digits = Just (-logBase 10 eps) } + Just n -> optfail ("Epsilon must be between 0 and 1 -- " ++ string ++ "\n") + _ -> optfail ("Invalid epsilon -- " ++ string ++ "\n") + + effort :: String -> Options -> IO Options + effort string o = + case parse_int string of + Just e | e > 0 -> return o { opt_effort = e } + Just e -> optfail ("Effort must be positive -- " ++ string ++ "\n") + _ -> optfail ("Invalid effort -- " ++ string ++ "\n") + + hex :: Options -> IO Options + hex o = return o { opt_hex = True } + + stats :: Options -> IO Options + stats o = return o { opt_stats = True } + + latex :: Options -> IO Options + latex o = return o { opt_latex = True } + + table :: Options -> IO Options + table o = return o { opt_table = True } + + count :: String -> Options -> IO Options + count string o = + case parse_int string of + Just n | n >= 1 -> return o { opt_count = Just n } + Just n -> optfail ("Invalid count, must be positive -- " ++ string ++ "\n") + _ -> optfail ("Invalid count -- " ++ string ++ "\n") + + rseed :: String -> Options -> IO Options + rseed string o = + case reads string of + [(g, "")] -> return o { opt_rseed = Just g } + _ -> optfail ("Invalid random seed -- " ++ string ++ "\n") + +-- | Process /argv/-style command line options into an 'Options' structure. +dopts :: [String] -> IO Options +dopts argv = do + let (o, args, errs) = getOpt Permute options argv + opts <- foldM (flip id) defaultOptions o + when (errs /= []) $ do + optfail (concat errs) + case args of + [] -> return opts + [string] -> do + case parse_SymReal string of + Just theta -> return opts { opt_theta = Just theta } + _ -> optfail ("Invalid theta -- " ++ string ++ "\n") + h1:h2:[] -> optfail ("Too many non-option arguments -- " ++ h1 ++ ", " ++ h2 ++ "\n") + h1:h2:_ -> optfail ("Too many non-option arguments -- " ++ h1 ++ ", " ++ h2 ++ "...\n") + +-- | Print usage message to 'stdout'. +usage :: IO () +usage = do + putStr (usageInfo header options) + where header = + "Usage: gridsynth [OPTION...] \n" + ++ "Arguments:\n" + ++ " z-rotation angle. May be symbolic, e.g. pi/128\n" + ++ "Options:" + +-- ---------------------------------------------------------------------- +-- * The main function + +-- | Main function: read options, then execute the appropriate tasks. +main :: IO() +main = do + -- Read options. + argv <- getArgs + options <- dopts argv + case opt_table options of + False -> main_default options + True -> main_maketable options + +-- ---------------------------------------------------------------------- +-- ** Default main + +-- | The default task for the main function: synthesize one angle, for +-- one given precision, possibly with outputting some statistics. +main_default :: Options -> IO() +main_default options = do + let digits = case opt_digits options of + Nothing -> 10 + Just d -> d + let prec = digits * logBase 2 10 + theta <- case opt_theta options of + Nothing -> optfail "Missing argument: theta.\n" + Just t -> return t + case opt_count options of + Nothing -> return () + Just c -> optfail "Option -c is only supported with --table.\n" + let exponent = ceiling digits + let l = opt_latex options + let effort = opt_effort options + + -- Set random seed. + g <- case opt_rseed options of + Nothing -> newStdGen + Just g -> return g + + -- Payload. + t0 <- getCurrentTime + let (m,err,cinfo) = gridsynth_stats g prec theta effort + gates = to_gates m + if opt_hex options then + printf "%x\n" (convert gates :: Integer) + else if opt_latex options then + putStrLn (if gates == [] then "I" else showlatex gates) + else + putStrLn (if gates == [] then "I" else convert gates) + t1 <- getCurrentTime + + -- Print optional statistics + let ct = length cinfo + let tcount = length $ filter (==T) gates + let ulower = last [ u | (u, status) <- cinfo, status /= Fail ] + let klower = fromInteger (denomexp ulower) + let tlower = if klower == 0 then 0 else 2*klower - 2 + let secs = diffUTCTime t1 t0 + let err_d = case err of + Nothing -> Nothing + Just x -> Just (x * logBase 10 2) + + when (opt_stats options) $ do + putStrLn ("Random seed: " ++ show g) + putStrLn ("T-count: " ++ show tcount) + putStrLn ("Lower bound on T-count: " ++ show tcount) + putStrLn ("Theta: " ++ showf l theta) + putStrLn ("Epsilon: " ++ showf_exp l 10 exponent (Just digits)) + putStrLn ("Matrix: " ++ showf l m) + putStrLn ("Actual error: " ++ showf_exp l 10 exponent err_d) + putStrLn ("Runtime: " ++ show secs) + putStrLn ("Candidates tried: " ++ show ct ++ " (" + ++ show (length [u | (u, Fail) <- cinfo]) ++ " failed, " + ++ show (length [u | (u, Timeout) <- cinfo]) ++ " timed out, " + ++ show (length [u | (u, Success) <- cinfo]) ++ " succeeded)") + putStrLn ("Time/candidate: " ++ show (secs / fromIntegral ct)) + +-- ---------------------------------------------------------------------- +-- ** Generate output in LaTeX table format + +-- | Run one instance of the algorithm, using the given theta, and +-- measuring various things including the running time. Note: here, +-- the precision is expressed in /decimal/, not binary, digits. +-- +-- The inputs are, respectively: a source of randomness, the angle θ, +-- the precision in decimal digits, and an amount of effort to spend +-- on factoring. The outputs are, respectively: the approximating +-- operator /U/; the approximating circuit, log[sub 0.5] of the actual +-- approximation error (or 'Nothing' if the error is 0), the number of +-- candidates tried, the /T/-count of /U/, the computed lower bound +-- for the /T/-count, and the runtime in seconds. +one_run :: (RandomGen g, Show g) => g -> SymReal -> Double -> Int -> IO (U2 DOmega, [Gate], Maybe Double, Int, Int, Int, Double) +one_run g theta prec_d effort = do + let !prec = prec_d * logBase 2 10 + let !exponent = floor prec_d + putStrLn ("% Epsilon: " ++ show_exp 10 exponent (Just prec_d)) + putStrLn ("% Theta: " ++ show theta) + putStrLn ("% Random seed: " ++ show g) + t0 <- getCurrentTime + let (op, err, cinfo) = gridsynth_stats g prec theta effort + circ = synthesis_u2 op + tcount = length $ filter (==T) circ + putStrLn ("% T-count: " ++ show tcount) + t1 <- getCurrentTime + let secs = diffUTCTime t1 t0 + ct = length cinfo + -- find the first candidate that *might* have succeeded - this gives + -- a lower bound on the shorest possible T-count. + ulower = last [ u | (u, status) <- cinfo, status /= Fail ] + klower = fromInteger (denomexp ulower) + tlower = if klower == 0 then 0 else 2*klower - 2 + ((u, _), (t, _)) = fromOperator op + let err_d = case err of + Nothing -> Nothing + Just x -> Just (x * logBase 10 2) + putStrLn ("% Lower bound on T-count: " ++ show tlower) + putStrLn ("% Circuit: " ++ if circ == [] then "I" else convert circ) + putStrLn ("% u: " ++ showlatex u) + putStrLn ("% t: " ++ showlatex t) + putStrLn ("% Actual error: " ++ show_exp 10 exponent err_d) + putStrLn ("% Runtime: " ++ show secs) + putStrLn ("% Candidates tried: " ++ show ct ++ " (" + ++ show (length [u | (u, Fail) <- cinfo]) ++ " failed, " + ++ show (length [u | (u, Timeout) <- cinfo]) ++ " timed out, " + ++ show (length [u | (u, Success) <- cinfo]) ++ " succeeded)") + putStrLn ("% Time/candidate: " ++ show (secs / fromIntegral ct)) + putStrLn "" + hFlush stdout + return (op, circ, err, ct, tcount, tlower, fromRational (toRational secs)) + +-- | Repeat the algorithm /n/ times with the same parameters but +-- random angles, to average things like running time. The inputs are, +-- respectively: a source of randomness, a repeat count, the precision +-- in decimal digits,, and an amount of effort to spend on factoring. +many_runs :: (RandomGen g, Show g) => g -> Int -> Double -> Int -> IO () +many_runs g n prec_d effort = do + let gs = take n $ expand_seed g + results <- sequence $ do + g <- gs + return $ do + let (theta', g') = randomR (0, 2047) g + let theta = fromInteger theta' * pi / 2048 :: SymReal + one_run g' theta prec_d effort + -- Output the LaTeX of one row of the table + let (_,_,err,_,tcount,tlower,_) = head results + total_time = sum [ t | (_,_,_,_,_,_,t) <- results ] + total_candidates = sum [ ct | (_,_,_,ct,_,_,_) <- results ] + avg_time = total_time / fromIntegral n + avg_candidates = fromIntegral total_candidates / fromIntegral n :: Double + time_per_candidate = total_time / fromIntegral total_candidates + err_d = case err of + Nothing -> Nothing + Just x -> Just (x * logBase 10 2) + exponent = floor prec_d + putStrPad 30 (showlatex_exp 5 exponent (Just prec_d) ++ " &") + putStrLn ("% Epsilon") + putStrPad 30 (show tcount ++ " &") + putStrLn ("% T-count") + putStrPad 30 ("\\geq " ++ show tlower ++ " &") + putStrLn ("% Lower bound on T-count") + putStrPad 30 (showlatex_exp 5 exponent err_d ++ " &") + putStrLn ("% Actual error") + putStrPad 30 (printf "%0.4fs" avg_time ++ " &") + putStrLn ("% Runtime, averaged over " ++ show n ++ " runs") + putStrPad 30 (printf "%0.1f" avg_candidates ++ " &") + putStrLn ("% Candidates tried, averaged over " ++ show n ++ " runs") + putStrPad 30 (printf "%0.4fs" time_per_candidate ++ " \\\\") + putStrLn ("% Time per candidate, averaged over " ++ show n ++ " runs") + putStrLn "" + putStrLn "% ----------------------------------------------------------------------" + putStrLn "" + hFlush stdout + return () + +-- | Generate the table of \"Experimental Results\" used in the +-- article. +main_maketable :: Options -> IO () +main_maketable options = do + -- Read some parameters. + let theta = case opt_theta options of + Nothing -> pi/128 + Just t -> t + let count = case opt_count options of + Nothing -> 50 + Just c -> c + let precisions = case opt_digits options of + Nothing -> [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 500, 1000] + Just d -> [d] + let effort = opt_effort options + + -- Set random seed. + g <- case opt_rseed options of + Nothing -> newStdGen + Just g -> return g + putStrLn ("% Initial random seed: " ++ show g) + putStrLn "" + + -- Expand random seed. + let gs = expand_seed g + + -- Payload. + sequence_ $ do + (prec_d, g) <- zip precisions gs + return $ do + let (g1, g2) = split g + one_run g1 theta prec_d effort + many_runs g2 count prec_d effort + +-- ---------------------------------------------------------------------- +-- * Miscellaneous + +-- | Round a 'RealFrac' value to the given number of decimals. +round_to :: (RealFrac r) => Integer -> r -> r +round_to n x = fromInteger (round (x * 10^n)) / 10^n + +-- | Show the number 10[sup -/x/] in the format 10^(-n) or +-- 1.23*10^(-n), with precision /d/ and exponent -/n/. A value of +-- 'Nothing' is treated as 0. +-- +-- For example, display @0.316*10^(-13)@ instead of @10^(-13.5)@. +show_exp :: (Show r, RealFrac r, Floating r, PrintfArg r) => Integer -> Integer -> Maybe r -> String +show_exp d n x + | y == 1 = "10^(" ++ show (-n) ++ ")" + | otherwise = printf ("%." ++ show d ++ "f") y ++ "*10^(" ++ show (-n) ++ ")" + where + y = case x of + Nothing -> 0 + Just x -> round_to d (10 ** (fromInteger n - x)) + +-- | Show the number 10[sup -/x/] in the format @10^{-n}@ or +-- @1.23\\cdot 10^{-n}@, with precision /d/ and exponent -/n/. A value +-- of 'Nothing' is treated as 0. +showlatex_exp :: (Show r, RealFrac r, Floating r, PrintfArg r) => Integer -> Integer -> Maybe r -> String +showlatex_exp d n x + | y == 1 = "10^{" ++ show (-n) ++ "}" + | otherwise = printf ("%." ++ show d ++ "f") y ++ "\\cdot 10^{" ++ show (-n) ++ "}" + where + y = case x of + Nothing -> 0 + Just x -> round_to d (10 ** (fromInteger n - x)) + +-- | Either 'show' or 'showlatex', depending on boolean flag. +showf :: (Show a, ShowLaTeX a) => Bool -> a -> String +showf True = showlatex +showf False = show + +-- | Either 'show_exp' or 'showlatex_exp', depending on boolean flag. +showf_exp :: (Show r, RealFrac r, Floating r, PrintfArg r) => Bool -> Integer -> Integer -> Maybe r -> String +showf_exp True = showlatex_exp +showf_exp False = show_exp + +-- | Expand a random seed /g/ into an infinite list of random seeds. +expand_seed :: (RandomGen g) => g -> [g] +expand_seed g = g1 : expand_seed g2 where + (g1,g2) = split g + +-- | Output the given string, right-padded to /n/ characters using spaces. +putStrPad :: Int -> String -> IO() +putStrPad n s = putStr (s ++ replicate (n-l) ' ') + where + l = length s diff --git a/programs/newsynth.hs b/programs/newsynth.hs deleted file mode 100644 --- a/programs/newsynth.hs +++ /dev/null @@ -1,270 +0,0 @@ --- | This module provides a command line interface to the --- decomposition library. - -module Main where - -import Quantum.Synthesis.Newsynth -import Quantum.Synthesis.SymReal -import Quantum.Synthesis.CliffordT -import Quantum.Synthesis.Ring -import Quantum.Synthesis.Matrix -import Quantum.Synthesis.LaTeX - -import CommandLine - --- import other stuff -import Control.Monad -import Data.Time -import System.Console.GetOpt -import System.Environment -import System.Exit -import System.Random -import Text.Printf - --- ---------------------------------------------------------------------- --- * Option processing - --- | A data type to hold values set by command line options. -data Options = Options { - opt_digits :: Double, -- ^ Requested precision in digits (default: 10). - opt_theta :: SymReal, -- ^ Angle to approximate. - opt_hex :: Bool, -- ^ Output operator in hex coding? (default: ASCII). - opt_stats :: Bool, -- ^ Output statistics? - opt_latex :: Bool, -- ^ Additional LaTeX output? - opt_count :: Integer, -- ^ Repetition count for stats (default: 1). - opt_rseed :: Maybe StdGen -- ^ An optional random seed. -} deriving Show - --- | The default options. -defaultOptions :: Options -defaultOptions = Options - { opt_digits = 10, - opt_theta = 0.0, - opt_hex = False, - opt_stats = False, - opt_latex = False, - opt_count = 1, - opt_rseed = Nothing - } - --- | The list of command line options, in the format required by 'getOpt'. -options :: [OptDescr (Options -> IO Options)] -options = - [ Option ['h'] ["help"] (NoArg help) "print usage info and exit", - Option ['d'] ["digits"] (ReqArg digits "") "set precision in decimal digits (default: 10)", - Option ['b'] ["bits"] (ReqArg bits "") "set precision in bits", - Option ['e'] ["epsilon"] (ReqArg epsilon "") "set precision as epsilon (default: 1e-10)", - Option ['x'] ["hex"] (NoArg hex) "output hexadecimal coding (default: ASCII)", - Option ['s'] ["stats"] (NoArg stats) "output statistics", - Option ['l'] ["latex"] (NoArg latex) "additional output in LaTeX format", - Option ['c'] ["count"] (ReqArg count "") "average statistics over runs (default: 1)", - Option ['r'] ["rseed"] (ReqArg rseed "\"\"") "set optional random seed (default: random)" - ] - where - help :: Options -> IO Options - help o = do - usage - exitSuccess - - digits :: String -> Options -> IO Options - digits string o = - case parse_double string of - Just n | n >= 0 -> return o { opt_digits = n } - Just n -> optfail ("Number of digits must not be negative -- " ++ string ++ "\n") - _ -> optfail ("Invalid digits -- " ++ string ++ "\n") - - bits :: String -> Options -> IO Options - bits string o = - case parse_double string of - Just n | n >= 0 -> return o { opt_digits = n * logBase 10 2 } - Just n -> optfail ("Number of bits must not be negative -- " ++ string ++ "\n") - _ -> optfail ("Invalid bits -- " ++ string ++ "\n") - - epsilon :: String -> Options -> IO Options - epsilon string o = - case parse_double string of - Just eps | eps < 1 && eps > 0 -> return o { opt_digits = -logBase 10 eps } - Just n -> optfail ("Epsilon must be between 0 and 1 -- " ++ string ++ "\n") - _ -> optfail ("Invalid epsilon -- " ++ string ++ "\n") - - hex :: Options -> IO Options - hex o = return o { opt_hex = True } - - stats :: Options -> IO Options - stats o = return o { opt_stats = True } - - latex :: Options -> IO Options - latex o = return o { opt_latex = True } - - count :: String -> Options -> IO Options - count string o = - case parse_int string of - Just n | n >= 1 -> return o { opt_count = n } - Just n -> optfail ("Invalid count, must be positive -- " ++ string ++ "\n") - _ -> optfail ("Invalid count -- " ++ string ++ "\n") - - rseed :: String -> Options -> IO Options - rseed string o = - case reads string of - [(g, "")] -> return o { opt_rseed = Just g } - _ -> optfail ("Invalid random seed -- " ++ string ++ "\n") - --- | Process /argv/-style command line options into an 'Options' structure. -dopts :: [String] -> IO Options -dopts argv = do - let (o, args, errs) = getOpt Permute options argv - opts <- foldM (flip id) defaultOptions o - when (errs /= []) $ do - optfail (concat errs) - case args of - [] -> optfail "Missing argument: theta.\n" - [string] -> do - case parse_SymReal string of - Just theta -> return opts { opt_theta = theta } - _ -> optfail ("Invalid theta -- " ++ string ++ "\n") - h1:h2:[] -> optfail ("Too many non-option arguments -- " ++ h1 ++ ", " ++ h2 ++ "\n") - h1:h2:_ -> optfail ("Too many non-option arguments -- " ++ h1 ++ ", " ++ h2 ++ "...\n") - --- | Print usage message to 'stdout'. -usage :: IO () -usage = do - putStr (usageInfo header options) - where header = - "Usage: newsynth [OPTION...] \n" - ++ "Arguments:\n" - ++ " z-rotation angle. May be symbolic, e.g. pi/128\n" - ++ "Options:" - --- ---------------------------------------------------------------------- --- * The main function - --- | Main function: read options, then execute the appropriate tasks. -main :: IO() -main = do - -- Read options. - argv <- getArgs - options <- dopts argv - let digits = opt_digits options - let prec = digits * logBase 2 10 - let theta = opt_theta options - let count = opt_count options - let exponent = ceiling digits - - -- Set random seed. - g <- case opt_rseed options of - Nothing -> newStdGen - Just g -> return g - - -- Expand random seed opt_count times. - let gs = expand_seed count g - - -- Do it for each element of gs. - stats <- sequence $ flip map (zip gs [1..]) $ \(g,n) -> do - - when (count > 1 && (opt_stats options || opt_latex options)) $ do - putStrLn ("Solution " ++ show n ++ ":") - - -- Payload. - t0 <- getCurrentTime - let (m,err,ct) = newsynth_stats prec theta g - gates = to_gates m - if opt_hex options then - printf "%x\n" (convert gates :: Integer) - else - putStrLn (if gates == [] then "I" else convert gates) - t1 <- getCurrentTime - - -- Print optional statistics - let tcount = length $ filter (==T) gates - let secs = diffUTCTime t1 t0 - - when (opt_stats options || opt_latex options) $ do - putStrLn ("Random seed: " ++ show g) - putStrLn ("T-count: " ++ show tcount) - - when (opt_stats options) $ do - putStrLn ("Theta: " ++ show theta) - putStrLn ("Epsilon: " ++ show_exp 10 exponent (Just digits)) - putStrLn ("Matrix: " ++ show m) - putStrLn ("Actual error: " ++ show_exp 10 exponent err) - putStrLn ("Runtime: " ++ show secs) - putStrLn ("Candidates tried: " ++ show ct) - putStrLn ("Time/candidate: " ++ show (secs / fromInteger ct)) - - -- Optional LaTeX output - when (opt_latex options) $ do - putStrLn ("LaTeX Gates: " ++ showlatex gates) - putStrLn ("LaTeX Theta: " ++ showlatex theta) - putStrLn ("LaTeX Epsilon: " ++ showlatex_exp 5 exponent (Just digits)) - putStrLn ("LaTeX Matrix: " ++ showlatex (convert gates :: U2 DOmega)) - putStrLn ("LaTeX Actual error: " ++ showlatex_exp 5 exponent err) - putStrLn ("LaTeX Runtime: " ++ show (round_to 2 secs)) - putStrLn ("LaTeX Candidates tried: " ++ show ct) - putStrLn ("LaTeX Time/candidate: " ++ show (round_to 4 (secs / fromInteger ct))) - - when (count > 1 && (opt_stats options || opt_latex options)) $ do - putStrLn "" - - return (ct, secs) - - -- If count > 1, show summary stats. - when (count > 1) $ do - let (cts, secss) = unzip stats - let ct_total = sum cts - let secs_total = sum secss - - when (opt_stats options || opt_latex options) $ do - putStrLn "Summary:" - putStrLn ("Number of runs: " ++ show count) - putStrLn ("Total runtime: " ++ show secs_total) - - when (opt_stats options) $ do - putStrLn ("Average runtime: " ++ show (secs_total / fromInteger count)) - putStrLn ("Average candidates tried: " ++ show (fromInteger ct_total / fromInteger count :: Double)) - putStrLn ("Average time/candidate: " ++ show (secs_total / fromInteger ct_total)) - - when (opt_latex options) $ do - putStrLn ("LaTeX Average runtime: " ++ show (round_to 2 (secs_total / fromInteger count))) - putStrLn ("LaTeX Average candidates tried: " ++ show (fromInteger ct_total / fromInteger count :: Double)) - putStrLn ("LaTeX Average time/candidate: " ++ show (round_to 4 (secs_total / fromInteger ct_total))) - --- ---------------------------------------------------------------------- --- * Miscellaneous - --- | Round a 'RealFrac' value to the given number of decimals. -round_to :: (RealFrac r) => Integer -> r -> r -round_to n x = fromInteger (round (x * 10^n)) / 10^n - --- | Show the number 10[sup -/x/] in the format 10^(-n) or --- 1.23*10^(-n), with precision /d/ and exponent -/n/. A value of --- 'Nothing' is treated as 0. --- --- For example, display @0.316*10^(-13)@ instead of @10^(-13.5)@. -show_exp :: (Show r, RealFrac r, Floating r, PrintfArg r) => Integer -> Integer -> Maybe r -> String -show_exp d n x - | y == 1 = "10^(" ++ show (-n) ++ ")" - | otherwise = printf ("%." ++ show d ++ "f") y ++ "*10^(" ++ show (-n) ++ ")" - where - y = case x of - Nothing -> 0 - Just x -> round_to d (10 ** (fromInteger n - x)) - --- | Show the number 10[sup -/x/] in the format @10^{-n}@ or --- @1.23\\cdot 10^{-n}@, with precision /d/ and exponent -/n/. A value --- of 'Nothing' is treated as 0. -showlatex_exp :: (Show r, RealFrac r, Floating r, PrintfArg r) => Integer -> Integer -> Maybe r -> String -showlatex_exp d n x - | y == 1 = "10^{" ++ show (-n) ++ "}" - | otherwise = printf ("%." ++ show d ++ "f") y ++ "\\cdot 10^{" ++ show (-n) ++ "}" - where - y = case x of - Nothing -> 0 - Just x -> round_to d (10 ** (fromInteger n - x)) - --- | Expand a random seed /g/ into a list [/g/[sub 1], …, --- /g/[sub /n/]] of /n/ random seeds. This is done in such a way that --- /g/[sub 1] = /g/. -expand_seed :: (RandomGen g) => Integer -> g -> [g] -expand_seed 0 g = [] -expand_seed n g = g:expand_seed (n-1) g' where - (g', _) = split g