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