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)

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

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)

Specific convex sets

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

The closed unit disk.

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_increasing :: (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 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.

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

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.

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.