Gridsynth: Ross-Selinger synthesis

The implementation of the Ross-Selinger (gridsynth) synthesis algorithm follows a structured six-step framework as summarized in the algorithm tab. We can divide them into two distinct processing sections: Phase 1 - Matrix Approximation and Phase 2 - Exact Synthesis.

Matrix Approximation via Geometric Search

Step 1: Mapping the Region

In the first phase, we convert the target rotation angle \theta and precision tolerance \epsilon into a geometric search task. To do so, we map the approximation problem to a two-dimensional grid problem, where we find a ring element u (i) that lies in a region A in the complex plane and (ii) whose \sqrt{2}-conjugate u^{\bullet} lies in a region B. Here A = R_\epsilon is the \epsilon-region, which is formed by the intersection of the unit disk D with the half-plane \vec u \cdot \vec z \geq 1 - \epsilon^2/2, where \vec z = e^{-i\theta/2} and the circular segment is of the width \approx \epsilon^2/2. In contrast, the constraint that acts on the conjugate copy of the plane is simply B=D [1].

Step 2: Grid Projection

As the target matrix must be built from the exact ring entries in \mathbb{Z}[\omega] scaled by 1/\sqrt{2}^{\,k}, the continuous region is intersected with a discrete scaled grid \mathrm{Grid}(B) = \{\, u \in \tfrac{1}{\sqrt{2}^{\,k}}\mathbb{Z}[\omega] : u^{\bullet} \in B \,\}. As the discreteness comes from the conjugate constraint u^{\bullet}\in B and not from truncating a dense set, the grid appears non-square and non-uniformly spaced. Collectively, this boils down the approximation task to locating the candidate lattice points that fall inside the \epsilon-region R_\epsilon. Internally, this region is tracked by a tight, skewed enclosing ellipse with semi-axes a \propto \epsilon and b \propto \epsilon^2.

Step 3: Lattice Reduction & Bounded Search

Searching an elongated, skewed region directly is computationally expensive as the candidate enumerator scans the axis-aligned bounding box of the region, and for a skewed needle, that box is mostly empty space. In the left panel below, the skewed region forces a tall bounding box to scan a thousand lattice points, yet only the thin diagonal of thirty points actually lies in the region. To solve this problem, the algorithm performs a grid transformation by applying a grid operator G. This is a special matrix with determinant \pm 1 built from integer shears and the scaling unit \lambda = 1+\sqrt{2}. Geometrically, the scaling unit acts as a squeeze mapping, stretching the ellipse in one direction while squishing it in the other, whereas the shears shift the coordinate space to counteract the extreme tilt.

Applying G makes the enclosing ellipse upright, and the algorithm repeatedly applies it to drive its uprightness above a threshold of 1/6 [1]. Because G is a lattice automorphism, the solutions in the region are preserved exactly, while the bounding box collapses, i.e., the upright box scans multifold fewer points for the same solutions.

The Step Lemma guarantees each application step shrinks the skew by at least 10\% [1], so this upright form is reached in \mathcal{O}(\log(1/\epsilon)) steps; solutions are then enumerated cheaply and mapped back by G^{-1}. For our example above, we can see it in action:

Exact Synthesis via Algebraic Conversion

Step 4: Diophantine Solving

At the end of the last step, for each grid candidate u, the number theoretic engine forms the following:

\xi = 2^{k} - u^{*} u \;\in\; \mathbb{Z}[\sqrt{2}],

and attempts to solve the Diophantine equation t^{*} t = \xi for t \in \mathbb{Z}[\omega]. The candidates are tested in the order produced by the grid enumeration at a current denominator exponent k, where k is incremented only when no candidate at that level yields a solution. The task of finding t can be mapped to factoring the integer norm n = \xi^{\bullet}\xi. A solution exists iff the prime factors of n satisfy strict algebraic constraints; if a candidate fails these constraints, it lacks a valid representation in \mathbb{Z}[\omega] and the engine skips to testing the next candidate.

Because the number-theoretic properties of the obtained candidates fluctuate unpredictably, this step reduces to a probabilistic factoring iteration. However, since a positive constant fraction of all candidates are mathematically guaranteed to be solvable, a hit typically lands within the first few tries. In practice, this allows the process to behave as a rapid, constant-time arithmetic test rather than a prolonged search bottleneck.

Step 5: Unitary Matrix Construction

Once both u and t are determined, the isolated ring elements are assembled into a mathematically exact unitary over \mathbb{D}[\omega] [1]:

U = \frac{1}{\sqrt{2}^{\,k}} \begin{pmatrix} u & -t^{*} \\ t & u^{*} \end{pmatrix},

which is a zero-error algebraic description of the approximating operator that preserves unitarity.

Step 6: Sequence Generation

To turn the exact unitary U into a Clifford+T gate sequence, it is mapped to its Matsumoto-Amano normal form:

(T \mid I) (HT\mid SHT)^{*} \mathcal{C}.

The number theoretic engine on this step works on the SO(3) representation and repeatedly checks divisibility by the ring factor \delta = 1+\omega, peeling off syllables (HT or SHT) and driving a recursive least denominator exponent (LDE) reduction [2]. Each iteration strips a power of \delta until the matrix reaches \mathrm{LDE} = 0 and collapses to a Clifford sequence \mathcal{C}, making the resulting T-count at most 2k.

This tree diagram illustrates the traceback pathway as a deterministic finite automaton: the final step transforms U to its SO(3) form and reduces its denominator exponent down to a base computational identity node based on the k-parities [3].

References

[1] N. J. Ross and Peter Selinger, Optimal ancilla-free Clifford+T approximation of z-rotations. arXiv preprint arXiv:1403.2975, (2014).

[2] Ken Matsumoto and Kazuyuki Amano, Representation of Quantum Circuits with Clifford and π/8 Gates. arXiv preprint arXiv:0806.3834, (2008).

[3] Brett Giles and Peter Selinger, Remarks on Matsumoto and Amano's normal form for single-qubit Clifford+T operators. arXiv preprint arXiv:1312.6584, (2013).