Quantum.Synthesis.GridProblems

This module provides functions for solving one- and two-dimensional grid problems.

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 α^• ∈ 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

gridpoints :: (RootTwoRing r, Fractional r, Floor r, Ord r) => (r, r) -> (r, r) -> [ZRootTwo]

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_parity :: (RootTwoRing r, Fractional r, Floor r, Ord r) => Integer -> (r, r) -> (r, r) -> [ZRootTwo]

Like gridpoints, but only produce solutions a + b√2 where a has the same parity as the given integer.

Randomized solutions

gridpoint_random :: (RootTwoRing r, Fractional r, Floor r, Ord r, RandomGen g) => (r, r) -> (r, r) -> g -> Maybe ZRootTwo

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_parity :: (RootTwoRing r, Fractional r, Floor r, Ord r, RandomGen g) => Integer -> (r, r) -> (r, r) -> g -> Maybe ZRootTwo

Like gridpoint_random, but only produce solutions a + b√2 where a has the same parity as the given integer.

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^k such that α ∈ A and α^• ∈ B.

gridpoints_scaled :: (RootTwoRing r, Fractional r, Floor r, Ord r) => (r, r) -> (r, r) -> Integer -> [DRootTwo]

Given intervals A = [x₀, x₁] and B = [y₀, y₁], and an integer k ≥ 0, output all solutions α ∈ ℤ[√2] / √2^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_parity :: (RootHalfRing r, Fractional r, Floor r, Ord r) => DRootTwo -> (r, r) -> (r, r) -> Integer -> [DRootTwo]

Like gridpoints_scaled, but assume k ≥ 1, take an additional parameter β ∈ ℤ[√2] / √2^k, and return only those α such that β − α ∈ ℤ[√2] / √2^k-1.

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^• ∈ 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.

type Point r = (r, r)

A point in the plane.

point_fromDRootTwo :: RootHalfRing r => Point DRootTwo -> Point r

Convert a point with coordinates in DRootTwo to a point with coordinates in any RootHalfRing.

type Operator a = Matrix Two Two a

An operator is a real 2×2-matrix.

data Ellipse r

An ellipse is given by an operator D and a center p; the ellipse in this case is

A = { v | (v-p)^† D (v-p) ≤ 1}.

Constructors

Ellipse (Operator r) (Point r)

Instances

Show r => Show (Ellipse r)

Methods

showsPrec :: Int -> Ellipse r -> ShowS

show :: Ellipse r -> String

showList :: [Ellipse r] -> ShowS

type CharFun = Point DRootTwo -> Bool

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 LineIntersector r = Point DRootTwo -> Point DRootTwo -> Maybe (r, r)

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 (an approximation of) t₀ and t₁ such that p(t) ∈ A iff t ∈ [t₀, t₁].

data ConvexSet r

A compact convex set is given by a bounding ellipse, a characteristic function, and a line intersector.

Constructors

ConvexSet (Ellipse r) CharFun (LineIntersector r)

Instances

Show r => Show (ConvexSet r)

Methods

showsPrec :: Int -> ConvexSet r -> ShowS

show :: ConvexSet r -> String

showList :: [ConvexSet r] -> ShowS

Specific convex sets

unitdisk :: (Fractional r, Ord r, RootHalfRing r, Quadratic QRootTwo r) => ConvexSet r

The closed unit disk.

disk :: (Fractional r, Ord r, RootHalfRing r, Quadratic QRootTwo r) => DRootTwo -> ConvexSet r

A closed disk of radius √s, centered at the origin. Assume s > 0.

rectangle :: (Fractional r, Ord r, RootHalfRing r) => (r, r) -> (r, r) -> ConvexSet r

A closed rectangle with the given dimensions.

General solutions

gridpoints2 :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> [DOmega]

Given bounded convex sets A and B, enumerate all solutions u ∈ ℤ[ω] of the 2-dimensional grid problem for A and B.

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^k such that u ∈ A and u^• ∈ B.

gridpoints2_scaled :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> Integer -> [DOmega]

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_with_gridop :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> Operator DRootTwo -> Integer -> [DOmega]

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_increasing :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> [(Integer, [DOmega])]

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. 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_with_gridop :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r, Floor r) => ConvexSet r -> ConvexSet r -> Operator DRootTwo -> [(Integer, [DOmega])]

Like gridpoints2_increasing, except that instead of performing a precomputation, we input the desired grid operator. It must make the two given sets upright.

Implementation details

One-dimensional grid problems

gridpoints_internal :: (RootTwoRing r, Fractional r, Floor r, Ord r) => (r, r) -> (r, r) -> [ZRootTwo]

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.

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

toOperator :: ((a, a), (a, a)) -> Operator a

Construct a 2×2-matrix, by rows.

fromOperator :: Operator a -> ((a, a), (a, a))

Extract the entries of a 2×2-matrix, by rows.

op_fromDRootTwo :: RootHalfRing r => Operator DRootTwo -> Operator r

Convert an operator with entries in DRootTwo to an operator with entries in any RootHalfRing.

operator_from_bz :: (RootTwoRing a, Floating a) => a -> a -> Operator a

The (b,z)-representation of a positive operator with determinant 1 is

where b, z ∈ ℝ and e > 0 with e² = b² + 1. Create such an operator from parameters b and z.

operator_to_bz :: (Fractional a, Real a, RootTwoRing a) => Operator a -> (a, Double)

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_bl2z :: (Floating a, Real a) => Operator a -> (a, a)

A version of operator_to_bz that returns (b, λ^2z) instead of (b, z). This is a critical optimization, as this function is called often, and logarithms are relatively expensive to compute.

det :: Ring a => Operator a -> a

The determinant of a 2×2-matrix.

operator_skew :: Ring a => Operator a -> a

Compute the skew of a positive operator of determinant 1. We define the skew of a positive definite real operator D to be

uprightness :: Floating a => Operator a -> a

Compute the uprightness of a positive operator D.

The uprightness of D is the ratio of the area of the ellipse E = {v | v^†Dv ≤ 1} to the area of its bounding box R. It is given by

States

type OperatorPair a = (Operator a, Operator a)

A state is a pair (D, Δ) of real positive definite matrices of determinant 1. It encodes a pair of ellipses.

skew :: Ring a => OperatorPair a -> a

The skew of a state is the sum of the skews of the two operators.

bias :: (Fractional a, Real a, RootTwoRing a) => OperatorPair a -> Double

The bias of a state is ζ - z.

Grid operators

Consider the set ℤ[ω] ⊆ ℂ. In identifying ℂ with ℝ², we can alternatively identify ℤ[ω] with the set of all vectors (x, y)^† ∈ ℝ² 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

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 D[√2], the automorphism x ↦ x^•, 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.

opR :: RootHalfRing r => Operator r

The special grid operator R: a clockwise rotation by 45°.

opA :: Ring r => Operator r

The special grid operator A: a clockwise shearing with offset 2, parallel to the x-axis.

opA_inv :: Ring r => Operator r

The special grid operator A⁻¹: a counterclockwise shearing with offset 2, parallel to the x-axis.

opA_power :: RootTwoRing r => Integer -> Operator r

The operator A^k.

opB :: RootTwoRing r => Operator r

The special grid operator B: a clockwise shearing with offset √2, parallel to the x-axis.

opB_inv :: RootTwoRing r => Operator r

The special grid operator B⁻¹: a counterclockwise shearing with offset √2, parallel to the x-axis.

opB_power :: RootTwoRing r => Integer -> Operator r

The operator B^k.

opK :: RootHalfRing r => Operator r

The special grid operator K.

opX :: Ring r => Operator r

The Pauli X operator is a special grid operator.

opZ :: Ring r => Operator r

The Pauli operator Z is a special grid operator.

opS :: RootTwoRing r => Operator r

The special grid operator S: a scaling by λ = 1+√2 in the x-direction, and by λ⁻¹ = -1+√2 in the y-direction.

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_inv :: RootTwoRing r => Operator r

The special grid operator S⁻¹, the inverse of opS.

opS_power :: RootTwoRing r => Integer -> Operator r

Return S^k.

Action of grid operators on states

action :: (RealFrac r, RootHalfRing r, Adjoint r) => (Operator r, Operator r) -> Operator DRootTwo -> (Operator r, Operator r)

Compute the right action of a grid operator G on a state (D, Δ). This is defined as:

(D, Δ) ⋅ G := (G^†DG, G^•TΔ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.

shift_sigma :: RootTwoRing a => Integer -> Operator a -> Operator a

Given an operator D, compute σ^kDσ^k.

shift_tau :: RootTwoRing a => Integer -> Operator a -> Operator a

Given an operator Δ, compute τ^kΔτ^k.

shift_state :: RootTwoRing a => Integer -> OperatorPair a -> OperatorPair a

Compute the k-shift of a state (D,Δ).

Skew reduction

lemma_A :: (RealFrac r, RootTwoRing r, Floating r) => r -> r -> Integer

An implementation of the A-Lemma. Given z and ζ, compute the integer m such that the operator A^m reduces the skew.

lemma_B :: (RealFrac r, RootTwoRing r, Floating r) => r -> r -> Integer

An implementation of the B-Lemma. Given z and ζ, compute the integer m such that the operator B^m reduces the skew.

lemma_A_l2 :: (RealFrac r, RootTwoRing r, Floating r) => r -> r -> Integer

A version of lemma_A that inputs λ^2z instead of z and λ^2ζ instead of ζ. Compute the constant m such that the operator A^m reduces the skew.

lemma_B_l2 :: (RealFrac r, RootTwoRing r, Floating r) => r -> r -> Integer

A version of lemma_B that inputs λ^2z instead of z and λ^2ζ instead of ζ. Compute the constant m such that the operator B^m reduces the skew.

step_lemma :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r) => OperatorPair r -> Maybe (Operator DRootTwo)

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.

reduction :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r) => OperatorPair r -> Operator DRootTwo

Repeatedly apply the Step Lemma to the given state, until the skew is 15 or less.

to_upright :: (RealFrac r, Floating r, Ord r, RootTwoRing r, RootHalfRing r, Adjoint r) => OperatorPair r -> Operator DRootTwo

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_sets :: (Adjoint r, RootHalfRing r, RealFrac r, Floating r) => ConvexSet r -> ConvexSet r -> Operator DRootTwo

Given a pair of convex sets, return a grid operator G making both sets upright.

Action of special grid operators on convex sets

point_transform :: Ring r => Operator r -> Point r -> Point r

Apply a linear transformation G to a point p.

ellipse_transform :: (Ring r, Adjoint r) => Operator r -> Ellipse r -> Ellipse r

Apply a special linear transformation G to an ellipse A. This results in the new ellipse G(A) = { G(z) | z ∈ A }.

charfun_transform :: Operator DRootTwo -> CharFun -> CharFun

Apply a special grid operator G to a characteristic function.

lineintersector_transform :: Ring r => Operator DRootTwo -> LineIntersector r -> LineIntersector r

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 }.

convex_transform :: (Ring r, Adjoint r, RootHalfRing r) => Operator DRootTwo -> ConvexSet r -> ConvexSet r

Apply a special linear transformation G to a convex set A. This results in the new convex set G(A) = { G(z) | z ∈ A }.

Bounding boxes

boundingbox_ellipse :: Floating r => Ellipse r -> ((r, r), (r, r))

Calculate the bounding box for an ellipse.

boundingbox :: Floating r => ConvexSet r -> ((r, r), (r, r))

Calculate a bounding box for a convex set. Returns ((x₀, x₁), (y₀, y₁)).

Auxiliary functions

within :: Ord a => a -> (a, a) -> Bool

We write x `within` (a,b) for a ≤ x ≤ b, or equivalently, x ∈ [a, b].

fatten_interval :: Fractional a => (a, a) -> (a, a)

Given an interval, return a slightly bigger one.

lambda :: RootTwoRing r => r

The constant λ = 1 + √2.

lambda_inv :: RootTwoRing r => r

The constant λ⁻¹ = √2 - 1.

lambdapower :: RootTwoRing r => Integer -> r

Return λ^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.

signpower :: Num r => Integer -> r

Return (-1)^k, where k ∈ ℤ.

floorlog :: (Fractional b, Ord b) => b -> b -> (Integer, b)

Given positive numbers b and x, return (n, r) such that

• x = r bⁿ and
• 1 ≤ r < b.

In other words, let n = ⌊log_b x⌋ and r = x b⁻ⁿ. 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.

logBase_double :: (Fractional a, Real a) => a -> a -> Double

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.

iprod :: Num r => Point r -> Point r -> r

The inner product of two points.

point_sub :: Num r => Point r -> Point r -> Point r

Subtract two points.

special_inverse :: Ring r => Operator r -> Operator r

Calculute the inverse of an operator of determinant 1. Note: this does not work correctly for operators whose determinant is not 1.