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

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

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.

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

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

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

A point in the plane.

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

An operator is a real 2×2-matrix.

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

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.

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₁].

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

The closed unit disk.

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

A closed rectangle with the given dimensions.

Given bounded convex sets A and B, enumerate all solutions u ∈ ℤ[ω] of the 2-dimensional grid problem for A and 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.

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

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.

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

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.

Construct a 2×2-matrix, by rows.

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

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

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.

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.

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.

The determinant of a 2×2-matrix.

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

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

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

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

The bias of a state is ζ - z.

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

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

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

The operator A^k.

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

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

The operator B^k.

The special grid operator K.

The Pauli X operator is a special grid operator.

The Pauli operator Z is a special grid operator.

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.

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

Return S^k.

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

(D, Δ) ⋅ G := (G^†DG, G^•TΔG^•).

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

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

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

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

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

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.

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.

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.

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

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.

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

Apply a linear transformation G to a point p.

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

Apply a special grid operator G to a characteristic function.

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

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

Calculate the bounding box for an ellipse.

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

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

Given an interval, return a slightly bigger one.

The constant λ = 1 + √2.

The constant λ⁻¹ = √2 - 1.

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.

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

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.

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.

The inner product of two points.

Subtract two points.

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