packages feed

newsynth 0.1.1.0 → 0.2

raw patch · 31 files changed

+2327/−706 lines, 31 filesdep +containersdep ~fixedprecnew-component:exe:gridsynthbinary-added

Dependencies added: containers

Dependency ranges changed: fixedprec

Files

ChangeLog view
@@ -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
+ Quantum/Synthesis/Diophantine.hs view
@@ -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))+++
Quantum/Synthesis/EuclideanDomain.hs view
@@ -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
+ Quantum/Synthesis/GridProblems.hs view
@@ -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\". <http://arxiv.org/abs/1403.2975>.++-- ----------------------------------------------------------------------+-- ** 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))+
+ Quantum/Synthesis/GridSynth.hs view
@@ -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\". <http://arxiv.org/abs/1403.2975>.+-- +-- 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)+
Quantum/Synthesis/LaTeX.hs view
@@ -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
Quantum/Synthesis/Matrix.hs view
@@ -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).
Quantum/Synthesis/Newsynth.hs view
@@ -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. <http://arxiv.org/abs/1212.6253>.+-- +-- 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
+ Quantum/Synthesis/QuadraticEquation.hs view
@@ -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
Quantum/Synthesis/Ring.hs view
@@ -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)+
+ Quantum/Synthesis/StepComp.hs view
@@ -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+
+ images/Re.png view

binary file changed (absent → 6028 bytes)

+ images/area.png view

binary file changed (absent → 2343 bytes)

+ images/bz.png view

binary file changed (absent → 1597 bytes)

+ images/ellipse-rectangle.png view

binary file changed (absent → 5850 bytes)

+ images/gridop-A.png view

binary file changed (absent → 1093 bytes)

+ images/gridop-Ai.png view

binary file changed (absent → 1188 bytes)

+ images/gridop-B.png view

binary file changed (absent → 1252 bytes)

+ images/gridop-Bi.png view

binary file changed (absent → 1373 bytes)

+ images/gridop-K.png view

binary file changed (absent → 1465 bytes)

+ images/gridop-R.png view

binary file changed (absent → 1286 bytes)

+ images/gridop-S.png view

binary file changed (absent → 1259 bytes)

+ images/gridop-Si.png view

binary file changed (absent → 1369 bytes)

+ images/gridop-V.png view

binary file changed (absent → 1391 bytes)

+ images/gridop-X.png view

binary file changed (absent → 1089 bytes)

+ images/gridop-Z.png view

binary file changed (absent → 1144 bytes)

+ images/gridop.png view

binary file changed (absent → 1795 bytes)

+ images/skew.png view

binary file changed (absent → 1346 bytes)

newsynth.cabal view
@@ -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\",-    <http://arxiv.org/abs/1212.6253>.+  * "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\", <http://arxiv.org/abs/1403.2975>.   .   * "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
+ programs/gridsynth.hs view
@@ -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 "<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 ['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",+    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 "<n>")      "repeat count for --table mode (default: 50)",+    Option ['r'] ["rseed"]   (ReqArg rseed "\"<s>\"")  "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...] <theta>\n"+            ++ "Arguments:\n"+            ++ " <theta>                   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
− programs/newsynth.hs 
@@ -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 "<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 ['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 "<n>")   "average statistics over <n> runs (default: 1)",-    Option ['r'] ["rseed"]   (ReqArg rseed "\"<s>\"") "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...] <theta>\n"-            ++ "Arguments:\n"-            ++ " <theta>                   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