newsynth 0.2.0.1 → 0.3
raw patch · 8 files changed
+434/−145 lines, 8 files
Files
- ChangeLog +7/−0
- Quantum/Synthesis/GridProblems.hs +155/−73
- Quantum/Synthesis/GridSynth.hs +161/−22
- Quantum/Synthesis/QuadraticEquation.hs +33/−9
- Quantum/Synthesis/Ring.hs +22/−10
- Quantum/Synthesis/Ring/FixedPrec.hs +6/−0
- newsynth.cabal +1/−1
- programs/gridsynth.hs +49/−30
ChangeLog view
@@ -1,6 +1,13 @@ ChangeLog +v0.3 2015/05/15+ (2015/04/15) NJR, PS1 - added a new --phase option for+ approximation up to a phase.+ (2015/04/15) NJR, PS1 - various bug fixes.+ v0.2.0.1 2014/10/08+ (2014/10/07) PS1 - added Applicative and Functor instances to+ silence compiler warnings. (2014/10/06) PS1 - updated dependencies for compatibility with base 4.7 and random 1.1. (2014/09/05) PS1 - fixed a bug where the actual T-count was output
Quantum/Synthesis/GridProblems.hs view
@@ -1,3 +1,5 @@+{-# LANGUAGE FlexibleContexts #-}+ -- | This module provides functions for solving one- and -- two-dimensional grid problems. @@ -7,6 +9,8 @@ import Quantum.Synthesis.Matrix import Quantum.Synthesis.QuadraticEquation +import Control.Monad+import Data.Maybe import System.Random -- ----------------------------------------------------------------------@@ -33,42 +37,25 @@ -- | 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 α.+-- 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 ] +gridpoints (x0, x1) (y0, y1) = do+ [ beta | beta' <- gridpoints_internal (x0', x1') (y0', y1'),+ let beta = beta' + alpha,+ test beta ] 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+ a = floor_of (x0 + y0) `div` 2+ b = floor_of (roottwo * (x0 - y0)) `div` 4+ alpha = RootTwo a b+ xoff = fromZRootTwo alpha+ yoff = fromZRootTwo (adj2 alpha)+ x0' = x0 - xoff + x1' = x1 - xoff + y0' = y0 - yoff + y1' = y1 - yoff + test x = fromZRootTwo x `within` (x0, x1) && fromZRootTwo (adj2 x) `within` (y0, y1)+ -- | 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]@@ -213,13 +200,10 @@ -- 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))+-- /v/ + /tw/, where /v/ and /w/ ≠ 0 are vectors. Given /v/ and /w/,+-- the line intersector returns (an approximation of) /t/₀ and /t/₁+-- such that /p/(/t/) ∈ /A/ iff /t/ ∈ [/t/₀, /t/₁].+type LineIntersector r = (Point DRootTwo -> Point DRootTwo -> Maybe (r, r)) -- | A compact convex set is given by a bounding ellipse, a -- characteristic function, and a line intersector.@@ -232,22 +216,39 @@ -- ** Specific convex sets -- | The closed unit disk.-unitdisk :: (Fractional r, Ord r, RootHalfRing r, Quadratic r) => ConvexSet r+unitdisk :: (Fractional r, Ord r, RootHalfRing r, Quadratic QRootTwo r) => ConvexSet r unitdisk = ConvexSet ell tst int where ell = Ellipse 1 (0,0) int p v- | q == Nothing = (1, 0)- | otherwise = (t0, t1)+ | q == Nothing = Nothing+ | otherwise = Just (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)+ q = quadratic (fromDRootTwo a :: QRootTwo) (fromDRootTwo b) (fromDRootTwo c) Just (t0, t1) = q tst (x,y) = x^2 + y^2 <= 1 +-- | A closed disk of radius √/s/, centered at the origin. Assume /s/ > 0.+disk :: (Fractional r, Ord r, RootHalfRing r, Quadratic QRootTwo r) => DRootTwo -> ConvexSet r+disk s = ConvexSet ell tst int where+ ell = Ellipse (1/fromDRootTwo s `scalarmult` 1) (0,0)++ int p v+ | q == Nothing = Nothing+ | otherwise = Just (t0, t1)+ where+ a = iprod v v+ b = 2 * iprod v p+ c = iprod p p - s+ q = quadratic (fromDRootTwo a :: QRootTwo) (fromDRootTwo b) (fromDRootTwo c)+ Just (t0, t1) = q+ + tst (x,y) = x^2 + y^2 <= s+ -- | 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@@ -259,11 +260,11 @@ 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)+ | vx == 0 && px `within` (x0, x1) = Just (min t0y t1y, max t0y t1y)+ | vx == 0 = Nothing+ | vy == 0 && py `within` (y0, y1) = Just (min t0x t1x, max t0x t1x)+ | vy == 0 = Nothing+ | otherwise = Just (t0, t1) where (px, py) = p (vx, vy) = v@@ -308,15 +309,20 @@ -- 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+gridpoints2_scaled setA setB = gridpoints2_scaled_with_gridop setA setB opG where+ opG = to_upright_sets setA setB++-- | Like 'gridpoints2_scaled', except that instead of performing a+-- precomputation, we input the desired grid operator. It must make+-- the two given sets upright.+gridpoints2_scaled_with_gridop :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> Operator DRootTwo -> Integer -> [DOmega]+gridpoints2_scaled_with_gridop setA setB opG = 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@@ -335,22 +341,26 @@ 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)+ let xs = gridpoints_scaled (x0A, x1A+lambda) (x0B, x1B+lambda) (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)+ let iA = intA' (x0, beta') (dx, 0)+ let iB = intB' (x0_bul, beta'_bul) (dx_bul, 0)+ guard (isJust iA)+ guard (isJust iB)+ let Just (t0A, t1A) = iA+ let Just (t0B, t1B) = iB -- 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)))- + let dtA = 10 / max 10 (2^k * (t1B - t0B))+ let dtB = 10 / max 10 (2^k * (t1A - t0A))+ -- Enumerate the solutions in the x-coordinate (ensuring correct -- parity to make sure it's a grid point) -- @@ -374,25 +384,86 @@ -- | 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 +-- in order of increasing /k/. The results are returned in the form+-- +-- > [ (0, l0), (1, l1), (2, l2), ... ],+-- +-- where /l0/ is a list of solutions for /k/=0, /l1/ is a list of+-- solutions for /k/=1, and so on.+gridpoints2_increasing :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> [(Integer, [DOmega])]+gridpoints2_increasing setA setB = gridpoints2_increasing_with_gridop setA setB opG where- solutions_fun = gridpoints2_scaled setA setB- solutions = solutions_fun 0 ++ additional_solutions 1- additional_solutions k = exact_solutions k ++ additional_solutions (k+1)+ opG = to_upright_sets setA setB++-- | Like 'gridpoints2_increasing', except that instead of performing+-- a precomputation, we input the desired grid operator. It must make+-- the two given sets upright.+gridpoints2_increasing_with_gridop :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> Operator DRootTwo -> [(Integer, [DOmega])]+gridpoints2_increasing_with_gridop setA setB opG = solutions + where+ solutions_fun = gridpoints2_scaled_with_gridop setA setB opG+ solutions = (0, solutions_fun 0) : additional_solutions 1+ additional_solutions k = (k, exact_solutions k) : additional_solutions (k+1) exact_solutions k = [ z | z <- solutions_fun k, denomexp z == k ] -- ---------------------------------------------------------------------- -- * Implementation details +-- ----------------------------------------------------------------------+-- ** One-dimensional grid problems++-- | Similar to 'gridpoints', except:+-- +-- 1. Assume that /x0/ and /y0/ are not too far from the origin (say,+-- between -10 and 10). This is to avoid problems with numeric+-- instability when /x0/ and /y0/ are much larger than /dx/ and /dy/,+-- respectively. /y0/ are not too far from the origin.+-- +-- 2. The function potentially returns some non-solutions, so the+-- caller should test for accuracy.+gridpoints_internal :: (RootTwoRing r, Fractional r, Floor r, Ord r) => (r, r) -> (r, r) -> [ZRootTwo]+gridpoints_internal (x0, x1) (y0, y1)+ | dy <= 0 && dx > 0 = + map adj2 $ gridpoints_internal (y0, y1) (x0, x1)+ | dy >= lambda && even n =+ map (lambda_inv_n *) $ gridpoints_internal (lambda_n*x0, lambda_n*x1) (lambda_bul_n*y0, lambda_bul_n*y1)+ | dy >= lambda && odd n =+ map (lambda_inv_n *) $ gridpoints_internal (lambda_n*x0, lambda_n*x1) (lambda_bul_n*y1, lambda_bul_n*y0)+ | dy > 0 && dy < 1 && even n = + map (lambda_m *) $ gridpoints_internal (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_internal (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] ] + 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)++-- ----------------------------------------------------------------------+-- ** Two-dimensional grid problems+ -- $ 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\". <http://arxiv.org/abs/1403.2975>. -- ------------------------------------------------------------------------- ** Positive operators and ellipses+-- *** Positive operators and ellipses -- | Construct a 2×2-matrix, by rows. toOperator :: ((a, a), (a, a)) -> Operator a@@ -467,7 +538,7 @@ ((a, b), (c, d)) = fromOperator m -- ------------------------------------------------------------------------- ** States+-- *** States -- | A state is a pair (/D/, Δ) of real positive definite matrices of -- determinant 1. It encodes a pair of ellipses.@@ -485,7 +556,7 @@ (beta, zeta) = operator_to_bz matB -- ------------------------------------------------------------------------- ** Grid operators+-- *** Grid operators -- $ Consider the set ℤ[ω] ⊆ ℂ. In identifying ℂ with ℝ², we can -- alternatively identify ℤ[ω] with the set of all vectors (/x/,@@ -611,7 +682,7 @@ | otherwise = opS_inv^(-k) -- ------------------------------------------------------------------------- ** Action of grid operators on states+-- *** Action of grid operators on states -- | Compute the right action of a grid operator /G/ on a state (/D/, -- Δ). This is defined as:@@ -625,7 +696,7 @@ g4 = op_fromDRootTwo (adj2 g) -- ------------------------------------------------------------------------- ** Shifts+-- *** Shifts -- $ A shift is not quite the application of a grid operator, because -- the shifts σ and τ actually involve a square root of λ. However,@@ -646,7 +717,7 @@ shift_state k (d,delta) = (shift_sigma k d, shift_tau k delta) -- ------------------------------------------------------------------------- ** Skew reduction+-- *** 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.@@ -779,9 +850,20 @@ a' = a `scalardiv` (sqrt (det a)) b' = b `scalardiv` (sqrt (det b)) opG = reduction (a',b')++-- | Given a pair of convex sets, return a grid operator /G/ making+-- both sets upright.+to_upright_sets :: (Adjoint r, RootHalfRing r, RealFrac r, Floating r) => ConvexSet r -> ConvexSet r -> Operator DRootTwo+to_upright_sets setA setB = opG where+ ConvexSet ellA tstA intA = setA+ ConvexSet ellB tstB intB = setB+ Ellipse matA ctrA = ellA+ Ellipse matB ctrB = ellB+ + opG = to_upright (matA, matB) -- ------------------------------------------------------------------------- ** Action of special grid operators on convex sets+-- *** 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@@ -827,7 +909,7 @@ tstA' = charfun_transform opG tstA -- ------------------------------------------------------------------------- ** Bounding boxes+-- *** Bounding boxes -- | Calculate the bounding box for an ellipse. boundingbox_ellipse :: (Floating r) => Ellipse r -> ((r, r), (r, r))
Quantum/Synthesis/GridSynth.hs view
@@ -1,4 +1,5 @@ {-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE FlexibleContexts #-} -- | This module implements the approximate single-qubit synthesis -- algorithm of@@ -24,6 +25,7 @@ import Quantum.Synthesis.QuadraticEquation import System.Random+import Data.Function import Data.Number.FixedPrec -- ----------------------------------------------------------------------@@ -67,15 +69,6 @@ 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;@@ -84,6 +77,23 @@ data DStatus = Success | Fail | Timeout deriving (Eq, Show) +-- | 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, Integer, DStatus)])+gridsynth_stats g prec theta effort = dynamic_fixedprec2 digits f prec theta where+ digits = ceiling (15 + 2.5 * prec * logBase 10 2) -- heuristic formula!+ f prec theta = gridsynth_internal g prec theta effort+ +-- | A version of 'gridsynth_stats' that returns the optimal operator+-- /up to a global phase/. (The default behavior is to return the+-- optimal operator exactly).+gridsynth_phase_stats :: (RandomGen g) => g -> Double -> SymReal -> Int -> (U2 DOmega, Maybe Double, [(DOmega, Integer, DStatus)])+gridsynth_phase_stats g prec theta effort = dynamic_fixedprec2 digits f prec theta where+ digits = ceiling (15 + 2.5 * prec * logBase 10 2) -- heuristic formula!+ f prec theta = gridsynth_phase_internal g prec theta effort+ -- ---------------------------------------------------------------------- -- * Implementation details @@ -97,7 +107,7 @@ -- -- \[center [image Re.png]] -epsilon_region :: (Floating r, Ord r, RootHalfRing r, Quadratic r) => r -> r -> ConvexSet r+epsilon_region :: (Floating r, Ord r, RootHalfRing r, Quadratic QRootTwo r) => r -> r -> ConvexSet r epsilon_region epsilon theta = ConvexSet ell tst int where -- A bounding ellipse for the ε-region.@@ -111,16 +121,16 @@ -- 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)+ | q == Nothing = Nothing+ | vz == 0 && rhs <= 0 = Just (t0, t1)+ | vz == 0 && otherwise = Nothing+ | vz > 0 = Just (max t0 t2, t1)+ | otherwise = Just (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)+ q = quadratic (fromDRootTwo a :: QRootTwo) (fromDRootTwo b) (fromDRootTwo c) Just (t0, t1) = q -- solve (p + tv) * z >= d@@ -137,6 +147,48 @@ d = 1 - epsilon^2/2 z = (zx, zy) +-- | The ε-/region/, scaled by an additional factor of √/s/, where /s/+-- > 0. The center of scaling is the origin.+epsilon_region_scaled :: (Floating r, Ord r, RootHalfRing r, Quadratic QRootTwo r) => DRootTwo -> r -> r -> ConvexSet r+epsilon_region_scaled s epsilon theta = ConvexSet ell tst int where+ + -- A bounding ellipse for the ε-region.+ ell = Ellipse mat ctr+ ctr = (rd*zx, rd*zy)+ mat = bmat * mmat * special_inverse bmat+ mmat = toOperator ((ev1, 0), (0, ev2))+ bmat = toOperator ((zx, -zy), (zy, zx))+ ev1 = 4 * (1 / epsilon)^4 / fromDRootTwo s+ ev2 = (1 / epsilon)^2 / fromDRootTwo s+ + -- A line intersector for the ε-region.+ int p v+ | q == Nothing = Nothing+ | vz == 0 && rhs <= 0 = Just (t0, t1)+ | vz == 0 && otherwise = Nothing+ | vz > 0 = Just (max t0 t2, t1)+ | otherwise = Just (t0, min t1 t2)+ where+ a = iprod v v+ b = 2 * iprod v p+ c = iprod p p - s+ q = quadratic (fromDRootTwo a :: QRootTwo) (fromDRootTwo b) (fromDRootTwo c)+ Just (t0, t1) = q+ + -- solve (p + tv) * z >= d * r+ -- equivalently, t * vz >= d * r - pz+ vz = iprod (point_fromDRootTwo v) z+ rhs = rd - iprod (point_fromDRootTwo p) z+ t2 = rhs / vz++ -- The characteristic function of the ε-region.+ tst (x, y) = x^2 + y^2 <= s && zx * fromDRootTwo x + zy * fromDRootTwo y >= rd+ + zx = cos (-theta/2)+ zy = sin (-theta/2)+ rd = (1 - epsilon^2/2) * sqrt (fromDRootTwo s)+ z = (zx, zy)+ -- ---------------------------------------------------------------------- -- ** Main algorithm implementation @@ -154,18 +206,23 @@ -- 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 :: forall r g.(RootHalfRing r, Ord r, Floating r, Adjoint r, Floor r, RealFrac r, Quadratic QRootTwo r, RandomGen g) => g -> r -> r -> Int -> (U2 DOmega, Maybe Double, [(DOmega, Integer, 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+ raw_candidates = gridpoints2_increasing region unitdisk+ candidates = [ (u, t) | (k, us) <- raw_candidates,+ let t = tcount k,+ u <- us ] (uU, log_err, candidate_info) = first_solvable [] g candidates + tcount k = if k > 0 then 2*k - 2 else 0+ 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+ first_solvable candidate_info g ((u, tcount) : us) = case answer_t of+ Just (Just t) -> let (uU, log_err) = with_successful_candidate u t in (uU, log_err, ((u, tcount, Success) : candidate_info))+ Just Nothing -> first_solvable ((u, tcount, Fail) : candidate_info) g2 us+ Nothing -> first_solvable ((u, tcount, Timeout) : candidate_info) g2 us where (g1, g2) = split g xi = real (1 - adj u * u)@@ -183,3 +240,85 @@ uU_fixed = matrix_map fromDOmega uU zrot_fixed = zrot (theta :: r) +-- | The internal implementation of the ellipse-based approximate+-- synthesis algorithm, up to a phase. The parameters are the same as+-- for 'gridsynth_internal'.+gridsynth_phase_internal :: forall r g.(RootHalfRing r, Ord r, Floating r, Adjoint r, Floor r, RealFrac r, Quadratic QRootTwo r, Quadratic r r, RandomGen g) => g -> r -> r -> Int -> (U2 DOmega, Maybe Double, [(DOmega, Integer, DStatus)])+gridsynth_phase_internal g prec theta effort = (uU, log_err, candidate_info) where+ epsilon = 2 ** (-prec)+ region0 = epsilon_region epsilon theta+ disk0 = unitdisk+ region1 = epsilon_region_scaled (2 + roottwo) epsilon theta+ disk1 = disk (2 - roottwo)+ opG = to_upright_sets region0 disk0+ raw_candidates0 = gridpoints2_increasing_with_gridop region0 disk0 opG+ raw_candidates1 = gridpoints2_increasing_with_gridop region1 disk1 opG+ candidates0 = [ (t, 0, us) | (k, us) <- raw_candidates0,+ let t = tcount k ]+ candidates1 = [ (t, 1, us') | (k, us) <- raw_candidates1,+ let t = 1 + tcount k,+ let us' = [ u * delta_inv | u <- us ] ]+ merged = mergeBy (compare `on` first) candidates0 candidates1+ candidates = [ (t, ph, u) | (t, ph, us) <- merged, u <- us ]+ (uU, log_err, candidate_info) = first_solvable [] g candidates+ + fabs (Cplx a b) = sqrt(a^2 + b^2)++ tcount k = if k > 0 then 2*k - 2 else 0++ first_solvable candidate_info g [] = error "gridsynth: internal error: finite list of candidates?"+ first_solvable candidate_info g ((tcount, phase, u) : us) = case answer_t of+ Just (Just t) -> + let (uU, log_err) = with_successful_candidate u t phase in + (uU, log_err, ((u, tcount, Success) : candidate_info))+ Just Nothing -> first_solvable ((u, tcount, Fail) : candidate_info) g2 us+ Nothing -> first_solvable ((u, tcount, 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 0 = (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)+ + with_successful_candidate u t 1 = (uU, log_err) where+ uU | denomexp (u + t) < denomexp (u + omega * t)+ = matrix2x2 (u, -(adj t) * omega_inv) (t, adj u * omega_inv)+ | otherwise+ = matrix2x2 (u, -(adj t)) (t * omega_inv, adj u * omega_inv)+ log_err + | err <= 0 = Nothing+ | otherwise = Just (logBase_double 0.5 err)+ err = sqrt (real (hs_sqnorm (sqrt_omega `scalarmult` uU_fixed - zrot_fixed)) / 2)+ uU_fixed = matrix_map fromDOmega uU+ zrot_fixed = zrot (theta :: r)+ sqrt_omega = Cplx (cos (pi/8)) (sin (pi/8)) + omega_inv = omega^7++ delta_inv = roothalf * (omega - i)++-- ----------------------------------------------------------------------+-- * Auxiliary functions++-- | Merge the elements of two lists in increasing order, assuming+-- that each of the lists is already sorted. The first argument is a+-- comparison function for elements.+mergeBy :: (a -> a -> Ordering) -> [a] -> [a] -> [a]+mergeBy c [] l2 = l2+mergeBy c l1 [] = l1+mergeBy c (h1:t1) (h2:t2)+ | c h1 h2 == LT = h1:(mergeBy c t1 (h2:t2))+ | otherwise = h2:(mergeBy c (h1:t1) t2)++-- | Return the first component of a triple.+first :: (a,b,c) -> a+first (a,b,c) = a
Quantum/Synthesis/QuadraticEquation.hs view
@@ -1,3 +1,6 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+ -- | This module provides a type class 'Quadratic', for solving -- quadratic equations. @@ -7,6 +10,7 @@ import Data.Number.FixedPrec import Quantum.Synthesis.Ring+import Quantum.Synthesis.ToReal -- | This type class provides a primitive method for solving quadratic -- equations. For many floating-point or fixed-precision@@ -14,14 +18,14 @@ -- 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+class Quadratic t a where+ -- | '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+ -- the coefficients /a/, /b/, and /c/ can be taken to be of an exact -- type; therefore instances have the opportunity to work with -- infinite precision.- quadratic :: QRootTwo -> QRootTwo -> QRootTwo -> Maybe (a, a)+ quadratic :: t -> t -> t -> Maybe (a, a) -- ---------------------------------------------------------------------- -- FixedPrec instance@@ -34,7 +38,7 @@ -- * 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 :: (Fractional t, Floor t, Ord t) => t -> t -> Maybe (Integer, Integer) int_quadratic b c | radix < 0 = Nothing | otherwise = Just (t0, t1)@@ -67,8 +71,8 @@ -- /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 +quadratic_fixedprec :: (Fractional t, Floor t, Ord t, Precision e) => t -> t -> t -> Maybe (FixedPrec e, FixedPrec e)+quadratic_fixedprec a b c | False = Just (r, r) | otherwise = do (x0, x1) <- int_quadratic b' c'@@ -82,5 +86,25 @@ c' = prec'^2 * c/a q = int_quadratic b' c' -instance (Precision e) => Quadratic (FixedPrec e) where- quadratic = qroottwo_quadratic_fixedprec+instance (Fractional t, Floor t, Ord t, Precision e) => Quadratic t (FixedPrec e) where+ quadratic = quadratic_fixedprec++-- ----------------------------------------------------------------------+-- Double instance++instance (ToReal t) => Quadratic t Double where+ quadratic a' b' c'+ | radix < 0 = Nothing+ | b >= 0 = Just (t1, t2)+ | otherwise = Just (t1', t2')+ where+ radix = b^2 - 4*a*c+ s1 = -b - sqrt radix+ s2 = -b + sqrt radix+ t1 = s1 / (2*a)+ t2 = (2*c) / s1+ t1' = (2*c) / s2+ t2' = s2 / (2*a)+ a = to_real a'+ b = to_real b'+ c = to_real c'
Quantum/Synthesis/Ring.hs view
@@ -79,9 +79,20 @@ -- ** Rings with √2 -- | A type class for rings that contain √2.+-- +-- Minimal complete definition: 'roottwo'. The default definition of+-- 'fromZRootTwo' uses the expression @x+roottwo*y@. However, this can+-- give potentially bad round-off errors for fixed-precision types,+-- where the expression @roottwo*y@ can be vastly inaccurate if @y@ is+-- large. For such rings, one should provide a custom definition. class (Ring a) => RootTwoRing a where -- | The square root of 2. roottwo :: a++ -- | The unique ring homomorphism from ℤ[√2] to any ring containing+ -- √2. This exists because ℤ[√2] is the free such ring.+ fromZRootTwo :: (RootTwoRing a) => ZRootTwo -> a+ fromZRootTwo (RootTwo x y) = fromInteger x + roottwo * fromInteger y instance RootTwoRing Double where roottwo = sqrt 2@@ -96,9 +107,20 @@ -- ** Rings with 1\/√2 -- | A type class for rings that contain 1\/√2.+-- +-- Minimal complete definition: 'roothalf'. The default definition of+-- 'fromDRootTwo' uses the expression @x+roottwo*y@. However, this can+-- give potentially bad round-off errors for fixed-precision types,+-- where the expression @roottwo*y@ can be vastly inaccurate if @y@ is+-- large. For such rings, one should provide a custom definition. class (HalfRing a, RootTwoRing a) => RootHalfRing a where -- | The square root of ½. roothalf :: a++ -- | The unique ring homomorphism from [bold D][√2] to any ring containing+ -- 1\/√2. This exists because [bold D][√2] = ℤ[1\/√2] is the free such ring.+ fromDRootTwo :: (RootHalfRing a) => DRootTwo -> a+ fromDRootTwo (RootTwo x y) = fromDyadic x + roottwo * fromDyadic y instance RootHalfRing Double where roothalf = sqrt 0.5@@ -448,11 +470,6 @@ -- | The ring ℤ[√2]. type ZRootTwo = RootTwo Integer --- | The unique ring homomorphism from ℤ[√2] to any ring containing--- √2. This exists because ℤ[√2] is the free such ring.-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@@ -479,11 +496,6 @@ -- | The ring [bold D][√2] = ℤ[1\/√2]. type DRootTwo = RootTwo Dyadic---- | The unique ring homomorphism from [bold D][√2] to any ring containing--- 1\/√2. This exists because [bold D][√2] = ℤ[1\/√2] is the free such ring.-fromDRootTwo :: (RootHalfRing a) => DRootTwo -> a-fromDRootTwo (RootTwo x y) = fromDyadic x + roottwo * fromDyadic y -- ---------------------------------------------------------------------- -- ** The field ℚ[√2]
Quantum/Synthesis/Ring/FixedPrec.hs view
@@ -8,9 +8,15 @@ instance Precision e => RootHalfRing (FixedPrec e) where roothalf = sqrt 0.5+ fromDRootTwo (RootTwo x y)+ | y >= 0 = fromDyadic x + sqrt (fromDyadic (2*y^2))+ | otherwise = fromDyadic x - sqrt (fromDyadic (2*y^2)) instance Precision e => RootTwoRing (FixedPrec e) where roottwo = sqrt 2+ fromZRootTwo (RootTwo x y)+ | y >= 0 = fromInteger x + sqrt (fromInteger (2*y^2))+ | otherwise = fromInteger x - sqrt (fromInteger (2*y^2)) instance Precision e => HalfRing (FixedPrec e) where half = 0.5
newsynth.cabal view
@@ -7,7 +7,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change-version: 0.2.0.1+version: 0.3 -- A short (one-line) description of the package. synopsis: Exact and approximate synthesis of quantum circuits
programs/gridsynth.hs view
@@ -32,6 +32,7 @@ data Options = Options { opt_digits :: Maybe Double, -- ^ Requested precision in decimal digits (default: 10). opt_theta :: Maybe SymReal, -- ^ The angle θ to approximate.+ opt_phase :: Bool, -- ^ Decompose up to a global phase? 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?@@ -46,6 +47,7 @@ defaultOptions = Options { opt_digits = Nothing, opt_theta = Nothing,+ opt_phase = False, opt_effort = 25, opt_hex = False, opt_stats = False,@@ -62,6 +64,7 @@ Option ['d'] ["digits"] (ReqArg digits "<n>") "set precision in decimal digits (default: 10)", Option ['b'] ["bits"] (ReqArg bits "<n>") "set precision in bits", Option ['e'] ["epsilon"] (ReqArg epsilon "<n>") "set precision as epsilon (default: 1e-10)",+ Option ['p'] ["phase"] (NoArg phase) "decompose up to a global phase (default: no)", Option ['f'] ["effort"] (ReqArg effort "\"<n>\"") "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",@@ -97,6 +100,9 @@ Just n -> optfail ("Epsilon must be between 0 and 1 -- " ++ string ++ "\n") _ -> optfail ("Invalid epsilon -- " ++ string ++ "\n") + phase :: Options -> IO Options+ phase o = return o { opt_phase = True }+ effort :: String -> Options -> IO Options effort string o = case parse_int string of@@ -188,6 +194,9 @@ let exponent = ceiling digits let l = opt_latex options let effort = opt_effort options+ let gridsynth_fun = case opt_phase options of+ False -> gridsynth_stats+ True -> gridsynth_phase_stats -- Set random seed. g <- case opt_rseed options of@@ -196,8 +205,10 @@ -- Payload. t0 <- getCurrentTime- let (m,err,cinfo) = gridsynth_stats g prec theta effort- gates = to_gates m+ let (m,err,cinfo) = gridsynth_fun g prec theta effort+ gates = case opt_phase options of+ False -> to_gates m+ True -> strip_phases (to_gates m) if opt_hex options then printf "%x\n" (convert gates :: Integer) else if opt_latex options then@@ -209,9 +220,7 @@ -- 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 tlower = last [ tcount | (u, tcount, status) <- cinfo, status /= Fail ] let secs = diffUTCTime t1 t0 let err_d = case err of Nothing -> Nothing@@ -227,9 +236,9 @@ 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)")+ ++ show (length [u | (u, tc, Fail) <- cinfo]) ++ " failed, "+ ++ show (length [u | (u, tc, Timeout) <- cinfo]) ++ " timed out, "+ ++ show (length [u | (u, tc, Success) <- cinfo]) ++ " succeeded)") putStrLn ("Time/candidate: " ++ show (secs / fromIntegral ct)) -- ----------------------------------------------------------------------@@ -240,32 +249,34 @@ -- 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+-- the precision in decimal digits, an amount of effort to spend on+-- factoring, and a boolean flag determining whether we should+-- decompose up to a global phase. 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 -> Bool -> IO (U2 DOmega, [Gate], Maybe Double, Int, Integer, Integer, Double)+one_run g theta prec_d effort phase = do+ let gridsynth_fun = case phase of+ False -> gridsynth_stats+ True -> gridsynth_phase_stats 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+ let (op, err, cinfo) = gridsynth_fun g prec theta effort circ = synthesis_u2 op- tcount = length $ filter (==T) circ+ tcount = fromIntegral $ 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+ tlower = last [ tcount | (u, tcount, status) <- cinfo, status /= Fail ] ((u, _), (t, _)) = fromOperator op let err_d = case err of Nothing -> Nothing@@ -277,9 +288,9 @@ 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)")+ ++ show (length [u | (u, tc, Fail) <- cinfo]) ++ " failed, "+ ++ show (length [u | (u, tc, Timeout) <- cinfo]) ++ " timed out, "+ ++ show (length [u | (u, tc, Success) <- cinfo]) ++ " succeeded)") putStrLn ("% Time/candidate: " ++ show (secs / fromIntegral ct)) putStrLn "" hFlush stdout@@ -288,16 +299,17 @@ -- | 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+-- in decimal digits, an amount of effort to spend on factoring, and a+-- flag that determines whether to factor up to a global phase.+many_runs :: (RandomGen g, Show g) => g -> Int -> Double -> Int -> Bool -> IO ()+many_runs g n prec_d effort phase = 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+ one_run g' theta prec_d effort phase -- Output the LaTeX of one row of the table let (_,_,err,_,tcount,tlower,_) = head results total_time = sum [ t | (_,_,_,_,_,_,t) <- results ]@@ -344,6 +356,7 @@ Nothing -> [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 500, 1000] Just d -> [d] let effort = opt_effort options+ let phase = opt_phase options -- Set random seed. g <- case opt_rseed options of@@ -360,8 +373,8 @@ (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+ one_run g1 theta prec_d effort phase+ many_runs g2 count prec_d effort phase -- ---------------------------------------------------------------------- -- * Miscellaneous@@ -416,3 +429,9 @@ putStrPad n s = putStr (s ++ replicate (n-l) ' ') where l = length s++-- | Strip global phase gates from a word.+strip_phases :: [Gate] -> [Gate]+strip_phases [] = []+strip_phases (W:xs) = strip_phases xs+strip_phases (x:xs) = x : strip_phases xs