# Synthesizing Toffoli-optimal quantum circuits for arbitrary multi-qubit
unitaries

Priyanka Mukhopadhyay<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper we study the Clifford+Toffoli universal fault-tolerant gate set.
We introduce a generating set in order to represent any unitary implementable
by this gate set and with this we derive a bound on the Toffoli-count of
arbitrary multi-qubit unitaries. We analyse the channel representation of the
generating set elements, with the help of which we infer
$|\mathcal{J}_{n}^{Tof}|<|\mathcal{J}_{n}^{T}|$, where $\mathcal{J}_{n}^{Tof}$
and $\mathcal{J}_{n}^{T}$ are the set of unitaries exactly implementable by
the Clifford+Toffoli and Clifford+T gate set, respectively. We develop
Toffoli-count optimal synthesis algorithms for both approximately and exactly
implementable multi-qubit unitaries. With the help of these we prove
$|\mathcal{J}_{n}^{Tof}|=|\mathcal{J}_{n}^{CS}|$, where $\mathcal{J}_{n}^{CS}$
is the set of unitaries exactly implementable by the Clifford+CS gate set.

##  1 Introduction

Quantum circuits is a popular model for describing and implementing quantum
algorithms, similar to the case of classical computation. These consist of a
series of elementary operations or gates belonging to a universal set, which
usually consists of Clifford group gates and at least one non-Clifford gate
[1, 2]. Example of some universal gate sets include Clifford+T, V-basis,
Clifford+CS, Toffoli+H. Often, choice of implementing gate set is determined
by the underlying hardware technology. An important component of any quantum
computer compilation process is quantum circuit synthesis and optimization.
Depending upon the application the goal is often to optimize the number of a
specific gate. For example, optimization of multi-qubit gates like CNOT is
more relevant for NISQ era, when their implementations are more error-prone
than the single qubit gates. Optimization of non-Clifford gates are more
relevant for the fault-tolerant regime because most error-correction schemes
implement Clifford gates transversally, allowing the logical operations to be
performed precisely and with time proportional to the physical gate time. The
non-Clifford gates, however, require large ancilla factories and additional
operations like gate teleportation and state distillation [3, 4]. These are
less accurate procedures which require both additional time and space compared
to a single physical gate. Thus the optimization of multi-qubit non-Clifford
gate like Toffoli is important for both the NISQ as well as fault-tolerant
regime. It is also worth noting that the minimum number of non-Clifford gates
required to implement certain unitaries is a quantifier of difficulty in many
algorithms that try to classically simulate qunatum computation [5, 6].

The Solovay-Kitaev algorithm [7, 8] guarantees that given an $n$-qubit unitary
$W$, we can generate a circuit with a ”discrete finite” universal gate set,
like Clifford+T, V-basis, Clifford+CS, Clifford+Toffoli, such that the unitary
$U$ implemented by the circuit is at most a certain distance from $W$. In
fact, in quantum computation a set of gates is said to be universal if any
quantum operation can be approximated to arbitrary accuracy by a quantum
circuit involving only these gates [2]. A unitary is called exactly
implementable by a gate set if there exists a quantum circuit with these
gates, that implements it (up to some global phase). Otherwise, it is
approximately implementable. Accordingly, a synthesis algorithm can be (a)
exact when $U=e^{i\phi}W$ ($\phi$ is the global phase); or (b) approximate
when $d(U,W)\leq\epsilon$ for some $\epsilon>0$. $d(.)$ is a distance metric
and for quantum circuit synthesis two popular metrics are the operator norm
[9, 10] and global phase invariant distance [11, 12, 13, 14, 15]. In this
paper we use the global phase invariant distance because it ignores the global
phase and hence avoids unnecessarily long approximating sequences that
achieves a specific global phase. This distance is composable and inequalities
relating it to operator norm can be found in [14]. It has also been used to
synthesize unitaries in other models like topological quantum computation [16,
17].

In this paper we consider the Clifford+Toffoli gate set, where the Toffoli
gate (TOF) is defined as follows.

$\displaystyle\text{TOF}=\begin{bmatrix}1&0&0&0&0&0&0&0\\\ 0&1&0&0&0&0&0&0\\\
0&0&1&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&1&0&0\\\
0&0&0&0&0&0&0&1\\\ 0&0&0&0&0&0&1&0\\\ \end{bmatrix}$ (1)

This gate set is universal for quantum computing, in fact universality has
been shown for a more restricted gate set i.e. the Toffoli+H [18, 19, 20]. The
Toffoli+H quantum circuits have been studied in the context of diagrammatic
calculi [21], path-sums [22] and quantum logic [23]. Fault-tolerant
implementations of the Toffoli gate has been shown in [24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34]. Toffoli gates occur frequently while synthesizing
circuits for Hamiltonian simulation [35, 36] and many other applications where
multi-controlled unitaries are used [37, 38, 39, 40]. Algorithms have been
developed to synthesize Toffoli+H circuits for exactly implementable
unitaries, but with no constraint on optimality of the non-Clifford Toffoli
gate [41]. Also, some of these algorithms only work for 3 qubit unitaries and
none work for approximately implementable unitaries.

###  1.1 Our contributions

In the following points we discuss briefly the main results in our paper and
with these we also try to shed light on the motivations for studying the
mentioned problems related to this gate set.

We already mentioned the importance of optimizing the Toffoli gate both for
the NISQ as well as fault-tolerant era. So one aim is to have some
representation or formalism to express the unitaries implementable by this
gate set and also to characterize them. With this we hope to compare the
different universal gate sets, specifically the size of the set of unitaries,
exactly implementable by these gate sets. Such questions are very relevant for
resource estimations of various algorithms. Due to various reasons like the
fidelity of gates, ease of implementation, etc, different hardware platforms
may prefer different universal gate sets. The ability to compare the relative
number of different non-Clifford gates required to implement a certain unitary
with different universal gate sets, will help in making wiser decisions about
which gate set to use [42, 43]. Also, in literature some popular unitaries
have implementations only in a specific gate set. We can then have a more
convenient estimate on the non-Clifford gate count simply by replacing each
such gate with the new gate set. Another important goal is to use these
formalisms in order to design Toffoli-optimal synthesis algorithms for
arbitrary multi-qubit unitaries. We develop both provable and heuristic
algorithms. The former type of algorithms have a rigorous proof of optimality
and the claimed complexity. For the latter type both the claimed optimality
and complexity depends on some unproven conjectures. The need of developing
heuristic algorithms usually arise due to the lack of efficient provable
algorithms. Then, depending on some observations some conjectures are coined,
which if/when proven, will prove the stated claims.

Algorithm |  Type of algorithm |  Type of unitary | Time complexity | Space complexity
---|---|---|---|---
APPROX-TOF-OPT |  Provable |  Approximately implementable | $O\left(n^{2m_{\epsilon}}2^{(6n-12)m_{\epsilon}}\right)$ | $O\left(n^{2}2^{6n-12}\right)$
Nested MITM |  Provable |  Exactly implementable | $O\left(n^{2(c-1)\left\lceil\frac{m}{c}\right\rceil}2^{(6n-12)(c-1)\left\lceil\frac{m}{c}\right\rceil}\right)$ | $O\left(n^{2\left\lceil\frac{m}{c}\right\rceil}2^{(6n-12)\left\lceil\frac{m}{c}\right\rceil}\right)$
EXACT-TOF-OPT |  Heuristic |  Exactly implementable | $\text{poly}(n^{2}2^{6n-12},m)$ | $\text{poly}(n^{2}2^{6n-12},m)$
Table 1: Summary of space and time complexity of the different algorithms
developed in this paper. Here $m$ and $m_{\epsilon}$ are the Toffoli-count and
$\epsilon$-Toffoli-count of the input $n$-qubit unitary ($n>0$), respectively.
$c\geq 2$ is the level of nesting.

  1. 1.

We define a generating set $\mathcal{G}_{Tof}$ such that any unitary exactly
implementable by the Clifford+Toffoli gate set can be written as product of
unitaries from this gate set and a Clifford (Section 3.1). We show that
$|\mathcal{G}_{Tof}|\in O\left(n^{2}2^{6n-12}\right)$.

  2. 2.

With the help of $\mathcal{G}_{Tof}$ we also derive a lower bound of
$\Omega\left(\log_{4}\left(\frac{1}{\epsilon}\right)\right)$ on the
$\epsilon$-Toffoli count of arbitrary multi-qubit unitaries (Section 3.2). For
a unitary $W$, its $\epsilon$-Toffoli count is the minimum number of Toffoli
gates required to implement any unitary that is within $\epsilon$ distance of
$W$.

Most of the previous works have proven bounds on the total gate count or
1-qubit non-Clifford gate count like T [13, 10] and V [44, 11, 12, 9]. Most of
these are upper bounds and work for 1-qubit gates, specifically z-rotation
gates. The bound on T-count derived in [13] is empirical, meaning it is
interpolated from experimental data and there are no rigorous analysis. In
fact, it is not clear how to extend the number-theoretic arguments used in
[11, 12, 9, 10] for multi-qubit gates. To the best of our knowledge, this is
the first paper that derives bound on the count of a multi-qubit non-Clifford
gate, while implementing arbitrary multi-qubit unitary. An upper bound on the
CS-count of 2-qubit unitaries, has been shown in [45].

For the special case of exactly implementable unitaries, an upper bound on the
Toffoli-count of $n$-qubit unitaries has been derived in [41]. Our lower bound
for this special case is quadratically smaller (in $n$) than this upper bound.

  3. 3.

We then derive the channel representation of unitaries in this set (Section
3.3). This representation was introduced in the context of Clifford+T gate set
in [46], and then extensively studied in [47, 48]. Its main use was in the
design of efficient algorithms for exactly implementable unitaries. In this
paper we have found another interesting aspect. It can also be used to
characterize unitaries that can be exactly implemented by these gate sets. For
example, we are able to derive interesting conclusions, like the impossibility
of the exact implementation of the T gate. While, we know that the Toffoli
gate can be implemented exactly with 7 T gates [46, 47] or 4 T gates [28] if
we allow extra ancillae and classical resources like measurement. Such
inferences are not possible with the results derived in earlier works [20,
41], that did not work with this representation.

This implies that
$\left|\mathcal{J}_{n}^{Tof}\right|<\left|\mathcal{J}_{n}^{T}\right|$, where
$\mathcal{J}_{n}^{Tof}$ and $\mathcal{J}_{n}^{T}$ are the set of unitaries
exactly implementable by the Clifford+Toffoli and Clifford+T gate set,
respectively. Using similar arguments we conclude that the V-basis gates
cannot be implemented exactly by the Clifford+Toffoli gate set. In fact,
Toffoli also cannot be implemented exactly by the V-basis gate set.

It is worth mentioning here that one of the advantages of using channel
representation in order to derive such exact implementability results is Fact
3.9, which ensures that a unitary is the channel representation of a Clifford
operator if and only if it has exactly one $\pm 1$ in each row and column.
This implies that to check implementability of a group of the form
Clifford+non-Clifford we need to focus on the channel representation of the
non-Clifford operators only. This is not possible if we work with the unitary
matrices only, as done in [20, 41, 49].

  4. 4.

The channel representation also inherits the advantages already established in
[47], for designing algorithms. First, channel representation of
$\mathcal{G}_{Tof}$ helps us develop Toffoli-count-optimal synthesis
algorithms, using the nested meet-in-the-middle (MITM) algorithm, introduced
in [47]. This algorithm was developed as a recursive generalization of the
MITM algorithm of [46]. The space and time complexity of our algorithm
(Section 4.2.1) is
$O\left(n^{2\lceil\frac{m}{c}\rceil}2^{(6n-12)\left\lceil\frac{m}{c}\right\rceil}\right)$
and
$O\left(n^{2(c-1)\lceil\frac{m}{c}\rceil}2^{(6n-12)(c-1)\left\lceil\frac{m}{c}\right\rceil}\right)$,
respectively. Here $m=\mathcal{T}^{of}(U)$ is the Toffoli-count of the input
unitary $U$ and $c\geq 2$ is the level of nesting used. Second, we define the
smallest denominator exponent ($\text{sde}_{2}$) with respect to 2 and discuss
how it can be used to design exponentially faster heuristic algorithms
(Section 4.2.2). The time and space complexity of our algorithm is
$\text{poly}(n^{2}2^{6n-12},\mathcal{T}^{of}(U))$. Third, we can develop
clever data structures, such that no floating point operations are required at
any step of our algorithms for exactly implementable unitaries. Fourth, we
design a fast multiplication algorithm to multiply a unitary with the channel
representation of a generating set element, which basically boils down to some
row additions involving part of the matrices. This is especially useful when
we have to perform many such multiplications.

To the best of our knowledge, these are the first Toffoli-count-optimal-
synthesis algorithms for exactly implementable unitaries. With the help of
these algorithms we prove that CS can be implemented optimally with 3 Toffoli
gates. In a separate paper we show that Toffoli can be optimally implemented
with 3 CS gates. This implies that
$\left|\mathcal{J}_{n}^{Tof}\right|=\left|\mathcal{J}_{n}^{CS}\right|$, where
$\mathcal{J}_{n}^{Tof}$ and $\mathcal{J}_{n}^{CS}$ are the set of unitaries
exactly implementable by the Clifford+Toffoli and Clifford+CS gate set,
respectively.

  5. 5.

We design a provable Toffoli-count-optimal synthesis algorithm for
approximately implementable unitaries (Section 4.1), analogous to the one in
[15] for Clifford+T. The time complexity of this algorithm is
$O\left(n^{2\mathcal{T}_{\epsilon}^{of}(W)}2^{(6n-12)\mathcal{T}_{\epsilon}^{of}(W)}\right)$,
where $\mathcal{T}_{\epsilon}^{of}(W)$ is the $\epsilon$-Toffoli-count of the
input unitary $W$. The space complexity is $O\left(n^{2}2^{6n-12}\right)$. To
the best of our knowledge, this is the first algorithm that synthesizes
Toffoli-count-optimal circuits for approximately implementable unitaries. We
have summarized the complexities of the various algorithms developed in this
paper in Table 1.

###  1.2 Relevant work

Much work has been done to characterize the unitaries exactly implementable by
the various universal gate sets [20] like Clifford+T [49, 46, 47], V-basis
[44], Clifford+CS [45] (2 qubits only), Toffoli+H [41]. Extensive work has
been done to synthesize arbitrary multi-qubit unitary [7, 50, 8, 51, 52, 53],
without optimality constraint on any particular gate. Some synthesis
algorithms have been developed for Toffoli+H [41] gate set, but these work
only for exactly implementable unitaries on 3 qubits, without any optimality
constraint. In [45] a CS-count-optimal synthesis algorithm has been developed
for exactly implementable 2-qubit unitaries. More work has been done for the
synthesis of V-count-optimal circuits [11, 12, 9], but these work only for
1-qubit unitaries, specifically for the z-rotations. Extensive work has been
done for T-count-optimal synthesis of exactly implementable multi-qubit
unitaries [54, 46, 47] and approximately implementable 1-qubit unitaries [55,
13, 10]. The algorithm in [15] is the only one that returns T-count-optimal
circuits for arbitrary multi-qubit unitaries. The algebraic procedures
introduced in this work circumvents the limitations of the number-theoretic
techniques in other optimal-synthesis papers, like [11, 12, 9, 13, 10], which
exploit the structure of 1-qubit unitaries and hence it is not known how to
generalize these for multi-qubit unitaries. In fact, they often guarantee the
optimal count for specific unitaries like z-rotations. Work has also been done
for T-depth-optimal synthesis [56, 48].

So far we only discussed about optimal-synthesis algorithms where the input is
an $n$-qubit unitary and output is a circuit. The complexity of these
algorithms cannot avoid the exponential dependence on $n$. If additional
optimality constraints are imposed then complexity increases even further, and
often it also exponentially depends on other parameters, for example the
minimum non-Clifford gate count. Heuristics have been designed that reduces
the dependence on these other parameters to polynomial [47]. However, with
exponential dependence on $n$, these algorithms become intractable on a
personal computer. So significant amount of work has been done to develop re-
synthesis algorithms where the input is usually an already synthesized circuit
and the output is a circuit with reduce gate count. Due to relaxed constraint
and more input information, the complexity of these algorithms are polynomial
in input size. However, in literature this reported complexity deos not
include the cost of generating the input circuit which itself is exponential
in $n$ [57]. Usually, optimal-synthesis algorithms are used to synthesize
smaller unitaries that are more frequent in quantum algorithms, for example
rotations. One can do piece-wise optimal synthesis of a larger circuit and
then estimate the non-Clifford gate cost.

##  2 Preliminaries

We denote the $n\times n$ identity matrix by $\mathbb{I}_{n}$ or $\mathbb{I}$
if dimension is clear from the context. We denote the $(i,j)^{th}$ entry of a
matrix $M$ by $M[i,j]$. Here we note that we denote the commutator bracket by
$[A,B]=AB-BA$, where $A,B$ are operators. Though we use the same kind of
brackets, the meaning should be clear from the context. We denote the $i^{th}$
row of $M$ by $M[i,.]$. We denote the set of $n$-qubit unitaries by
$\mathcal{U}_{n}$. The size of an $n$-qubit unitary is $N\times N$ where
$N=2^{n}$. The qubits on which a gate acts is mentioned in the subscript with
brackets. For example, $X_{(q)}$ implies an X gate acting on qubit $q$. For
multi-qubit controlled gate the control qubit(s) is (are) mentioned first,
followed by the target qubit(s), separated by semi-colon. For example,
$\text{CNOT}_{(i;j)}$ denotes CNOT gate controlled on qubit $i$ and target on
qubit $j$ and $\text{TOF}_{(i,j;k)}$ denotes the Toffoli gate with controls on
qubits $i,j$ and target on qubit $k$. For symmetric multi-qubit gates like CS,
where the unitary does not change if we interchange the control and target
qubit, we replace the semi-colon with a comma. For convenience, we skip the
subscript, when it is clear from the context. We have given detail description
about the n-qubit Pauli operators ($\mathcal{P}_{n}$) and the Clifford group
($\mathcal{C}_{n}$) in Appendix A.1. The Clifford+Toffoli group is generated
by $\mathcal{C}_{n}$ and the Toffoli gate (Equation 1).

A unitary $U$ is exactly implementable by the Clifford+Toffoli gate set if
there exists an implementing circuit (up to some global phase) consisting of
these gates, else it is approximately implementable by this gate set. We
denote the set of $n$-qubit exactly implementable (by Clifford+Toffoli)
unitaries by $\mathcal{J}_{n}^{Tof}$. Specifically, we say $W$ is
$\epsilon$-approximately implementable by the Clifford+Toffoli gate set if
there exists an exactly implementable unitary $U$ such that
$d(U,W)\leq\epsilon$. The Solovay-Kitaev algorithm [7, 8] guarantees that any
unitary is $\epsilon$-approximately implementable, for arbitrary precision
$\epsilon\geq 0$. In this paper we use the following distance measure
$d(.,.)$, which has been used in previous works like [51, 13, 15] (qubit based
computing), [16, 17] (topological quantum computing).

######  Definition 2.1 (Global phase invariant distance).

Given two unitaries $U,W\in\mathcal{U}_{n}$, we define the global phase
invariant distance between them as follows.

$\displaystyle
d(U,W)=\sqrt{1-\frac{\left|\text{Tr}\left(U^{\dagger}W\right)\right|}{N}}$

Like the operator norm, this distance is composable and inequalities relating
these two distance metrics have been derived in [14].

###  2.1 Toffoli-count of circuits and unitaries

#### Toffoli-count of circuits

The _Toffoli-count of a circuit_ is the number of Toffoli-gates in it.

#### Toffoli-count of exactly implementable unitaries

The _Toffoli-count of an exactly implementable unitary_ $U$, denoted by
$\mathcal{T}^{of}(U)$, is the minimum number of Toffoli-gates required to
implement it (up to a global phase) with a Clifford+Toffoli circuit.

####  $\epsilon$-Toffoli-count of approximately implementable unitaries

Let $W\in\mathcal{U}_{n}$ be an approximately implementable unitary. The
_$\epsilon$ -Toffoli-count_ of $W$, denoted by
$\mathcal{T}_{\epsilon}^{of}(W)$, is the minimum number of Toffoli gates
required to implement any unitary that is within distance $\epsilon$ of $W$.
In other words $\mathcal{T}_{\epsilon}^{of}(W)=\mathcal{T}^{of}(U)$, where
$U\in\mathcal{J}_{n}^{Tof}$, $d(U,W)\leq\epsilon$ and
$\mathcal{T}^{of}(U)\leq\mathcal{T}^{of}(U^{\prime})$ for any exactly
implementable unitary $U^{\prime}$ within distance $\epsilon$ of $W$. We call
a Toffoli-count-optimal circuit for any such $U$ as the _$\epsilon$ -Toffoli-
count-optimal_ circuit for $W$.

It is not hard to see that the above definition is very general and can be
applied to any unitary $W\in\mathcal{U}_{n}$, exactly or approximately
implementable. If a unitary is exactly implementable then $\epsilon=0$.

##  3 Results and Discussions

###  3.1 Generating set

Any unitary $U$, exactly implementable by Clifford+Toffoli can be written as
alternative product of a Clifford and a Toffoli gate. We can alternatively
express $U$ as product of unitaries that are Clifford conjugation of Toffoli
gates. We remember that Toffoli is a self-inverse unitary and so we need not
consider its inverse.

$\displaystyle U$ $\displaystyle=$ $\displaystyle
e^{i\phi}C_{m}(\text{TOF})C_{m-1}(\text{TOF})\ldots
C_{1}(\text{TOF})C_{0};\qquad\left[C_{i}\in\mathcal{C}_{n},\phi\in[0,2\pi)\right]$
(2) $\displaystyle=$ $\displaystyle
e^{i\phi}\left(C_{m}(\text{TOF})C_{m}^{\dagger}\right)\left(C_{m}C_{m-1}(\text{TOF})C_{m-1}^{\dagger}C_{m}^{\dagger}\right)\ldots\left(C_{m}\ldots
C_{1}(\text{TOF})C_{1}^{\dagger}\ldots C_{m}^{\dagger}\right)C_{m}\ldots
C_{1}C_{0}$ $\displaystyle=$ $\displaystyle
e^{i\phi}\left(\prod_{j=m}^{1}C_{j}^{\prime}(\text{TOF})C_{j}^{\prime\dagger}\right)C_{0}^{\prime}\qquad[C_{j}^{\prime}\in\mathcal{C}_{n};j=0,\ldots,m]$

Let $\text{TOF}_{(i,j;k)}$ be a Toffoli gate with controls on qubits $i,j$ and
target on qubit $k$. The following Pauli expansion can be obtained by
observing that any $n$-qubit unitary $U$ can be expressed in the Pauli basis
as $U=\sum_{i}\alpha_{i}P_{i}$, where $P_{i}\in\mathcal{P}_{n}$ and
$\alpha_{i}=\frac{1}{2^{n}}\text{Tr}(UP_{i})$.

$\displaystyle\text{TOF}_{(i,j;k)}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(\text{X}_{(i)}+\text{Z}_{(j)}+\text{Z}_{(k)}-\text{X}_{(i)}\text{Z}_{(j)}-\text{Z}_{(j)}\text{Z}_{(k)}-\text{Z}_{(k)}\text{X}_{(i)}+\text{X}_{(i)}\text{Z}_{(j)}\text{Z}_{(k)}\right)$
(3)

If $C\in\mathcal{C}_{n}$ is an $n$-qubit Clifford then,

$\displaystyle C(\text{TOF}_{(i,j;k)})C^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(P_{1}+P_{2}+P_{3}-P_{1}P_{2}-P_{2}P_{3}-P_{3}P_{1}+P_{1}P_{2}P_{3}\right):=G_{P_{1},P_{2},P_{3}}$
(4)

where $P_{1}=C\text{X}_{(i)}C^{\dagger}$, $P_{2}=C\text{Z}_{(j)}C^{\dagger}$
and $P_{3}=C\text{Z}_{(k)}C^{\dagger}$. Since
$[\text{Z}_{(j)},\text{Z}_{(k)}]=[\text{X}_{(i)},\text{Z}_{(j)}]=0$ when
$i\neq j\neq k$, so using Fact A.4 we have
$[P_{1},P_{2}]=[P_{2},P_{3}]=[P_{3},P_{1}]=0$. Also, $P_{1}\neq P_{2}$,
$P_{2}\neq P_{3}$, $P_{3}\neq P_{1}$ and $P_{1},P_{2},P_{3}\neq\mathbb{I}$. In
Appendix B (Lemma B.1) we have given detailed proof about the following
properties of the unitaries $G_{P_{1},P_{2},P_{3}}$.

######  Lemma 3.1.

If $P_{1},P_{2},P_{3}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}\\}$ pair-wise
commutes, then

  1. 1.

$G_{P_{1},P_{2},P_{3}}$ remains same for any permutation of the Paulis;

  2. 2.

$G_{P_{1},-P_{1}P_{2},P_{3}}=G_{P_{1},P_{2},P_{3}}$;

  3. 3.

$G_{P_{1},-P_{2},P_{3}}=G_{P_{1},P_{2},P_{3}}C$, for some
$C\in\mathcal{C}_{n}$;

  4. 4.

$G_{P_{1},P_{2},P_{1}P_{2}}=\mathbb{I}$;

  5. 5.

$G_{P_{1},P_{2},P_{1}P_{2}P_{3}}=G_{P_{1},P_{2},P_{3}}$.

From Equation 2 and using the above lemma we can define the following set,
which we call a generating set (modulo Clifford), because we can express any
exactly implementable unitary (up to a global phase) as product of unitaries
from this set and a Clifford.

$\displaystyle\mathcal{G}_{Tof}$ $\displaystyle=$
$\displaystyle\left\\{G_{P_{1},P_{2},P_{3}}:P_{1},P_{2},P_{3}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}\\};\quad
P_{1}\neq P_{2}\neq
P_{3};\quad[P_{1},P_{2}]=[P_{2},P_{3}]=[P_{3},P_{1}]=0;\right.$ (5)
$\displaystyle\left.P_{3}\neq
P_{1}P_{2};\quad(P_{1},P_{2},P_{3})\equiv(P_{1},\pm
P_{1}P_{2},P_{3})\equiv(P_{1},P_{2},P_{1}P_{2}P_{3})\equiv\pi_{3}(P_{1},P_{2},P_{3});\right.$
$\displaystyle\left.\pi_{3}\text{ is a permutation of }3\text{
elements}.\right\\}$

We use
$(P_{1},P_{2},P_{3})\equiv(P_{1}^{\prime},P_{2}^{\prime},P_{3}^{\prime})$ to
imply that only one pair is included in the set. In Algorithm 1 (Appendix E)
we have given the pseudocode for constructing this set. Specifically, we have
the following result.

######  Theorem 3.2.

Any unitary $U$ that is exactly implementable by the Clifford+Toffoli gate set
can be expressed as,

$\displaystyle
U=e^{i\phi}\left(\prod_{j=m}^{1}G_{P_{1j},P_{2j},P_{3j}}\right)C_{0};\qquad[C_{0}\in\mathcal{C}_{n};\phi\in[0,2\pi)];$

where $G_{P_{1j},P_{2j},P_{3j}}\in\mathcal{G}_{Tof}$, a set defined in
Equation 5.

######  Remark 3.1.

Let us consider what happens if we change the position of the target in the
Toffoli unitary described in Equation 3.

$\displaystyle\text{TOF}_{i,k;j}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(\text{X}_{(j)}+\text{Z}_{(i)}+\text{Z}_{(k)}-\text{Z}_{(i)}\text{X}_{(j)}-\text{Z}_{(k)}\text{X}_{(j)}-\text{Z}_{(i)}\text{Z}_{(k)}+\text{X}_{(j)}\text{Z}_{(i)}\text{Z}_{(k)}\right)$
$\displaystyle\text{TOF}_{j,k;i}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(\text{X}_{(k)}+\text{Z}_{(i)}+\text{Z}_{(j)}-\text{Z}_{(j)}\text{X}_{(k)}-\text{Z}_{(i)}\text{X}_{(k)}-\text{Z}_{(i)}\text{Z}_{(j)}+\text{X}_{(k)}\text{Z}_{(i)}\text{Z}_{(j)}\right)$
(6)

So, changing the controls and targets, simply permutes the position of the
Pauli operators. We can use the following Clifford conjugation equalities in
order to change the position of the operators.

$\displaystyle\text{SWAP}(\mathbb{I}\otimes\text{X})\text{SWAP}=\text{X}\otimes\mathbb{I};\qquad\text{SWAP}(\mathbb{I}\otimes\text{Z})\text{SWAP}=\text{Z}\otimes\mathbb{I};\qquad\text{H}\text{Z}\text{H}=\text{X};\qquad\text{H}\text{X}\text{H}=\text{Z}$
(7)

Thus in the generating set $\mathcal{G}_{Tof}$ we need not cosider Clifford
conjugations of the unitaries in Equation 6, that are obtained by changing the
position of the target of the unitary in Equation 3.

##### Cardinality of $\mathcal{G}_{Tof}$ :

We prove the following bound on $|\mathcal{G}_{Tof}|$. The proof has been
given in Appendix B (Theorem B.3).

######  Theorem 3.3.

$\displaystyle|\mathcal{G}_{Tof}|$ $\displaystyle\leq$
$\displaystyle\frac{1}{384}[64^{n}-52^{n}-17(16^{n}-4^{n})]+\frac{1}{576}[48^{n}-2\cdot
24^{n}]-\frac{1}{36}[12^{n}-2\cdot 6^{n}]+\frac{1}{24}(13^{n}-1)$
$\displaystyle\in$ $\displaystyle O\left(n^{2}64^{n-2}\right)$

In Table 2 we have given the cardinality of the generating set obtained in
practice, for different number of qubits, using Algorithm 1 (Appendix E).

$n$ | 3 | 4
---|---|---
$|\mathcal{G}_{tof}|$ | 129 | 7024
Time | 1.087 sec | 1 hr 44 secs
Table 2: Cardinality of $\mathcal{G}_{tof}$ for different number of qubits
$n$.

##### Circuit construction for $G_{P_{1},P_{2},P_{3}}$ :

We can construct a circuit implementing $G_{P_{1},P_{2},P_{3}}$ by deriving
the conjugating Clifford and determining the control and target qubits of the
Toffoli gate. Given 3 commuting non-identity Paulis $P_{1},P_{2},P_{3}$, we
use the algorithm in [58] in order to derive a conjugating Clifford
$C^{\prime}\in\mathcal{C}_{n}$ such that
$C^{\prime}P_{1}C^{\prime\dagger}=P_{1}^{\prime}$,
$C^{\prime}P_{2}C^{\prime\dagger}=P_{2}^{\prime}$ and
$C^{\prime}P_{3}C^{\prime\dagger}=P_{3}^{\prime}$,, where
$P_{1}^{\prime}=\bigotimes_{j=1}^{n}Q_{j}$,
$P_{2}^{\prime}=\bigotimes_{j=1}^{n}R_{j}$,
$P_{3}^{\prime}=\bigotimes_{j=1}^{n}S_{j}$ and
$Q_{j},R_{j},S_{j}\in\\{\mathbb{I},\text{Z}\\}$. That is, the output of this
algorithm is a triple of Z-operators i.e. $n$-qubit Paulis that are tensor
product of either $\mathbb{I}$ or Z. Then we use the following conjugation
relations,

$\displaystyle\text{SWAP}(\mathbb{I}\otimes\text{Z})\text{SWAP}=\text{Z}\otimes\mathbb{I};\quad\text{CNOT}_{(j;k)}(\mathbb{I}_{(j)}\otimes\text{Z}_{(k)})\text{CNOT}_{(j;k)}=\text{Z}_{(j)}\otimes\text{Z}_{(k)}$
$\displaystyle\text{CNOT}_{(j;k)}(\text{Z}_{(j)}\otimes\mathbb{I}_{(k)})\text{CNOT}_{(j;k)}=\text{Z}_{(j)}\otimes\mathbb{I}_{(k)};$

in order to derive Clifford $C^{\prime\prime}\in\mathcal{C}_{n}$ such that
$C^{\prime\prime}P_{1}^{\prime}C^{\prime\prime\dagger}=\text{Z}_{(a)}$,
$C^{\prime\prime}P_{2}^{\prime}C^{\prime\prime\dagger}=\text{Z}_{(b)}$ and
$C^{\prime\prime}P_{3}^{\prime}C^{\prime\prime\dagger}=\text{Z}_{(c)}$, where
$1\leq a,b,c\leq n$ and $a\neq b\neq c$. Now we use H gate (Equation 7) in
order to apply X in any one of the qubits. Suppose we apply $\text{H}_{(a)}$.
If $C=C^{\prime}C^{\prime\prime}\text{H}_{(a)}$, then we get
$C^{\prime}C^{\prime\prime}\text{H}_{(a)}\text{X}_{(a)}\text{H}_{(a)}C^{\prime\prime\dagger}C^{\prime\dagger}=P_{1}$,
$C^{\prime}C^{\prime\prime}\text{H}_{(a)}\text{Z}_{(b)}\text{H}_{(a)}C^{\prime\prime\dagger}C^{\prime\dagger}=P_{2}$
and
$C^{\prime}C^{\prime\prime}\text{H}_{(a)}\text{Z}_{(c)}\text{H}_{(a)}C^{\prime\prime\dagger}C^{\prime\dagger}=P_{3}$.
Thus a circuit for $G_{P_{1},P_{2},P_{3}}$ consists of $\text{TOF}_{(a,b;c)}$,
conjugated by Clifford $C=C^{\prime}C^{\prime\prime}\text{H}_{(a)}$. From
Remark 3.1 we can infer that interchanging the position of the target and
control qubits will change the conjugating Clifford.

###  3.2 Bound on Toffoli-count of unitaries

In this section we derive lower bound on the Toffoli-count of arbitrary
$n$-qubit unitaries. We know that any unitary can be expanded in the Pauli
basis. Let $W$ be an $n$-qubit unitary and

$\displaystyle
W=\sum_{i=1}^{4^{n}}\overline{\alpha_{i}}P_{i}\qquad[P_{i}\in\mathcal{P}_{n}].$
(8)

##### Not exactly implementable :

First we consider the case when $W$ is strictly approximately implementable.
Let
$U=\left(\prod_{j=1}^{m}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\right)C_{0}e^{i\phi}$
is a unitary such that $W=UE$, for some $C_{0}\in\mathcal{C}_{n}$,
$\phi\in[0,2\pi)$ and unitary $E$. Also $d(U,W)\leq\epsilon$, implying
$|\text{Tr}(E)|\geq N(1-\epsilon^{2})$, where $N=2^{n}$ and
$m=\mathcal{T}_{\epsilon}^{of}(W)$. Let
$\widetilde{U}=\prod_{j=1}^{m}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}$. Then,

$\displaystyle\left|\text{Tr}\left(\widetilde{U}W^{\dagger}\right)\right|=\left|\text{Tr}(EC_{0})\right|$
(9)

and if $C_{0}=\sum_{P\in\mathcal{P}_{n}}r_{P}P$, $M=\left|\left\\{P:r_{P}\neq
0\right\\}\right|$, then from Theorem 3.1 of [15] we have the following.

$\displaystyle\frac{1-\epsilon^{2}}{\sqrt{M}}-\sqrt{M(2\epsilon^{2}-\epsilon^{4})}\leq|\text{Tr}(EC_{0}P^{\prime}/N)|\leq\frac{1}{\sqrt{M}}+\sqrt{M(2\epsilon^{2}-\epsilon^{4})}\quad[\text{if
}\quad r_{P^{\prime}}\neq 0]$ (10) $\displaystyle
0\leq|\text{Tr}(EC_{0}P^{\prime}/N)|\leq\sqrt{M(2\epsilon^{2}-\epsilon^{4})}\quad[\text{if
}\quad r_{P^{\prime}}=0]$ (11)

Here $1\leq M\leq N^{2}$. We use this in order to lower bound
$m=\mathcal{T}_{\epsilon}^{of}(W)$ as a function of $\epsilon$. First we
observe that we can expand $\widetilde{U}$ as in the following lemma, the
proof of which has been given in Appendix C.

######  Lemma 3.4.

If $\widetilde{U}=\prod_{j=1}^{m}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}$, where
$G_{P_{1_{j}},P_{2_{j},P_{3_{j}}}}=a\mathbb{I}+bQ_{j}$,
$Q_{j}=P_{1_{j}}+P_{2_{j}}+P_{3_{j}}-P_{1_{j}}P_{2_{j}}-P_{2_{j}}P_{3_{j}}-P_{3_{j}}P_{1_{j}}+P_{1_{j}}P_{2_{j}}P_{3_{j}}$,
$a=\frac{3}{4}$ and $b=\frac{1}{4}$, then

$\displaystyle\widetilde{U}$ $\displaystyle=$ $\displaystyle
a^{m}\mathbb{I}+a^{m-1}b\left(\sum_{j}Q_{j}\right)+a^{m-2}b^{2}\left(\sum_{j_{1}<j_{2}}Q_{j_{1}}Q_{j_{2}}\right)+a^{m-3}b^{3}\left(\sum_{j_{1}<j_{2}<j_{3}}Q_{j_{1}}Q_{j_{2}}Q_{j_{3}}\right)+\cdots$
$\displaystyle\cdots+ab^{m-1}\left(\sum_{j_{1}<j_{2}<\cdots<j_{m-1}}Q_{j_{1}}Q_{j_{2}}\cdots
Q_{j_{m-1}}\right)+b^{m}\prod_{j=1}^{m}Q_{j}.$

Now from Equations 8 and 9 we have

$\displaystyle|\text{Tr}(EC_{0})|=\left|\text{Tr}\left(\widetilde{U}W^{\dagger}\right)\right|=\left|\sum_{i=1}^{4^{n}}\alpha_{i}\text{Tr}\left(\widetilde{U}P_{i}\right)\right|,$

and thus using triangle inequality and Equation 10, we have

$\displaystyle\sum_{i=1}^{4^{n}}|\alpha_{i}||\text{Tr}(\widetilde{U}P_{i}P^{\prime}/N)|$
$\displaystyle\leq$
$\displaystyle\frac{1}{\sqrt{M}}+\sqrt{M(2\epsilon^{2}-\epsilon^{4})}\leq\frac{1}{\sqrt{M}}+\sqrt{2M}\epsilon$
(12)
$\displaystyle\left||\alpha_{k}||\text{Tr}(\widetilde{U}P_{k}P^{\prime}/N)|-\sum_{i\neq
k}|\alpha_{i}||\text{Tr}(\widetilde{U}P_{i}P^{\prime}/N)|\right|$
$\displaystyle\geq$
$\displaystyle\frac{1-\epsilon^{2}}{\sqrt{M}}-\sqrt{M(2\epsilon^{2}-\epsilon^{4})}\geq\frac{1-\epsilon^{2}}{\sqrt{M}}-\sqrt{2M}\epsilon$
$\displaystyle-\left||\alpha_{k}||\text{Tr}(\widetilde{U}P_{k}P^{\prime}/N)|-\sum_{i\neq
k}|\alpha_{i}||\text{Tr}(\widetilde{U}P_{i}P^{\prime}/N)|\right|$
$\displaystyle\leq$
$\displaystyle\frac{-1+\epsilon^{2}}{\sqrt{M}}+\sqrt{2M}\epsilon.$ (13)

Without loss of generality let us assume that
$|\alpha_{k}||\text{Tr}(\widetilde{U}P_{k}P^{\prime}/N)|-\sum_{i\neq
k}|\alpha_{i}||\text{Tr}(\widetilde{U}P_{i}P^{\prime}/N)|>0$. So, adding the
above inequalities we get

$\displaystyle\sum_{i\neq
k}|\alpha_{i}||\text{Tr}(\widetilde{U}P_{i}P^{\prime}/N)|\leq\frac{\epsilon^{2}}{2\sqrt{M}}+\sqrt{2M}\epsilon.$

Now $\sum_{i\neq
k}|\alpha_{i}||\text{Tr}(\widetilde{U}P_{i}P^{\prime}/N)|\geq|\alpha_{\ell}||\text{Tr}(\widetilde{U}P_{\ell}P^{\prime}/N)|$,
for some $P_{\ell}$ such that the absolute value of coefficients on the right
is non-zero. For $M\geq 1$ and for non-Clifford $\widetilde{U}$, $EC_{0}$
there exists at least one such $P_{\ell}$. So,

$\displaystyle|\alpha_{\ell}||\text{Tr}(\widetilde{U}P_{\ell}P^{\prime}/N)|\leq\frac{\epsilon^{2}}{2\sqrt{M}}+\sqrt{2M}\epsilon.$
(14)

Now we observe that $\text{Tr}(\widetilde{U}P_{\ell}P^{\prime}/N)$ is the
coefficient of $P_{\ell}P^{\prime}$ in the Pauli expansion of $\widetilde{U}$,
as given in Lemma 3.4. Consider group of r-product terms of the form
$a^{m-r}b^{r}\left(\sum_{j_{1}<j_{2}<\cdots<j_{r}}Q_{j_{1}}Q_{j_{2}}\ldots
Q_{j_{r}}\right)$. We remember that each $Q_{j}$ is a linear combination of 7
Paulis and $Q_{j}^{2}=7\mathbb{I}-6Q_{j}$ (Lemma D.1 in Appendix D), which
implies that $Q^{m}$ is a linear combination of 8 Paulis, for any $m\geq 2$.
Thus each summand in a r-product group is a product of at most $r$ Paulis and
no two summands have the same combination of Paulis. We remember that we are
interested in the minimum $m$ for Toffoli-count. The expression of
$\widetilde{U}$ is sum of $m+1$ such groups of $r$-products, where each
combination is distinct from the other. Overall, there can be at least a
constant number of terms that are $P_{\ell}P^{\prime}$. Since $|b|<|a|<1$ so
minimum value of the coefficient can be $c|b|^{m}$, where $c$ is a constant.
Thus from Equation 14,

$\displaystyle|\alpha_{\ell}|c|b|^{m}$ $\displaystyle\leq$
$\displaystyle|\alpha_{\ell}||\text{Tr}(\widetilde{U}P_{\ell}P^{\prime}/N)|\leq\frac{\epsilon^{2}}{\sqrt{2M}}+\sqrt{2M}\epsilon=\epsilon\left(\frac{\epsilon}{2\sqrt{M}}+\sqrt{2M}\right).$
$\displaystyle|b|^{m}$ $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{c|\alpha_{\ell}|}\left(\frac{\epsilon}{2\sqrt{M}}+\sqrt{2M}\right)$
$\displaystyle m\log\frac{1}{|b|}$ $\displaystyle\geq$
$\displaystyle\log\frac{1}{\epsilon}+\log\left(\frac{c|\alpha_{\ell}|}{\frac{\epsilon}{2\sqrt{M}}+\sqrt{2M}}\right)$
$\displaystyle m$ $\displaystyle\geq$
$\displaystyle\log_{\frac{1}{|b|}}\left(\frac{1}{\epsilon}\right)-\log_{\frac{1}{|b|}}\left(\frac{\frac{\epsilon}{2\sqrt{M}}+\sqrt{2M}}{c|\alpha_{\ell}|}\right).$

Since $|b|=\frac{1}{4}$, we have the following asymptotic lower bound, with
respect to $\epsilon$.

######  Theorem 3.5.

Let $W=\sum_{i}\alpha_{i}P_{i}$ is an $n$-qubit unitary, approximately
implementable by the Clifford+Toffoli gate set. Here
$P_{i}\in\mathcal{P}_{n}$. Then for some constant $c$,

$\displaystyle\mathcal{T}_{\epsilon}^{of}\left(W\right)$ $\displaystyle\geq$
$\displaystyle\log_{4}\left(\frac{1}{\epsilon}\right)-\log_{4}\left(\frac{\frac{\epsilon}{2\sqrt{M}}+\sqrt{2M}}{c|\alpha_{max}|}\right)\in\Omega\left(\log_{4}\left(\frac{1}{\epsilon}\right)\right).$

##### Exactly implementable :

Now we consider the case when $W$ is exactly implementable i.e. $E=\mathbb{I}$
and $\epsilon=0$. Then, from Fact A.3 in Appendix A ([59]) or plugging in
$\epsilon=0$ in Equation 10 ,we obtain

$\displaystyle\sum_{i=1}^{4^{n}}|\alpha_{i}||\text{Tr}(\widetilde{U}P_{i}P^{\prime}/N)|$
$\displaystyle\leq$ $\displaystyle\frac{1}{\sqrt{M}}.$

Now,
$|\alpha_{\ell}||\text{Tr}(\widetilde{U}P_{\ell}P^{\prime}/N)|\leq\sum_{i=1}^{4^{n}}|\alpha_{i}||\text{Tr}(\widetilde{U}P_{i}P^{\prime}/N)|$
, where $|\alpha_{\ell}|\neq 0$. As previously argued
$c|b|^{m}\leq|\text{Tr}(\widetilde{U}P_{\ell}P^{\prime}/N)|$, for some
constant $c$. So,

$\displaystyle|\alpha_{\ell}|c|b|^{m}$ $\displaystyle\leq$
$\displaystyle\frac{1}{\sqrt{M}}$
$\displaystyle\left(\frac{1}{|b|}\right)^{m}$ $\displaystyle\geq$
$\displaystyle|\alpha_{\ell}|c\sqrt{M}$ $\displaystyle\implies m$
$\displaystyle\geq$
$\displaystyle\log_{\frac{1}{|b|}}\left(|\alpha_{max}|c\sqrt{M}\right)$ (15)

Thus we have the following result.

######  Theorem 3.6.

Let $W=\sum_{i}\alpha_{i}P_{i}$ is an $n$-qubit unitary, exactly implementable
by the Clifford+Toffoli gate set. Here $P_{i}\in\mathcal{P}_{n}$. Then for
some constant $c$ and $1\leq M\leq 4^{n}$,

$\displaystyle\mathcal{T}^{of}\left(W\right)$ $\displaystyle\geq$
$\displaystyle\log_{4}\left(|\alpha_{max}|c\sqrt{M}\right).$

##### Upper bound :

An upper bound on the Toffoli-count can be derived from the Solovay-Kitaev
theorem [7, 50, 8], which gives an upper bound of
$O\left(\log^{c}\left(\frac{1}{\epsilon}\right)\right)$ on the total number of
gates, where $3\leq c\leq 4$ is a constant. This trivially implies an upper
bound on the non-Clifford gate count, as well. Our lower bound in Theorem 3.5
is polynomially smaller than this upper bound.

For exactly implementable unitaries over the Toffoli-Hadamard gate set, an
upper bound of $O(n^{2}\log(n)k)$ (k is a non-negative integer) on the total
gate count was shown in [41]. This implies an upper bound on the Toffoli-
count, though for a restricted gate set (not Clifford+Toffoli). If we plug in
$M=4^{n}$ then our lower bound in Theorem 3.6 is quadratically smaller than
this upper bound.

####  3.2.1 Illustration : Toffoli-count of rotation gates

As examples, we consider the following rotation gates. In Appendix C.1 we have
explicitly derived the Pauli basis expansions that we show below, for each
unitary.

$\displaystyle R_{z}(\theta)$ $\displaystyle=$ $\displaystyle
e^{-i\frac{\theta}{2}}\begin{bmatrix}1&0\\\
0&e^{i\theta}\end{bmatrix}=e^{-i\frac{\theta}{2}}\left(\alpha_{1}\mathbb{I}+\alpha_{2}\text{Z}\right)$
$\displaystyle\mathbb{I}\otimes\mathbb{I}\otimes R_{z}(\theta)$
$\displaystyle=$ $\displaystyle
e^{-i\frac{\theta}{2}}\left(\alpha_{1}\mathbb{I}\mathbb{I}\mathbb{I}+\alpha_{2}\mathbb{I}\mathbb{I}\text{Z}\right)$
(16)

Here $\alpha_{1}=\frac{1+e^{i\theta}}{2}$,
$\alpha_{2}=\frac{1-e^{i\theta}}{2}$ and so

$\displaystyle|\alpha_{1}|$ $\displaystyle=$
$\displaystyle\left|\frac{1}{2}\left(1+\cos\theta+i\sin\theta\right)\right|=\frac{1}{2}\sqrt{((1+\cos\theta)^{2}+\sin^{2}\theta)}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{2(1+\cos\theta)}=\frac{1}{2}\sqrt{4\cos^{2}\frac{\theta}{2}}=\cos\frac{\theta}{2}$
$\displaystyle|\alpha_{2}|$ $\displaystyle=$
$\displaystyle\left|\frac{1}{2}\left(1-\cos\theta+i\sin\theta\right)\right|=\frac{1}{2}\sqrt{((1-\cos\theta)^{2}+\sin^{2}\theta)}$
(17) $\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{2(1-\cos\theta)}=\frac{1}{2}\sqrt{4\sin^{2}\frac{\theta}{2}}=\sin\frac{\theta}{2}.$

Next, we consider a controlled $R_{z}(\theta)$ unitary, as defined below.

$\displaystyle cR_{z}(\theta)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&e^{-i\theta/2}&0\\\
0&0&0&e^{i\theta/2}\end{bmatrix}=\alpha_{1}^{\prime}\mathbb{I}\mathbb{I}+\alpha_{2}^{\prime}\text{Z}\mathbb{I}+\alpha_{3}^{\prime}\mathbb{I}\text{Z}+\alpha_{4}^{\prime}\text{Z}\text{Z}$
$\displaystyle\mathbb{I}\otimes cR_{z}(\theta)$ $\displaystyle=$
$\displaystyle\alpha_{1}^{\prime}\mathbb{I}\mathbb{I}\mathbb{I}+\alpha_{2}^{\prime}\mathbb{I}\text{Z}\mathbb{I}+\alpha_{3}^{\prime}\mathbb{I}\mathbb{I}\text{Z}+\alpha_{4}^{\prime}\mathbb{I}\text{Z}\text{Z}$
(18)

Here $\alpha_{1}^{\prime}=\frac{1+\cos\frac{\theta}{2}}{2}$,
$\alpha_{2}^{\prime}=\frac{1-\cos\frac{\theta}{2}}{2}$,
$\alpha_{3}^{\prime}=\frac{i}{2}\sin\frac{\theta}{2}$,
$\alpha_{4}^{\prime}=-\alpha_{3}^{\prime}$. So,

$\displaystyle|\alpha_{1}^{\prime}|$ $\displaystyle=$
$\displaystyle\frac{1+\cos\frac{\theta}{2}}{2};\qquad|\alpha_{2}^{\prime}|=\frac{1-\cos\frac{\theta}{2}}{2}$
$\displaystyle|\alpha_{3}^{\prime}|=|\alpha_{4}^{\prime}|$ $\displaystyle=$
$\displaystyle\frac{\sin\frac{\theta}{2}}{2}.$ (19)

Next, we consider the following controlled rotation gate, that is widely used
in Quantum Fourier Transform [2].

$\displaystyle cR_{n}(\theta)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&e^{i\theta}\end{bmatrix}=\alpha_{1}^{\prime\prime}\mathbb{I}\mathbb{I}+\alpha_{2}^{\prime\prime}\mathbb{I}\text{Z}+\alpha_{3}^{\prime\prime}\text{Z}\mathbb{I}+\alpha_{4}^{\prime\prime}\text{Z}\text{Z}$
$\displaystyle\mathbb{I}\otimes cR_{n}(\theta)$ $\displaystyle=$
$\displaystyle\alpha_{1}^{\prime\prime}\mathbb{I}\mathbb{I}\mathbb{I}+\alpha_{2}^{\prime\prime}\mathbb{I}\mathbb{I}\text{Z}+\alpha_{3}^{\prime\prime}\mathbb{I}\text{Z}\mathbb{I}+\alpha_{4}^{\prime\prime}\mathbb{I}\text{Z}\text{Z}$
(20)

Here $\alpha_{1}^{\prime\prime}=\frac{3+e^{i\theta}}{4}$,
$\alpha_{2}^{\prime\prime}=\alpha_{3}^{\prime\prime}=\frac{1-e^{i\theta}}{4}$,
$\alpha_{4}^{\prime\prime}=-\alpha_{2}^{\prime\prime}$. So,

$\displaystyle|\alpha_{1}^{\prime\prime}|$ $\displaystyle=$
$\displaystyle\frac{1}{4}\sqrt{(3+\cos\theta)^{2}+\sin^{2}\theta}=\frac{1}{4}\sqrt{10+6\cos\theta}=\frac{1}{2}\sqrt{1+3\cos^{2}\frac{\theta}{2}}$
$\displaystyle|\alpha_{2}^{\prime\prime}|=|\alpha_{3}^{\prime\prime}|=|\alpha_{4}^{\prime\prime}|$
$\displaystyle=$
$\displaystyle\frac{1}{4}\sqrt{(1-\cos\theta)^{2}+\sin^{2}\theta}=\frac{1}{4}\sqrt{2-2\cos\theta}=\frac{1}{2}\sin\frac{\theta}{2}.$
(21)

Next, we consider the 2-qubit Given’s rotation unitary, that is used in
quantum chemistry [60, 61], quantum simulations [62, 63], quantum machine
learning [64], variational quantum algorithms [65].

$\displaystyle Givens(\theta)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}1&0&0&0\\\ 0&\cos(\theta)&-\sin(\theta)&0\\\
0&\sin(\theta)&\cos(\theta)&0\\\
0&0&0&1\end{bmatrix}=\alpha_{1}^{\prime\prime\prime}\mathbb{I}\mathbb{I}+\alpha_{2}^{\prime\prime\prime}\text{X}\text{Y}+\alpha_{3}^{\prime\prime\prime}\text{Y}\text{X}+\alpha_{4}^{\prime\prime\prime}\text{Z}\text{Z}$
$\displaystyle\mathbb{I}\otimes Givens(\theta)$ $\displaystyle=$
$\displaystyle\alpha_{1}^{\prime\prime\prime}\mathbb{I}\mathbb{I}\mathbb{I}+\alpha_{2}^{\prime\prime\prime}\mathbb{I}\text{X}\text{Y}+\alpha_{3}^{\prime\prime\prime}\mathbb{I}\text{Y}\text{X}+\alpha_{4}^{\prime\prime\prime}\mathbb{I}\text{Z}\text{Z}$
(22)

Here $\alpha_{1}^{\prime\prime\prime}=\frac{1+\cos\theta}{2}$,
$\alpha_{2}^{\prime\prime\prime}=i\frac{\sin\theta}{2}$,
$\alpha_{3}^{\prime\prime\prime}=-i\frac{\sin\theta}{2}$
$\alpha_{4}^{\prime\prime\prime}=\frac{1-\cos\theta}{2}$. So,

$\displaystyle|\alpha_{1}^{\prime\prime\prime}|$ $\displaystyle=$
$\displaystyle\frac{1+\cos\theta}{2};\qquad|\alpha_{4}^{\prime\prime\prime}|=\frac{1-\cos\theta}{2}$
$\displaystyle|\alpha_{2}^{\prime\prime\prime}|=|\alpha_{3}^{\prime\prime\prime}|$
$\displaystyle=$ $\displaystyle\frac{\sin\theta}{2}.$ (23)

We use $\text{diag}(d_{1},d_{2},\ldots,d_{N})$ to denote a $N\times N$
diagonal matrix with $d_{1},d_{2},\ldots,d_{N}$ as diagonal entries, and 0 in
the remaining places. Now we consider the following rotation gates, controlled
on 2 qubits.

$\displaystyle ccR_{n}(\theta)$ $\displaystyle=$
$\displaystyle\text{diag}\left(1,1,1,1,1,1,1,e^{i\theta}\right)$ (24)
$\displaystyle=$
$\displaystyle\widehat{\alpha_{1}}\mathbb{I}\mathbb{I}\mathbb{I}+\widehat{\alpha_{2}}\text{Z}\mathbb{I}\mathbb{I}+\widehat{\alpha_{3}}\mathbb{I}\text{Z}\mathbb{I}+\widehat{\alpha_{4}}\mathbb{I}\mathbb{I}\text{Z}+\widehat{\alpha_{5}}\text{Z}\text{Z}\mathbb{I}+\widehat{\alpha_{6}}\mathbb{I}\text{Z}\text{Z}+\widehat{\alpha_{7}}\text{Z}\mathbb{I}\text{Z}+\widehat{\alpha_{8}}\text{Z}\text{Z}\text{Z}$

Here $\widehat{\alpha_{1}}=\left(\frac{7+e^{i\theta}}{8}\right)$,
$\widehat{\alpha_{2}}=\widehat{\alpha_{3}}=\widehat{\alpha_{4}}=\widehat{\alpha_{8}}=\left(\frac{1-e^{i\theta}}{8}\right)$,
$\widehat{\alpha_{5}}=\widehat{\alpha_{6}}=\widehat{\alpha_{7}}=-\widehat{\alpha_{8}}$
and so,

$\displaystyle|\widehat{\alpha_{j}}|$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sqrt{(7+\cos\theta)^{2}+\sin^{2}\theta}=\frac{1}{4}\sqrt{9+7\cos^{2}\frac{\theta}{2}}\qquad[j=1]$
(25) $\displaystyle=$
$\displaystyle\frac{1}{8}\sqrt{(1-\cos\theta)^{2}+\sin^{2}\theta}=\frac{1}{4}\sin\frac{\theta}{2}\qquad[j\neq
1].$

Finally, we consider the following rotation gate.

$\displaystyle ccR_{z}(\theta)$ $\displaystyle=$
$\displaystyle\text{diag}\left(1,1,1,1,1,1,e^{-i\theta/2},e^{i\theta/2}\right)$
(26) $\displaystyle=$
$\displaystyle\widetilde{\alpha_{1}}\mathbb{I}\mathbb{I}\mathbb{I}+\widetilde{\alpha_{2}}\text{Z}\mathbb{I}\mathbb{I}+\widetilde{\alpha_{3}}\mathbb{I}\text{Z}\mathbb{I}+\widetilde{\alpha_{4}}\mathbb{I}\mathbb{I}\text{Z}+\widetilde{\alpha_{5}}\text{Z}\text{Z}\mathbb{I}+\widetilde{\alpha_{6}}\mathbb{I}\text{Z}\text{Z}+\widetilde{\alpha_{7}}\text{Z}\mathbb{I}\text{Z}+\widetilde{\alpha_{8}}\text{Z}\text{Z}\text{Z}$

Here $\widetilde{\alpha_{1}}=\frac{3+\cos\frac{\theta}{2}}{4}$,
$\widetilde{\alpha_{2}}=\widetilde{\alpha_{3}}=\frac{1-\cos\frac{\theta}{2}}{4}$,
$\widetilde{\alpha_{5}}=-\widetilde{\alpha_{2}}$,
$\widetilde{\alpha_{6}}=\widetilde{\alpha_{7}}=i\frac{\sin\frac{\theta}{2}}{4}$
and $\widetilde{\alpha_{4}}=\widetilde{\alpha_{8}}=-\widetilde{\alpha_{6}}$,
and so,

$\displaystyle|\widetilde{\alpha_{j}}|$ $\displaystyle=$
$\displaystyle\frac{3+\cos\frac{\theta}{2}}{4}\qquad[j=1]$ $\displaystyle=$
$\displaystyle\frac{1-\cos\frac{\theta}{2}}{4}\qquad[j=2,3,5]$
$\displaystyle=$
$\displaystyle\frac{\sin\frac{\theta}{2}}{4}\qquad[j=4,6,7,8].$

Since each of these are approximately implementable unitaries, so we can
obtain a lower bound on their $\epsilon$-Toffoli-count by plugging in the
respective values in Theorem 3.5.

###  3.3 Channel representation

An $n$-qubit unitary $U$ can be completely determined by considering its
action on a Pauli $P_{s}\in\mathcal{P}_{n}:P_{s}\rightarrow
UP_{s}U^{\dagger}$. The set of all such operators (with
$P_{s}\in\mathcal{P}_{n}$) completely determines $U$ up to a global phase.
Since $\mathcal{P}_{n}$ is a basis for the space of all Hermitian $2^{n}\times
2^{n}$ matrices we can write

$\displaystyle UP_{s}U^{\dagger}$ $\displaystyle=$
$\displaystyle\sum_{P_{r}\in\mathcal{P}_{n}}\widehat{U}_{rs}P_{r},\qquad\text{where
}\quad\widehat{U}_{rs}=\frac{1}{2^{n}}\text{Tr}\left(P_{r}UP_{s}U^{\dagger}\right).$
(27)

This defines a $4^{n}\times 4^{n}$ matrix $\widehat{U}$ with rows and columns
indexed by Paulis $P_{r},P_{s}\in\mathcal{P}_{n}$. We refer to $\widehat{U}$
as the channel representation of $U$. By Hermitian conjugation each entry of
$\widehat{U}$ is real. The channel representation respects matrix
multiplication and tensor product i.e. $\widehat{UW}=\widehat{U}\widehat{W}$
and $\left(\widehat{U\otimes W}\right)=\widehat{U}\otimes\widehat{W}$. Setting
$V=U^{\dagger}$ it follows that
$\widehat{U^{\dagger}}=\left(\widehat{U}\right)^{\dagger}$, and we see that
the channel representation $\widehat{U}$ is unitary.

######  Fact 3.7.

The channel representation inherits the decomposition from Theorem 3.2 (and in
this decomposition there is no global phase factor).

$\displaystyle\widehat{U}=\left(\prod_{j=m}^{1}\widehat{G_{P_{1j},P_{2j},P_{3j}}}\right)\widehat{C_{0}}\qquad[G_{P_{1j},P_{2j},P_{3j}}\in\mathcal{G}_{Tof}]$

A decomposition in which $m=\mathcal{T}^{of}(U)$ is called a Toffoli-count-
optimal decomposition. The channel representation of the unitaries in the
generating set $\mathcal{G}_{Tof}$ has some important properties, which we
summarize in the following theorem. These facilitate the design of proper data
structure and efficient algorithms for matrix operations like multiplication,
inverse, involving only integer arithmetic. These properties also give the
intuition to develop clever heuristics for Toffoli-count-optimal synthesis
algorithms. Further, we can draw inferences about the exact implementability
of certain unitaries.

######  Theorem 3.8.

Let $G_{P_{1},P_{2},P_{3}}\in\mathcal{G}_{Tof}$, where
$P_{1},P_{2},P_{3}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}\\}$. Then its
channel representation $\widehat{G_{P_{1},P_{2},P_{3}}}$ has the following
properties.

  1. 1.

The diagonal entries are $1$ or $\frac{1}{2}$.

  2. 2.

If a diagonal entry is 1 then all other entries in the corresponding row and
column is 0.

  3. 3.

If a diagonal entry is $\frac{1}{2}$ then there exists three entries in the
corresponding row and column that are equal to $\pm\frac{1}{2}$, rest is 0.

  4. 4.

Exactly $2^{2n-3}$ i.e. $\frac{1}{8}^{th}$ of the diagonal elements are $1$,
while the remaining are $\frac{1}{2}$.

The proof has been given in Appendix D. Since Cliffords map Paulis to Paulis,
up to a possible phase factor of -1, so we have the following.

######  Fact 3.9 ([46]).

Let $\widehat{\mathcal{C}_{n}}=\\{\widehat{C}:C\in\mathcal{C}_{n}\\}$. A
unitary $Q\in\widehat{\mathcal{C}_{n}}$ if and only if it has one non-zero
entry in each row and column, equal to $\pm 1$.

The above fact, along with Theorem 3.8 and Fact 3.7, implies the following.

######  Lemma 3.10.

If $U$ is a unitary exactly implementable by the Clifford+Toffoli gate set,
then each entry of $\widehat{U}$ belongs to the ring
$\mathbb{Z}\left[\frac{1}{2}\right]=\\{\frac{a}{2^{k}}:a\in
2\mathbb{Z}+1,k\in\operatorname{\mathbb{N}}\\}$.

Since the channel representation of the non-Clifford T gate do not belong to
this ring [47], so it cannot be exactly implemented by Clifford+Toffoli. On
the other hand, we know that Toffoli can be implemented exactly with at most 7
T gates [47] or with 4 T gates [28] if extra ancilla and classical resources
like measurement is used. Thus every unitary exactly implementable by the
Clifford+Toffoli gate set can be exactly implemented by the Clifford+T gate
set. But there exists unitaries exactly implementable by the Clifford+T gate
set that cannot be exactly implemented by the Clifford+Toffoli gate set. So
the following theorem follows.

######  Theorem 3.11.

Let $\mathcal{J}_{n}^{T}$ and $\mathcal{J}_{n}^{Tof}$ be the set of unitaries
exactly implementable by the Clifford+T and Clifford+Toffoli gate set,
respectively. Then $|\mathcal{J}_{n}^{Tof}|<|\mathcal{J}_{n}^{T}|.$

In a separate paper we prove that the channel representation of unitaries
exactly implementable by the V-basis gate set, has entries in the ring
$\mathbb{Z}\left[\frac{1}{\sqrt{5}}\right]$. Thus these gates are not exactly
implementable by the Clifford+Toffoli gate set. In fact, Toffoli also cannot
be exactly implemented by the V-basis gate set. Using Theorem 3.5 we can say
that we require $\Omega\left(\log_{4}\left(\frac{1}{\epsilon}\right)\right)$
Toffoli gates in order to implement an $\epsilon$-approximation of T and
V-basis gates. In Section 3.5 we show that the CS gate unitary, which is
defined as,

$\displaystyle\text{CS}=\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
0&0&0&i\end{bmatrix},$ (28)

can be implemented optimally with 3 Toffoli gates. In a separate paper we show
that Toffoli can be implemented with 3 CS gates. Thus we have the following
theorem.

######  Theorem 3.12.

Let $\mathcal{J}_{n}^{CS}$ and $\mathcal{J}_{n}^{Tof}$ be the set of unitaries
exactly implementable by the Clifford+CS and Clifford+Toffoli gate set,
respectively. Then $|\mathcal{J}_{n}^{Tof}|=|\mathcal{J}_{n}^{CS}|.$

Here we define the following.

######  Definition 3.13.

For any non-zero $v\in\mathbb{Z}\left[\frac{1}{2}\right]$ the smallest
2-denominator exponent, denoted by $\text{sde}_{2}$, is the smallest
$k\in\operatorname{\mathbb{N}}$ for which $v=\frac{a}{2^{k}}$, where $a\in
2\mathbb{Z}+1$.

We define $\text{sde}_{2}(0)=0$. For a $d_{1}\times d_{2}$ matrix $M$ with
entries over this ring we define

$\displaystyle\text{sde}_{2}(M)$ $\displaystyle=$
$\displaystyle\max_{a\in[d_{1}],b\in[d_{2}]}\text{sde}_{2}(M_{ab}).$

######  Fact 3.14.

Let $v_{1},v_{2},v_{3},v_{4}\in\mathbb{Z}\left[\frac{1}{2}\right]$ such that
$\max\\{\text{sde}_{2}(v_{1}),\text{sde}_{2}(v_{2}),\text{sde}_{2}(v_{3}),\text{sde}_{2}(v_{4})\\}=k$.
Then $\text{sde}_{2}\left(\frac{1}{2}(v_{1}\pm v_{2}\pm v_{3}\pm
v_{4})\right)=k+1$ or $\leq k$.

###### Proof.

Let $v_{1}=\frac{a}{2^{w}}$, $v_{2}=\frac{b}{2^{x}}$, $v_{3}=\frac{c}{2^{y}}$
and $v_{4}=\frac{d}{2^{z}}$, where $a,b,c,d\in 2\mathbb{Z}+1$ and
$w,x,y,z\in\operatorname{\mathbb{N}}$. Without loss of generality let
$\max\\{w,x,y,z\\}=x$

$\displaystyle\frac{1}{2}(v_{1}\pm v_{2}\pm v_{3}\pm v_{4})$ $\displaystyle=$
$\displaystyle\frac{1}{2^{x+1}}\left(a\cdot 2^{x-w}\pm b\pm c\cdot 2^{x-y}\pm
d\cdot 2^{x-z}\right)$

If the numerator is odd then there is no further reduction of the fraction and
$\text{sde}_{2}$ of the resultant sum is $k+1$. If the numerator is even then
there is a reduction of the fraction and the resultant $\text{sde}_{2}$ is at
most $k$. ∎

We represent each element
$v=\frac{a}{2^{k}}\in\mathbb{Z}\left[\frac{1}{2}\right]$ by a pair $[a,k]$ of
integers. In this way, all arithmetic operations involving elements within
this ring, for example, matrix multiplication can be done with integer
arithmetic. We have given the pseudocode for some of these operations in
Appendix E (Algorithms 3 and 4). This also implies that we do not require any
floating point operation for our Toffoli-count-optimal synthesis algorithms
for exactly implementable unitaries.

##### Compact representation of $\widehat{G_{P_{1},P_{2},P_{3}}}$ :

It is sufficient to represent the $2^{2n}\times 2^{2n}$ matrix
$\widehat{G_{P_{1},P_{2},P_{3}}}$ with an array of length $7\cdot 2^{2n-3}$.
For brevity, we write $\widehat{G_{P_{1},P_{2},P_{3}}}[P_{r},P_{s}]$ as
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,s]$. Each entry of this array is of the
form $(r,\pm s,\pm s^{\prime},\pm s^{\prime\prime})$, which implies
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,r]=\frac{1}{2}$,
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,s]=\pm\frac{1}{2}$,
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,s^{\prime}]=\pm\frac{1}{2}$ and
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,s^{\prime\prime}]=\pm\frac{1}{2}$. The
remaining entries of the matrix can be easily deduced, using Theorem 3.8. In
Algorithm 2 (Appendix E) we have provided a pseudocode for computing the
channel representation of these basis elements and storing them using this
compact representation.

##### Channel representation of $G_{P_{1},P_{2},P_{3}}^{\dagger}$ :

Since Toffoli is a self-inverse unitary, so

$\widehat{G_{P_{1},P_{2},P_{3}}}^{\dagger}=\widehat{G_{P_{1},P_{2},P_{3}}}.$

###  3.4 Multiplication by $\widehat{G_{P_{1},P_{2},P_{3}}}$

Let $U_{p}=\widehat{G_{P_{1},P_{2},P_{3}}}U$, where $U$ is another matrix of
dimension $2^{2n}\times 2^{2n}$. Then,

$\displaystyle
U_{p}[r,j]=\sum_{k=1}^{2^{2n}}\widehat{G_{P_{1},P_{2},P_{3}}}[r,k]U[k,j]$

and it can be computed very efficiently with the following observations.

  1. 1.

Suppose $\widehat{G_{P_{1},P_{2},P_{3}}}[r,r]=1$. From Theorem 3.8,
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,s]=0$ for each $s\neq r$. Then,

$\displaystyle
U_{p}[r,j]=\widehat{G_{P_{1},P_{2},P_{3}}}[r,r]U[r,j]=U[r,j]\qquad[\forall
j\in\\{1,\ldots,2^{2n}\\}]$

and so $U_{p}[r,.]\leftarrow U[r,.]$ i.e. the $r^{th}$ row of $U$ gets copied
into the $r^{th}$ row of $U_{p}$.

  2. 2.

Let $\widehat{G_{P_{1},P_{2},P_{3}}}[r,r]=\frac{1}{2}$. From Theorem 3.8, we
know there are 3 other non-zero off-diagonal elements. Let
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,s]=\pm\frac{1}{2}$,
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,s^{\prime}]=\pm\frac{1}{2}$ and
$\widehat{G_{P_{1},P_{2},P_{3}}}[r,s^{\prime\prime}]=\pm\frac{1}{2}$.

$\displaystyle U_{p}[r,j]$ $\displaystyle=$
$\displaystyle\widehat{G_{P_{1},P_{2},P_{3}}}[r,r]U[r,j]+\widehat{G_{P_{1},P_{2},P_{3}}}[r,s]U[s,j]+\widehat{G_{P_{1},P_{2},P_{3}}}[r,s^{\prime}]U[s^{\prime},j]$
$\displaystyle+\widehat{G_{P_{1},P_{2},P_{3}}}[r,s^{\prime\prime}]U[s^{\prime\prime},j]=\frac{1}{2}\left(U[r,j]\pm
U[s,j]\pm U[s^{\prime},j]\pm U[s^{\prime\prime},j]\right)$

Thus we see that the $r^{th}$ row of $U_{p}$ is a linear combination of the
$r^{th},s^{th},s^{\prime th},s^{\prime\prime th}$ rows of $U$, multiplied by
$\frac{1}{2}$. i.e. $U_{p}[r,.]\leftarrow\frac{1}{2}\left(U[r,.]\pm U[s,.]\pm
U[s^{\prime},.]\pm U[s^{\prime\prime},.]\right)$.

We have given a pseudocode of this procedure in Algorithm 5 (Appendix E).
Since $\frac{1}{8}^{th}$ of the diagonals in $\widehat{G_{P_{1},P_{2},P_{3}}}$
is 1 (Theorem 3.8), so we get the following result.

######  Theorem 3.15.

$U_{p}=\widehat{G_{P_{1},P_{2},P_{3}}}U$ can be computed in time
$O\left(7\cdot 2^{4n-3}\right)$.

Currently the fastest algorithm to multiply two $2^{2n}\times 2^{2n}$ matrices
has a time complexity of $2^{4.7457278n}$ [66]. This modest asymptotic
improvement in the complexity has a significant impact on the actual running
time, especially when many such matrix multiplications are performed, as in
our algorithms for exactly implementable unitaries. Also, using the
representations discussed earlier, we only work with integers. This is not
only faster, but also overcomes other floating point arithmetic issues like
precision, etc.

#####  $\text{sde}_{2}$ of product matrix :

Let $U_{p}=\widehat{G_{P_{1},P_{2},P_{3}}}\widehat{U}$, where
$\widehat{U}=\prod_{j}\widehat{G_{P_{1j},P_{2j},P_{3j}}}\widehat{C}$ is an
exactly implementable unitary and
$P_{1},P_{2},P_{3},P_{1j},P_{2j},P_{3j}\in\mathcal{P}_{n}$,
$C\in\mathcal{C}_{n}$. From Lemma 3.10,
$\widehat{U}[k,\ell]\in\mathbb{Z}\left[\frac{1}{2}\right]$, where
$k,\ell\in\\{1,\ldots,2^{2n}\\}$. Thus entries of $U_{p}$ are of the form
$\frac{1}{2}\left(v_{1}\pm v_{2}\pm v_{3}\pm v_{4}\right)$, where
$v_{1},v_{2},v_{3},v_{4}\in\mathbb{Z}\left[\frac{1}{2}\right]$. We can apply
Fact 3.14 and conclude that the resulting $\text{sde}_{2}$ can increase by 1,
remain unchanged or decrease.

Now, we know that when we multiply by inverse then the change in sde should be
reversed. Thus sde can decrease by at most 1 because it can increase by at
most 1.

######  Lemma 3.16.

Let $U_{p}=\widehat{G_{P_{1},P_{2},P_{3}}}\widehat{U}$, where
$\widehat{U}=\prod_{j}\widehat{G_{P_{1j},P_{2j},P_{3j}}}\widehat{C}$ and
$P_{1},P_{2},P_{3},P_{1j},P_{2j},P_{3j}\in\mathcal{P}_{n}$,
$C\in\mathcal{C}_{n}$. Then
$\text{sde}_{2}(U_{p})-\text{sde}_{2}(\widehat{U})\in\\{\pm 1,0\\}$.

###  3.5 Implementation results

We implemented our heuristic algorithm EXACT-TOF-OPT in Python on a machine
with Intel(R) Core(TM) i7-5500K CPU at 2.4GHz, having 8 GB RAM and running
Windows 10. This algorithm has been explained in Section 4.2.2. There are two
other algorithms described in Section 4 but their complexity is exponential in
number of qubits as well as Toffoli-count. So they become intractable on a
personal laptop, running on 1 core. In this heuristic algorithm we basically
try to guess a sequence
$\widetilde{U}=\prod_{j=1}^{m}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}$ such that
$U\widetilde{U}^{\dagger}=C_{0}\in\mathcal{C}_{n}$, where $U$ is the input
unitary. We work with the channel representation of unitaries. We prune the
search space depending upon the sde and Hamming weight of the intermediate
product unitaries. We have probed two rules in order to select the unitaries -
Rule A, where we group according to sde and Rule B, where we group according
to sde and Hamming weight. More details can be found in Section 4.2.2.

To test the heuristic algorithm we generated a number of random 3-qubit
unitaries with Toffoli-count at most 5 and 10. The procedure for generating
the random unitaries has been explained in Appendix E (RANDOM-CHAN-REP,
Algorithm 12). We found that using Rule A, the output Toffoli-count is at most
the input Toffoli-count for all the cases, but the time consumed is much
higher. In fact, the space utilized also becomes too high. Using Rule B, the
output is at most the input for nearly 80$\%$ of the cases, but the time
consumed is much less. Here, we mention that we terminate the algorithm if it
runs for more than 1 hour. In Table 3 we have given the implementation
results. We could not apply Rule A on input unitaries with 10 Toffolis because
it soon ran out of space. The average is computed over 10 input unitaries for
each group. These include only those whose implementation was not aborted due
to time being more than 1 hour or space consumption being too high. We have
also mentioned the maximum number of unitaries that need to be stored while
executing the program. This gives a sense about the space complexity.

We also implemented the unitary $\mathbb{I}\otimes\text{CS}$ using both the
rules and found that its Toffoli-count is 3. It took roughly 19 secs to get a
decomposition and maximum number of unitaries stored is at most 3. For this
particular case we did an exhaustive search in order to find a lower Toffoli-
count decomposition. But we found that this is the optimal Toffoli-count for
CS gate. A circuit implementing $\mathbb{I}\otimes\text{CS}$ has been shown in
Figure 1. More specifically, using our heuristic algorithm EXACT-TOF-OPT,
using both the rules we found that we obtain,

$\displaystyle\widehat{\mathbb{I}\otimes\text{CS}}=\left(\prod_{j=1}^{3}\widehat{G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}}\right)\widehat{C_{0}},\qquad[C_{0}\in\mathcal{C}_{3}]$

where
$(P_{1_{1}},P_{2_{1}},P_{3_{1}})=(\mathbb{I}\mathbb{I}\text{Z},\mathbb{I}\text{Z}\mathbb{I},\text{Z}\mathbb{I}\mathbb{I})$,
$(P_{1_{2}},P_{2_{2}},P_{3_{2}})=(\mathbb{I}\mathbb{I}\text{Y},\mathbb{I}\text{Z}\mathbb{I},\text{Z}\mathbb{I}\mathbb{I})$
and
$(P_{1_{3}},P_{2_{3}},P_{3_{3}})=(\mathbb{I}\mathbb{I}\text{X},\mathbb{I}\text{Z}\mathbb{I},\text{Z}\mathbb{I}\mathbb{I})$.
It is straightforward to verify the above equation. Now, as discussed in
Section 3.1, for each generating set element
$G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}$ we find conjugating Cliffords
$C_{j}\in\mathcal{C}_{3}$ such that
$C_{j}(\mathbb{I}\mathbb{I}\text{X})C_{j}^{\dagger}=P_{1_{j}}$,
$C_{j}(\mathbb{I}\text{Z}\mathbb{I})C_{j}^{\dagger}=P_{2_{j}}$ and
$C_{j}(\text{Z}\mathbb{I}\mathbb{I})C_{j}^{\dagger}=P_{3_{j}}$. Since
$\text{H}\text{X}\text{H}=\text{Z}$ and
$\text{S}\text{X}\text{S}^{\dagger}=\text{Y}$, so we obtain the circuit as
shown in Figure 1. We have boxed the corresponding parts of the circuit that
implement each generating set unitary.

Input count | Avg. time (A) | Avg. time (B) | Max $\\#$unitaries (A) | Max $\\#$unitaries (B)
---|---|---|---|---
5 |  $\approx$ 46 mins |  $\approx$ 100 s | 902 | 7
10 | - |  $\approx$ 250s | - | 10
Table 3: Implementation results of EXACT-TOF-COUNT on random 3-qubit
unitaries. In the first column we have the number of Toffolis in the input
unitary. The second and third column gives the average time using Divide-and-
Select Rule A and B respectively. The last two columns give the maximum number
of unitaries that need to be stored using rules A and B respectively. Figure
1: An implementation of $\mathbb{I}\otimes\text{CS}$ with 3 Toffoli gates.
$C_{0}$ is a 3-qubit Clifford. The boxed parts implement the three generating
set unitaries. The first box implements
$G_{\mathbb{I}\mathbb{I}\text{X},\mathbb{I}\text{Z}\mathbb{I},\text{Z}\mathbb{I}\mathbb{I}}$,
the second one implements
$G_{\mathbb{I}\mathbb{I}\text{Y},\mathbb{I}\text{Z}\mathbb{I},\text{Z}\mathbb{I}\mathbb{I}}$
and the third box implements
$G_{\mathbb{I}\mathbb{I}\text{Z},\mathbb{I}\text{Z}\mathbb{I},\text{Z}\mathbb{I}\mathbb{I}}$.

###  3.6 Discussions

In this paper we develop some theoretical concepts to better understand the
unitaries implementable by the Clifford+Toffoli gate set. First is the
generating set with which we can represent any unitary exactly or
approximately implmentable by this gate set. We derive a bound on the
cardinality of this generating set. Second is the analysis of the channel
representation of these generating set elements. With the help of these we are
able to derive some interesting results like a bound on the Toffoli-count of
arbitrary unitaries, impossibility of the exact implementation of T gate,
though we know that Toffoli can be exactly implemented by the Clifford+T gate
set. This clearly shows that the cardinality of the set of unitaries exactly
implementable by these gate sets is different. Such conclusions will help in
resource estimates in various applications.

With these theoretical concepts we develop Toffoli-count-optimal synthesis
algorithms for arbitrary multi-qubit unitaries. Such algorithms were
previously developed in the context of T-count-optimal synthesis. Complexity
of these algorithms depend on the cardinality of the generating set.
$|\mathcal{G}_{Tof}|$ is exponentially more than the size of the generating
set for Clifford+T. Also, when we work with channel representation the size of
the unitaries increases quadratically. So these algorithms soon become
intractable on a personal laptop. Thus we could not implement any of the
provable algorithms, because complexity of these depend exponentially on the
Toffoli-count. The heuristic algorithm with polynomial dependence on the
Toffoli-count, is much faster but it is not always accurate, especially with
Rule B. The accuracy depends on the characteristics of the channel
representation of the generating set. Thus inspite of using similar heuristics
the performance varied for Clifford+T and Clifford+Toffoli gate set. As
mentioned, due to the high complexity of the provable algorithms we could not
verify if the count returned by the heuristics were optimal. Since output
Toffoli-count was at most the Toffoli-count of the input circuit, so we
conjecture that they return optimal decomposition.

Inspite of this, we did get some nice result by implementing the heuristic
algorithm, for example, the exact implementation of CS gate. Due to its small
Toffoli count we could verify that in this case the heuristics did return the
optimal count. Since in a separate paper we show that Toffoli can be exactly
implemented by Clifford+CS, so it implies that the set of exactly
implementable unitaries by these two universal gate sets is the same. Thus
these heuristics are worth probing and can throw some interesting light on the
mathematical properties of the unitaries involved. More parameters can be
taken into account while dividing and selecting, the selecting criteria maybe
changed,etc.

As mentioned the theoretical concepts and techniques developed in this paper
helped us achieve some nice results and develop algorithms in the context of
Clifford+Toffoli, that could not be done before. These techniques have a
braoder applicability and can be used for other universal gate sets. More
efficient algorithms can be designed, possibly by sacrificing the optimality
criterion. These concepts can also be used to design Toffoli-depth-optimal
quantum circuits, analogous to the work in [48] for T-depth-optimality. We
leave these for future research.

##  4 Method

In this section we describe algorithms for Toffoli-count-optimal synthesis of
multi-qubit unitaries, both exactly and approximately implementable by
Clifford+Toffoli gate set.

###  4.1 Toffoli-count-optimal synthesis algorithm for approximately
implementable unitaries

Now we describe an algorithm for determining the $\epsilon$-Toffoli-count of
an $n$-qubit unitary $W$. Basically we use the procedure in [15], but with the
generating set $\mathcal{G}_{Tof}$. This algorithm has space and time
complexity $O\left(|\mathcal{G}_{Tof}|\right)$ and
$O\left(|\mathcal{G}_{Tof}|^{\mathcal{T}_{\epsilon}^{of}}\right)$
respectively. We briefly describe our algorithm here. More detail pseudocode
can be found in Appendix E. The optimization version, APPROX-TOF-OPT
(Algorithm 6), that finds the optimal Toffoli-count is an exhaustive search
procedure. In each iteration it calls the decision version APPROX-TOF-DECIDE
(Algorithm 7) in order to decide if $\mathcal{T}_{\epsilon}^{of}(W)=m$, for
increasing values of a variable $m$.

The main idea to solve the decision version is as follows. Suppose we have to
test if $\mathcal{T}_{\epsilon}^{of}(W)=m$. From the definitions given in
Section 2, we know that if this is true then $\exists
U\in\mathcal{J}_{n}^{Tof}$ such that $\mathcal{T}^{of}(U)=m$ and
$d(U,W)\leq\epsilon$. Let
$U=\left(\prod_{j=m}^{1}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\right)C_{0}$ where
$C_{0}\in\mathcal{C}_{n}$ and
$G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\in\mathcal{G}_{Tof}$. Also, we can write
$W=UE$, for some unitary $E$. Let
$\widetilde{U}=\prod_{j=m}^{1}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}$.

$\displaystyle
d(W,U)=\sqrt{1-\frac{|\text{Tr}(W^{\dagger}U)|}{N}}\leq\epsilon\implies|\text{Tr}(W^{\dagger}U)|\geq
N(1-\epsilon^{2})$ (29)
$\displaystyle\left|\text{Tr}\left(E^{\dagger}\right)\right|=\left|\text{Tr}\left(E\right)\right|\geq
N(1-\epsilon^{2})$ (30)

The above implies that $d(E,\mathbb{I})\leq\epsilon$. We have,

$\displaystyle\left|\text{Tr}\left(W^{\dagger}\left(\prod_{j=m}^{1}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\right)\right)\right|$
$\displaystyle=$
$\displaystyle\left|\text{Tr}\left(E^{\dagger}U^{\dagger}\left(\prod_{j=m}^{1}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\right)\right)\right|$
$\displaystyle=$
$\displaystyle\left|\text{Tr}\left(E^{\dagger}e^{-i\phi}C_{0}^{\dagger}\left(\prod_{j=1}^{m}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}^{\dagger}\right)\left(\prod_{j=m}^{1}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\right)\right)\right|$
$\displaystyle=$
$\displaystyle\left|\text{Tr}\left(E^{\dagger}C_{0}^{\dagger}\right)\right|=\left|\text{Tr}\left(EC_{0}\right)\right|.$

To test if we have guessed the correct $\widetilde{U}$, we perform the
following tests.

I. Amplitude test :  We use the following theorem, which says that if $E$ is
close to identity then distribution of absolute-value-coefficients of $EC_{0}$
and $C_{0}$ in the Pauli basis expansion, is almost similar.

######  Theorem 4.1 ([15]).

Let $E\in\mathcal{U}_{n}$ be such that
$\left|\text{Tr}\left(E\right)\right|\geq N\left(1-\epsilon^{2}\right)$, for
some $\epsilon\geq 0$. $C_{0}=\sum_{P\in\mathcal{P}_{n}}r_{P}P$ is an
$n$-qubit Clifford. If $\left|\left\\{P:r_{P}\neq 0\right\\}\right|=M$ then

$\displaystyle\frac{1-\epsilon^{2}}{\sqrt{M}}-\sqrt{M}\sqrt{2\epsilon^{2}-\epsilon^{4}}$
$\displaystyle\leq$
$\displaystyle\left|\text{Tr}\left(EC_{0}P_{1}\right)/N\right|\leq\frac{1}{\sqrt{M}}+\sqrt{M}\sqrt{2\epsilon^{2}-\epsilon^{4}}\quad[\text{if
}r_{P_{1}}\neq 0]$ (31) $\displaystyle\text{and}\quad 0$ $\displaystyle\leq$
$\displaystyle\left|\text{Tr}\left(EC_{0}P_{1}\right)/N\right|\leq\sqrt{M}\sqrt{2\epsilon^{2}-\epsilon^{4}}\quad[\text{if
}r_{P_{1}}=0]$ (32)

In fact, we can have a more general theorem as follows.

######  Theorem 4.2 ([15]).

Let $E\in\mathcal{U}_{n}$ be such that
$\left|\text{Tr}\left(E\right)\right|\geq N\left(1-\epsilon^{2}\right)$, for
some $\epsilon\geq 0$. $Q=\sum_{P\in\mathcal{P}_{n}}q_{P}P$ is an $n$-qubit
unitary. Then for each $P_{1}\in\mathcal{P}_{n}$,

$\displaystyle(1-\epsilon^{2})|q_{P_{1}}|-\sum_{P\in\mathcal{P}_{n}\setminus\\{P_{1}\\}}|q_{P}|\sqrt{2\epsilon^{2}-\epsilon^{4}}$
$\displaystyle\leq$
$\displaystyle\left|\text{Tr}\left(EQP_{1}\right)/N\right|\leq|q_{P_{1}}|+\sum_{P\in\mathcal{P}_{n}\setminus\\{P_{1}\\}}|q_{P}|\sqrt{2\epsilon^{2}-\epsilon^{4}}.$

So, we calculate $W^{\prime}=W^{\dagger}\widetilde{U}$ and then calculate the
set
$\mathcal{S}_{c}=\left\\{|\text{Tr}(W^{\prime}P)/N|:P\in\mathcal{P}_{n}\right\\}$,
of coefficients. Then we check if we can distinguish the following two subsets
$\mathcal{S}_{0}$ and $\mathcal{S}_{1}$, for some $1\leq M\leq N^{2}$. Further
details have been given in Algorithm 7.

$\displaystyle\mathcal{S}_{1}$ $\displaystyle=$
$\displaystyle\left\\{\left|\text{Tr}\left(EC_{0}P_{1}\right)/N\right|:r_{P_{1}}\neq
0\right\\}$ (33) $\displaystyle\mathcal{S}_{0}$ $\displaystyle=$
$\displaystyle\left\\{\left|\text{Tr}\left(EC_{0}P_{1}\right)/N\right|:r_{P_{1}}=0\right\\}$
(34)

After passing this test we have a unitary of the form $E^{\dagger}Q$ where $Q$
is a unitary with _equal_ or _nearly equal_ amplitudes or coefficients
(absolute value) at some points and zero or nearly zero at other points.

II. Conjugation test :  To ensure $Q$ is a Clifford i.e.
$W^{\prime}=E^{\dagger}C_{0}$ for some $C_{0}\in\mathcal{C}_{n}$, we perform
the conjugation test (Algorithm 8) for further verification. We use the
following theorem and its corollary.

######  Theorem 4.3 ([15]).

Let $E,Q\in\mathcal{U}_{n}$ such that $d(E,\mathbb{I})\leq\epsilon$.
$P^{\prime}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}\\}$ such that
$QP^{\prime}Q^{\dagger}=\sum_{P}\alpha_{P}P$, where $\alpha_{P}\in\mathbb{C}$.
Then for each $P^{\prime\prime}\in\mathcal{P}_{n}$,

$\displaystyle\min\\{0,|\alpha_{P^{\prime\prime}}|(1-4\epsilon^{2}+2\epsilon^{4})-2\epsilon\sum_{P\neq
P^{\prime\prime}}|\alpha_{P}|\\}$ $\displaystyle\leq$
$\displaystyle\left|\text{Tr}\left((E^{\dagger}QP^{\prime}Q^{\dagger}E)P^{\prime\prime}\right)\right|/N$
$\displaystyle\leq$
$\displaystyle\max\\{|\alpha_{P^{\prime\prime}}|+2\epsilon\sum_{P\neq
P^{\prime\prime}}|\alpha_{P}|,1\\}$

######  Corollary 4.4 ([15]).

Let $C_{0}\in\mathcal{C}_{n}$ and $P^{\prime}\in\mathcal{P}_{n}$ such that
$C_{0}P^{\prime}C_{0}^{\dagger}=\tilde{P}\in\mathcal{P}_{n}$. If
$E\in\mathcal{U}_{n}$ such that $d(E,\mathbb{I})\leq\epsilon$, then

$\displaystyle(1-4\epsilon^{2}+2\epsilon^{4})$ $\displaystyle\leq$
$\displaystyle\left|\text{Tr}\left((E^{\dagger}C_{0}P^{\prime}C_{0}^{\dagger}E)P^{\prime\prime}\right)\right|/N\leq
1\qquad\text{if }P^{\prime\prime}=\tilde{P}$ $\displaystyle 0$
$\displaystyle\leq$
$\displaystyle\left|\text{Tr}\left((E^{\dagger}C_{0}P^{\prime}C_{0}^{\dagger}E)P^{\prime\prime}\right)\right|/N\leq
2\epsilon\qquad\text{else. }$

The above theorem and corollary basically says that $EC_{0}$ (or
$E^{\dagger}C_{0}$) _approximately_ inherits the conjugation property of
$C_{0}$. For each $P^{\prime}\in\mathcal{P}_{n}$, if we expand
$C_{0}P^{\prime}C_{0}^{\dagger}$ in the Pauli basis then the absolute value of
the coefficients has value 1 at one point, 0 in the rest. If we expand
$EC_{0}P^{\prime}C_{0}^{\dagger}E^{\dagger}$ in the Pauli basis then one of
the coefficients (absolute value) will be almost 1 and the rest will be almost
0. From Theorem 4.3 this pattern will not show for at least one Pauli
$P^{\prime\prime\prime}\in\mathcal{P}_{n}$ if we have $E^{\dagger}Q$, where
$Q\notin\mathcal{C}_{n}$. If we expand
$EQP^{\prime\prime\prime}Q^{\dagger}E^{\dagger}$ or
$E^{\dagger}QP^{\prime\prime\prime}Q^{\dagger}E$ in the Pauli basis then the
_spike in the amplitudes_ will be in at least two points.

####  4.1.1 Synthesizing Toffoli-count-optimal circuits

From APPROX-TOF-DECIDE (Algorithm 7) we can obtain a sequence
$\\{U_{m},\ldots,U_{1}\\}$ of unitaries such that
$U=\left(\prod_{j=m}^{1}U_{j}\right)C_{0}e^{i\phi}$, where
$U_{j}=G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\in\mathcal{G}_{Tof}$,
$m=\mathcal{T}^{of}(U)=\mathcal{T}_{\epsilon}^{of}(W)$ and
$C_{0}\in\mathcal{C}_{n}$. We can efficiently construct circuits for each
$U_{j}\in\mathcal{G}_{Tof}$, as discussed earlier. In order to determine the
trailing $C_{0}$, we observe that from Algorithm 7 we can also obtain the
following information : (1) set $\mathcal{S}_{1}$, as defined in Equation 33,
(2) $M=|\mathcal{S}_{1}|$. Thus we can calculate $r=\frac{1}{\sqrt{M}}$ and
from step 7 we can actually calculate the set
$\widetilde{\mathcal{S}_{1}}=\left\\{(t_{P},P):|t_{P}|=|\text{Tr}(W^{\prime}P)/N|\in\mathcal{S}_{1}\right\\}$,
where
$W^{\prime}=W^{\dagger}\widetilde{U}=e^{-i\phi}E^{\dagger}C_{0}^{\dagger}$. We
perform the following steps.

  1. 1.

Calculate
$a_{P}=\frac{t_{P}}{t_{\mathbb{I}}}=\frac{\text{Tr}(W^{\prime}P)}{\text{Tr}(W^{\prime})}$,
where $(t_{P},P)\in\widetilde{\mathcal{S}_{1}}$ (or equivalently
$|t_{P}|\in\mathcal{S}_{1}$). We must remember that
$(t_{\mathbb{I}},\mathbb{I})\in\widetilde{\mathcal{S}_{1}}$.

It has been argued in [15] that for small enough $\epsilon$ (say
$\ll\frac{1}{M}$),
$a_{P}\approx\frac{\overline{r_{P}}}{\overline{r_{\mathbb{I}}}}$. From Fact
A.3 we know that
$\frac{|\overline{r_{P}}|}{|\overline{r_{\mathbb{I}}}|}=\frac{|r_{P}|}{|r_{\mathbb{I}}|}=1$.
So we adjust the fractions such that their absolute value is $1$. For small
enough $\epsilon$ this adjustment is not much and so with a slight abuse we
use the same notation for the adjusted values.

  2. 2.

Select $c,d\in\operatorname{\mathbb{R}}$ such that $c^{2}+d^{2}=r^{2}$. Let
$\widetilde{r_{\mathbb{I}}}=c+di$. Then we claim that the Clifford
$\widetilde{C_{0}}=\widetilde{r_{\mathbb{I}}}\sum_{P:r_{P}\neq
0}\overline{a_{P}}P$ is sufficient for our purpose.

It is not hard to see that $\widetilde{C_{0}}=e^{i\phi^{\prime}}C_{0}$ for
some $\phi^{\prime}\in[0,2\pi)$. Thus if
$U^{\prime}=\left(\prod_{j=m}^{1}U_{j}\right)\widetilde{C_{0}}$, then
$\mathcal{T}^{of}(U^{\prime})=\mathcal{T}^{of}(U)$ and
$d(U^{\prime},W)\leq\epsilon$. After determining the Clifford, we can
efficiently construct a circuit for it, for example, by using results in [67].

####  4.1.2 Complexity

Now, we analyse the time and space complexity of the above algorithm.

#### Time complexity

We first analyse the time complexity of $\mathcal{A}_{CONJ}$. The outer and
inner loop at steps 2 and 6, respectively, can run at most
$2^{2n}*2^{2n}=2^{4n}$ times, where $N=2^{n}$. The complexity of the inner
loop is dominated by the multiplications at step 13, which takes $O(2^{2n})$
time. Thus overall complexity of this sub-routine is $O(2^{6n})$.

Now we analyse the time complexity of APPROX-TOF-DECIDE, when testing for a
particular Toffoli-count $m$. The algorithm loops over all possible products
of $m$ unitaries $G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\in\mathcal{G}_{Tof}$. Thus
number of iterations (steps 2-14) is at most $|\mathcal{G}_{Tof}|^{m}\in
O\left(n^{2m}64^{(n-2)m}\right)\in O\left(n^{2m}2^{(6n-12)m}\right)$, from
Theorem 3.3 and this is the most expensive part of this algorithm. The $m$
multiplications at step 2 and the multiplications, trace and sorting at step 4
takes $O(m2^{2n}+2^{4n})$ time. The inner loop 5-13 happens $O(2^{2n})$ times.
Each of the list elements is checked and so step 8 has complexity $O(2^{2n})$.
Now let within the inner loop the conjugation test is called $k_{1}$ times. So
the loop 5-13 incurs a complexity $O(k_{1}\cdot 2^{6n}+(2^{2n}-k_{1})2^{2n})$,
when $k_{1}>0$, else it is $O(2^{4n})$. Let $k$ is the number of outer loops
(steps 2-14) for which conjugation test is invoked in the inner loop 5-13 and
$k_{1}$ is the maximum number of times this test is called within any 5-13
inner loop. Then the overall complexity of APPROX-TOF-DECIDE is
$O((n^{2m}2^{(6n-12)m}-k)\cdot(m2^{2n}+2^{4n})+k\cdot(m2^{2n}+2^{4n}+k_{1}2^{6n}+(2^{2n}-k_{1})2^{2n}))$.

The conjugation test is invoked only if a unitary passes the amplitude test.
We assume that $m<2^{2n}$ and the occurrence of non-Clifford unitaries with
equal amplitude is not so frequent such that $kk_{1}<n^{2m}2^{(6n-12)m}-k$.
Thus APPROX-TOF-DECIDE has a complexity of $O(n^{2m}2^{(6n-12)m})$, for one
particular value of $m$. Hence, the overall algorithm APPROX-TOF-OPT has a
time complexity
$O\left(n^{2\mathcal{T}_{\epsilon}^{of}(W)}2^{(6n-12)\mathcal{T}_{\epsilon}^{of}(W)}\right)$.

### Space complexity

The input to our algorithm is a $2^{n}\times 2^{n}$ unitary, with space
complexity $O(2^{2n})$. We need to store $\mathcal{G}_{Tof}$, either as
triples of Paulis or as matrices. It takes space
$O\left(|\mathcal{G}_{Tof}|\right)\in O\left(n^{2}2^{6n-12}\right)$ (Theorem
3.3). All the operations in APPROX-TOF-DECIDE or $\mathcal{A}_{CONJ}$ either
involve matrix multiplications or storing at most $2^{2n}$ numbers. Hence the
overall space complexity is $O\left(n^{2}2^{6n-12}\right)$.

###  4.2 Toffoli-count-optimal synthesis algorithms for exactly implementable
unitaries

In this section we describe two algorithms for synthesizing Toffoli-count-
optimal circuits of exactly implementable unitaries. For the first one we can
prove a rigorous bound on the optimality as well as space and time complexity,
both of whom are exponential in the Toffoli-count. The exponentially better
second algorithm has both space and time complexity polynomial in the Toffoli-
count, but the claimed optimality as well as complexity depends on some
conjectures.

####  4.2.1 Nested MITM

In this section we describe the nested meet-in-the-middle (MITM) framework,
introduced in [47], that can be used to synthesize Toffoli-count-optimal
circuits. It is based on the following result that can be proved with similar
arguments as in Lemma 3.1 of [48].

######  Lemma 4.5.

Let $S_{i}\subset U(2^{n})$ be the set of all unitaries implementable in
Toffoli-count $i$ over the Clifford+Toffoli gate set. Given a unitary $U$,
there exists a Clifford+Toffoli circuit of Toffoli-count $(m_{1}+m_{2})$
implementing $U$ if and only if $S_{m_{1}}^{\dagger}U\bigcap
S_{m_{2}}\neq\emptyset$.

The pseudocode of the algorithm has been given in Appendix E (Algorithm 9). We
briefly describe this algorithm here. The input consists of the unitary $U$,
generating set $\mathcal{G}_{Tof}$, test Toffoli-count $m$ and $c\geq 2$ that
indicates the extent of nesting or recursion we want in our meet-in-the-middle
approach. If $U$ is of Toffoli-count at most $m$ then the output consists of a
decomposition of $U$ into smaller Toffoli-count unitaries, else the algorithm
indicates that $U$ has Toffoli-count more than $m$.

The algorithm consists of $\left\lceil\frac{m}{c}\right\rceil$ iterations and
in the $i^{th}$ such iteration we generate circuits of depth $i$ ($S_{i}$) by
extending the circuits of depth $i-1$ ($S_{i-1}$) by one more level. Then we
use these two sets to search for circuits of depth at most $ci$. In order to
increase the efficiency of the search procedure we define the following. Let
$\widehat{\mathcal{J}_{n}^{Tof}}=\\{\widehat{W}:W\in\mathcal{J}_{n}^{Tof}\\}$.

######  Definition 4.6 (Coset label [46]).

Let $\widehat{W}\in\widehat{\mathcal{J}_{n}^{Tof}}$. Its coset label
$\widehat{W}^{(co)}$ is the matrix obtained by the following procedure. (1)
Rewrite $\widehat{W}$ so that each nonzero entry has a common denominator,
equal to $\sqrt{2}^{\text{sde}(\widehat{W})}$. (2) For each column of
$\widehat{W}$, look at the first non-zero entry (from top to bottom) which we
write as $v=\frac{a}{2^{\text{sde}(\widehat{W})}}$. If $a<0$, or if $a=0$ and
$b<0$, multiply every element of the column by $-1$. Otherwise, if $a>0$, or
$a=0$ and $b>0$, do nothing and move on to the next column. (3) After
performing this step on all columns, permute the columns so that they are
ordered lexicographically from left to right.

The following result shows that this is a kind of equivalence function.

######  Theorem 4.7 (Proposition 2 in [46]).

Let $\widehat{W},\widehat{V}\in\widehat{\mathcal{J}_{n}^{Tof}}$. Then
$\widehat{W}^{(co)}=\widehat{V}^{(co)}$ if and only if
$\widehat{W}=\widehat{V}\widehat{C}$ for some $C\in\mathcal{C}_{n}$.

In our algorithm the search is performed iteratively where in the $k^{th}$
($1\leq k\leq c-1$) round we generate unitaries of Toffoli-count at most $ki$
by taking $k$ unitaries
$\widehat{W_{1}},\widehat{W_{2}},\ldots,\widehat{W_{k}}$ where
$\widehat{W_{i}}\in S_{i}$ or $\widehat{W_{i}}\in S_{i-1}$. Let
$\widehat{W}=\widehat{W_{1}}\widehat{W_{2}}\ldots\widehat{W_{k}}$ and its
Toffoli-count is $k^{\prime}\leq ki$. We search for a unitary
$\widehat{W^{\prime}}$ in $S_{i}$ or $S_{i-1}$ such that
$\left(\widehat{W}^{\dagger}\widehat{U}\right)^{(co)}=\widehat{W^{\prime}}$.
By Lemma 4.5 if we find such a unitary it would imply that Toffoli-count of
$U$ is $k^{\prime}+i$ or $k^{\prime}+i-1$ respectively. In the other direction
if the Toffoli-count of $U$ is either $k^{\prime}+i$ or $k^{\prime}+i-1$ then
there should exist such a unitary $\widehat{W^{\prime}}$ in $S_{i}$ or
$S_{i-1}$ respectively. Thus if the Toffoli-count of $U$ is at most $m$ then
the algorithm terminates in one such iteration and returns a decomposition of
$U$. This proves the correctness of this algorithm.

##### Time and space complexity

With the concept of coset label we can impose a strict lexicographic ordering
on unitaries such that a set $S_{i}$ can be sorted with respect to this
ordering in $O\left(|S_{i}|\log|S_{i}|\right)$ time and we can search for an
element in this set in $O\left(\log|S_{i}|\right)$ time. Now consider the
$k^{th}$ round of the $i^{th}$ iteration (steps 9-9 of Algorithm 9 in Appendix
E). We build unitaries $\widehat{W}$ of Toffoli-count at most $ki$ using
elements from $S_{i}$ or $S_{i-1}$. Number of such unitaries is at most
$|S_{i}|^{k}$. Given a $\widehat{W}$, time taken to search for $W^{\prime}$ in
$S_{i}$ or $S_{i-1}$ such that
$\left(\widehat{W}^{\dagger}\widehat{U}\right)^{(co)}=\widehat{W^{\prime}}$ is
$O\left(\log|S_{i}|\right)$. Since $|S_{j}|\leq|\mathcal{G}_{Tof}|^{j}$, so
the $k^{th}$ iteration of the for loop within the $i^{th}$ iteration of the
while loop, takes time
$O\left(|\mathcal{G}_{Tof}|^{(c-1)i}\log|\mathcal{G}_{Tof}|\right)$. Thus the
time taken by the algorithm is
$O\left(|\mathcal{G}_{Tof}|^{(c-1)\left\lceil\frac{m}{c}\right\rceil}\log|\mathcal{G}_{Tof}|\right)$.

In the algorithm we store unitaries of depth at most
$\left\lceil\frac{m}{c}\right\rceil$. So the space complexity of the algorithm
is $O\left(|\mathcal{G}_{Tof}|^{\left\lceil\frac{m}{c}\right\rceil}\right)$.
Since $|\mathcal{G}_{Tof}|\in O\left(n^{2}2^{6n-12}\right)$ (Theorem 3.3), so
we have an algorithm with space complexity
$O\left(n^{2\lceil\frac{m}{c}\rceil}2^{(6n-12)\left\lceil\frac{m}{c}\right\rceil}\right)$
and time complexity
$O\left(n^{2(c-1)\lceil\frac{m}{c}\rceil}2^{(6n-12)(c-1)\left\lceil\frac{m}{c}\right\rceil}\right)$.

####  4.2.2 A polynomial complexity heuristic algorithm

In this section we describe a heuristic algorithm that returns a Toffoli-
count-optimal circuit with time and space complexity that are polynomial in
the Toffoli-count. Our algorithm is motivated by the polynomial complexity
T-count-optimal-synthesis algorithm in [47]. The conjectured exponential
advantage comes from an efficient way of pruning the search space. The
intermediate unitaries are grouped according to some properties and then one
such group is selected according to some criteria.

In the previous provable nested meet-in-the-middle algorithm we search for a
set of $\widehat{G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}}$ such that
$\widehat{U}^{{\dagger}}\prod_{j=\mathcal{T}^{of}(U)}^{1}\widehat{G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}}$
is $\widehat{C_{0}}$ for some $C_{0}\in\mathcal{C}_{n}$. Alternatively, we can
also search for a set of $\widehat{G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}}^{-1}$
such that
$\left(\prod_{j=1}^{\mathcal{T}^{of}(U)}\widehat{G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}}^{-1}\right)\widehat{U}$
is $\widehat{C_{0}}$ for some $C_{0}\in\mathcal{C}_{n}$, which is the approach
taken by our heuristic algorithm. The optimization version EXACT-TOF-OPT
(Algorithm 10 in Appendix E) is solved by repeatedly calling EXACT-TOF-DECIDE
(Algorithm 11 in Appendix E), the decision version, which tests whether the
Toffoli-count of an input unitary is at most an integer $m$ .

It will be useful to compare the procedure EXACT-TOF-DECIDE with building a
tree, whose nodes store some unitary and the edges represent
$\widehat{G_{P_{1},P_{2},P_{3}}}^{-1}$. Thus the unitary in a child node is
obtained by multiplying the unitary of the parent node with the
$\widehat{G_{P_{1},P_{2},P_{3}}}^{-1}$ represented by the edge. Number of
children nodes of any parent node is at most $|\mathcal{G}_{Tof}|-1$. The root
stores $\widehat{U}$ and we assume it is at depth 0 (Figure 2). We stop
building this tree the moment we reach a node which stores $\widehat{C_{0}}$
for some $C_{0}\in\mathcal{C}_{n}$. This implies that the path from the root
to this leaf gives a decomposition of $\widehat{U}$. In these kinds of
searches the size of the tree is one of the important factors that determine
the complexity of the algorithm. To reduce the complexity we prune this tree
(Figure 2).

At each level we group the nodes according to some “properties” or
“parameters” of the unitaries stored in them. We hope that these parameters
will “distinguish” the “correct” nodes at each level or depth of the tree and
thus we would get a decomposition. Note there can be more than one
decomposition of $\widehat{U}$ with the same or different Toffoli-count. By
“correct” nodes we mean those that occur in a Toffoli-count-optimal
decomposition of $\widehat{U}$. If the parameters always select only the
correct nodes then we expect to get much fewer nodes at each level of the tree
and the number of levels we have to build is $\mathcal{T}^{of}(U)$. But the
parameters we selected did not always distinguish the correct nodes and there
were some false positives. In order for the algorithm to succeed we have to be
careful so that we do not lose all correct nodes at any level and to make it
efficient we have to ensure that the number of false positives are not too
large and they eventually get eliminated.

We selected two parameters - sde and Hamming weight of the unitaries. We know
from Lemma 3.16 that sde of a child node unitary can differ by at most 1 from
its parent node unitary. While building a unitary we start with
$\widehat{\mathbb{I}_{n}}$ and multiply by subsequent
$\widehat{G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}}$ till we reach $\widehat{U}$. We
have observed that in most of these multiplications the sde increases by 1 and
the Hamming weight also gradually increases until it (Hamming weight) reaches
the maximum. So while doing the inverse operations i.e. decomposing
$\widehat{U}$ we expect that in most of the steps sde will decrease by 1 and
as we get close to the identity, the Hamming weight will also likely decrease.
If we multiply by a “wrong” $\widehat{G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}}^{-1}$
we expect to see the same changes with much less probability, which is the
probability of the false positives. This helps us distinguish the “correct”
and “wrong” nodes.

Specifically, at each level we divide the set $S$ of nodes into some subsets
and select one of them. Below are two possible ways to divide the nodes that
we have found effective. Suppose in one instance of EXACT-TOF-DECIDE, $m$ is
the target depth (maximum depth of the tree to be built) and we have built the
tree till depth $i$.

Divide-and-Select Rules :

  1. A

We divide into two sets - $S_{0}$ (sde increase) and $S_{1}$ (sde non-
increase). We select the set with the minimum cardinality such that the sde of
the unitaries in this set is at most $m-i$.

  2. B.

We divide into 9 sets - $S_{00}$ (both sde and Hamming weight increase),
$S_{01}$ (sde increase but Hamming weight decrease), $S_{02}$ (sde increase
but Hamming weight same), $S_{10}$ (sde decrease, Hamming weight increase) and
$S_{11}$ (both sde and Hamming weight decrease), $S_{12}$ (sde decrease but
Hamming weight same), $S_{20}$ (sde same but Hamming weight increase),
$S_{21}$ (sde same but Hamming weight decrease), $S_{22}$ (both sde and
Hamming weight same). We select the set with the minimum cardinality such that
the sde of the unitaries in this set is at most $m-i$.

We follow any one of the above methods of divide-and-select throughout the
algorithm. We note that in each of the above methods, if the cardinality of
the unitaries in the selected set is more than $m-i$ then it implies we cannot
get sde $0$ nodes within the next few levels.

Figure 2: The tree built in the heuristic procedure. At each level we select a
set of nodes according to some changes in the properties of the child
unitaries with respect to their parents, for example change in sde and Hamming
weight. Unitaries in the next level are generated from the selected set (black
nodes). We stop building the tree as soon as we reach a unitary with sde 0.
The path length to this node (in this case 3) is the Toffoli-count of the
unitary $U$. It also gives us a Toffoli-count-optimal decomposition of $U$.

##### Space and time complexity

The analysis of space and time complexity of the algorithm EXACT-TOF-OPT is
based on the following assumptions.

######  Conjecture 4.8.

Let $\widetilde{U}$ is the set of intermediate unitaries stored at any time,
or equivalently it is the set of selected children nodes at any level of the
above-mentioned tree. The cardinality of the set $\widetilde{U}$ in each
iteration of EXACT-TOF-DECIDE is at most $\text{poly}(n^{2}2^{6n-12})$, when
either method A or B of divide-and-select is applied. We get at least one
Toffoli-count-optimal decomposition.

Consider the sub-routine EXACT-TOF-DECIDE. There are $|\mathcal{G}_{Tof}|-1$
multiplications by $\widehat{G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}}^{-1}$ in each
iteration for each unitary in $\widetilde{U}$. And by the above conjecture
$|\widetilde{U}|\in\text{poly}(n^{2}2^{6n-12})$. Thus both the time and space
complexity of EXACT-TOF-DECIDE are in $\text{poly}(n^{2}2^{6n-12},m)$. We call
EXACT-TOF-DECIDE at most $\mathcal{T}^{of}(U)$ times to solve EXACT-TOF-OPT.
Thus space and time complexity of EXACT-TOF-OPT are in
$\text{poly}(n^{2}2^{6n-12},\mathcal{T}^{of}(U))$.

## Acknowledgements

P.Mukhopadhyay thanks Nathan Wiebe for helpful discussions.

## Code and data availability

All codes can be found online at the public repository :
https://github.com/PriyankaMukhopadhyay/Toffoli-opt

## References

  * [1] Daniel Gottesman.  The Heisenberg representation of quantum computers, talk at.  In International Conference on Group Theoretic Methods in Physics. Citeseer, 1998.
  * [2] Michael A Nielsen and Isaac L Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, 2010.
  * [3] Sergey Bravyi and Alexei Kitaev.  Universal quantum computation with ideal Clifford gates and noisy ancillas.  Physical Review A, 71(2):022316, 2005.
  * [4] Daniel Gottesman and Isaac L Chuang.  Quantum teleportation is a universal computational primitive.  arXiv preprint quant-ph/9908010, 1999.
  * [5] Sergey Bravyi and David Gosset.  Improved classical simulation of quantum circuits dominated by Clifford gates.  Physical Review Letters, 116(25):250501, 2016.
  * [6] Sergey Bravyi, Graeme Smith, and John A Smolin.  Trading classical and quantum computational resources.  Physical Review X, 6(2):021043, 2016.
  * [7] Aleksei Yur’evich Kitaev.  Quantum computations: algorithms and error correction.  Uspekhi Matematicheskikh Nauk, 52(6):53–112, 1997.
  * [8] CM Dawson and MA Nielsen.  The Solovay-Kitaev algorithm.  Quantum Information and Computation, 6(1):81–95, 2006.
  * [9] Neil J Ross.  Optimal ancilla-free Clifford+V approximation of Z-rotations.  Quantum Information & Computation, 15(11-12):932–950, 2015.
  * [10] Neil J Ross and Peter Selinger.  Optimal ancilla-free Clifford+T approximation of Z-rotations.  Quantum Information & Computation, 16(11-12):901–953, 2016.
  * [11] Alex Bocharov, Yuri Gurevich, and Krysta M Svore.  Efficient decomposition of single-qubit gates into V basis circuits.  Physical Review A, 88(1):012313, 2013.
  * [12] Andreas Blass, Alex Bocharov, and Yuri Gurevich.  Optimal ancilla-free Pauli+V circuits for axial rotations.  Journal of Mathematical Physics, 56(12):122201, 2015.
  * [13] Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca.  Practical approximation of single-qubit unitaries by single-qubit quantum Clifford and T circuits.  IEEE Transactions on Computers, 65(1):161–172, 2015.
  * [14] Priyanka Mukhopadhyay.  Composability of global phase invariant distance and its application to approximation error management.  Journal of Physics Communications, 5(11):115017, 2021.
  * [15] Vlad Gheorghiu, Michele Mosca, and Priyanka Mukhopadhyay.  T-count and T-depth of any multi-qubit unitary.  npj Quantum Information, 8(1):141, 2022.
  * [16] Vadym Kliuchnikov, Alex Bocharov, and Krysta M Svore.  Asymptotically optimal topological quantum compiling.  Physical Review Letters, 112(14):140504, 2014.
  * [17] Emil Génetay Johansen and Tapio Simula.  Fibonacci anyons versus Majorana fermions: A Monte Carlo approach to the compilation of braid circuits in $SU(2)_{k}$ anyon models.  PRX Quantum, 2(1):010334, 2021.
  * [18] Dorit Aharonov.  A simple proof that Toffoli and Hadamard are quantum universal.  arXiv preprint quant-ph/0301040, 2003.
  * [19] Yaoyun Shi.  Both Toffoli and controlled-NOT need little help to do universal quantum computing.  Quantum Information & Computation, 3(1):84–92, 2003.
  * [20] Matthew Amy, Andrew N Glaudell, and Neil J Ross.  Number-theoretic characterizations of some restricted Clifford+T circuits.  Quantum, 4:252, 2020.
  * [21] Renaud Vilmart.  A ZX-calculus with triangles for Toffoli-Hadamard, Clifford+T, and beyond.  arXiv preprint arXiv:1804.03084, 2018.
  * [22] Renaud Vilmart.  Rewriting and Completeness of Sum-Over-Paths in Dyadic Fragments of Quantum Computing.  arXiv preprint arXiv:2307.14223, 2023.
  * [23] Maria Luisa Dalla Chiara, Antonio Ledda, Giuseppe Sergioli, and Roberto Giuntini.  The Toffoli-Hadamard gate system: an algebraic approach.  Journal of philosophical logic, 42:467–481, 2013.
  * [24] Daniel Gottesman.  Stabilizer codes and quantum error correction.  arXiv preprint quant-ph/9705052, 1997.
  * [25] David G Cory, MD Price, W Maas, Emanuel Knill, Raymond Laflamme, Wojciech H Zurek, Timothy F Havel, and Shyamal S Somaroo.  Experimental quantum error correction.  Physical Review Letters, 81(10):2152, 1998.
  * [26] Arkady Fedorov, Lars Steffen, Matthias Baur, Marcus P da Silva, and Andreas Wallraff.  Implementation of a Toffoli gate with superconducting circuits.  Nature, 481(7380):170–172, 2012.
  * [27] Matthew D Reed, Leonardo DiCarlo, Simon E Nigg, Luyan Sun, Luigi Frunzio, Steven M Girvin, and Robert J Schoelkopf.  Realization of three-qubit quantum error correction with superconducting circuits.  Nature, 482(7385):382–385, 2012.
  * [28] Cody Jones.  Low-overhead constructions for the fault-tolerant Toffoli gate.  Physical Review A, 87(2):022328, 2013.
  * [29] Cody Jones.  Composite Toffoli gate with two-round error detection.  Physical Review A, 87(5):052334, 2013.
  * [30] Adam Paetznick and Ben W Reichardt.  Universal fault-tolerant quantum computation with only transversal gates and error correction.  Physical Review Letters, 111(9):090505, 2013.
  * [31] Theodore J Yoder.  Universal fault-tolerant quantum computation with Bacon-Shor codes.  arXiv preprint arXiv:1705.01686, 2017.
  * [32] Matthew B Hastings and Jeongwan Haah.  Distillation with sublogarithmic overhead.  Physical Review Letters, 120(5):050504, 2018.
  * [33] Laxmidhar Biswal, Debjyoti Bhattacharjee, Anupam Chattopadhyay, and Hafizur Rahaman.  Techniques for fault-tolerant decomposition of a multicontrolled Toffoli gate.  Physical Review A, 100(6):062326, 2019.
  * [34] Michael Beverland, Earl Campbell, Mark Howard, and Vadym Kliuchnikov.  Lower bounds on the non-Clifford resources for quantum computations.  Quantum Science and Technology, 5(3):035009, 2020.
  * [35] Priyanka Mukhopadhyay, Torin F Stetina, and Nathan Wiebe.  Quantum simulation of the first-quantized Pauli-Fierz hamiltonian.  arXiv preprint arXiv:2306.11198, 2023.
  * [36] Priyanka Mukhopadhyay, Nathan Wiebe, and Hong Tao Zhang.  Synthesizing efficient circuits for Hamiltonian simulation.  npj Quantum Information, 9(1):31, 2023.
  * [37] D Michael Miller, Robert Wille, and Zahra Sasanian.  Elementary quantum gate realizations for multiple-control Toffoli gates.  In 2011 41st IEEE International Symposium on Multiple-Valued Logic, pages 288–293. IEEE, 2011.
  * [38] HR Bhagyalakshmi and MK Venkatesha.  Optimized multiplier using reversible multi-control input Toffoli gates.  International Journal of VLSI Design & Communication Systems, 3(6):27, 2012.
  * [39] Rafaella Vale, Thiago Melo D Azevedo, Ismael CS Araújo, Israel F Araujo, and Adenilton J da Silva.  Circuit decomposition of multi-controlled special unitary single-qubit gates.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 2023.
  * [40] Mojtaba Noorallahzadeh, Mohammad Mosleh, Neeraj Kumar Misra, and Amin Mehranzadeh.  A novel design of reversible quantum multiplier based on multiple-control Toffoli synthesis.  Quantum Information Processing, 22(4):167, 2023.
  * [41] Matthew Amy, Andrew N Glaudell, Sarah Meng Li, and Neil J Ross.  Improved synthesis of Toffoli-Hadamard circuits.  arXiv preprint arXiv:2305.11305, 2023.
  * [42] Chia-Chun Lin, Amlan Chakrabarti, and Niraj K Jha.  Optimized quantum gate library for various physical machine descriptions.  IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 21(11):2055–2068, 2013.
  * [43] Marie Salm, Johanna Barzen, Frank Leymann, Benjamin Weder, and Karoline Wild.  Automating the comparison of quantum compilers for quantum circuits.  In Symposium and Summer School on Service-Oriented Computing, pages 64–80. Springer, 2021.
  * [44] Aram W Harrow, Benjamin Recht, and Isaac L Chuang.  Efficient discrete approximations of quantum gates.  Journal of Mathematical Physics, 43(9):4445–4451, 2002.
  * [45] Andrew N Glaudell, Neil J Ross, and Jacob M Taylor.  Optimal two-qubit circuits for universal fault-tolerant quantum computation.  npj Quantum Information, 7(1):1–11, 2021.
  * [46] David Gosset, Vadym Kliuchnikov, Michele Mosca, and Vincent Russo.  An algorithm for the T-count.  Quantum Information & Computation, 14(15-16):1261–1276, 2014.
  * [47] Michele Mosca and Priyanka Mukhopadhyay.  A polynomial time and space heuristic algorithm for T-count.  Quantum Science and Technology, 2021.
  * [48] Vlad Gheorghiu, Michele Mosca, and Priyanka Mukhopadhyay.  A (quasi-) polynomial time heuristic algorithm for synthesizing T-depth optimal circuits.  npj Quantum Information, 8(1):110, 2022.
  * [49] Brett Giles and Peter Selinger.  Exact synthesis of multiqubit Clifford+T circuits.  Physical Review A, 87(3):032332, 2013.
  * [50] Alexei Yu Kitaev, Alexander Shen, Mikhail N Vyalyi, and Mikhail N Vyalyi.  Classical and quantum computation.  Number 47. American Mathematical Soc., 2002.
  * [51] Austin G Fowler.  Constructing arbitrary Steane code single logical qubit fault-tolerant gates.  Quantum Information & Computation, 11(9-10):867–873, 2011.
  * [52] Timothée Goubault de Brugière, Marc Baboulin, Benoît Valiron, and Cyril Allouche.  Quantum circuits synthesis using Householder transformations.  Computer Physics Communications, 248:107001, 2020.
  * [53] Emanuel Malvetti, Raban Iten, and Roger Colbeck.  Quantum circuits for sparse isometries.  Quantum, 5:412, 2021.
  * [54] Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca.  Fast and efficient exact synthesis of single-qubit unitaries generated by Clifford and T gates.  Quantum Information & Computation, 13(7-8):607–630, 2013.
  * [55] Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca.  Asymptotically optimal approximation of single qubit unitaries by Clifford and T circuits using a constant number of ancillary qubits.  Physical Review Letters, 110(19):190502, 2013.
  * [56] Matthew Amy, Dmitri Maslov, Michele Mosca, and Martin Roetteler.  A meet-in-the-middle algorithm for fast synthesis of depth-optimal quantum circuits.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 32(6):818–830, 2013.
  * [57] Matthew Amy, Dmitri Maslov, and Michele Mosca.  Polynomial-time T-depth optimization of Clifford+T circuits via matroid partitioning.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 33(10):1476–1489, 2014.
  * [58] Ewout Van Den Berg and Kristan Temme.  Circuit optimization of Hamiltonian simulation by simultaneous diagonalization of pauli clusters.  Quantum, 4:322, 2020.
  * [59] Salman Beigi and Peter W Shor.  C_3, semi-Clifford and genralized semi-Clifford operations.  Quantum Information and Computation, 10(1-2):41–59, 2010.
  * [60] Juha J Vartiainen, Mikko Möttönen, and Martti M Salomaa.  Efficient decomposition of quantum gates.  Physical Review Letters, 92(17):177902, 2004.
  * [61] Juan Miguel Arrazola, Olivia Di Matteo, Nicolás Quesada, Soran Jahangiri, Alain Delgado, and Nathan Killoran.  Universal quantum circuits for quantum chemistry.  Quantum, 6:742, 2022.
  * [62] Ian D Kivlichan, Jarrod McClean, Nathan Wiebe, Craig Gidney, Alán Aspuru-Guzik, Garnet Kin-Lic Chan, and Ryan Babbush.  Quantum simulation of electronic structure with linear depth and connectivity.  Physical Review Letters, 120(11):110501, 2018.
  * [63] Google AI Quantum, Collaborators*â€ , Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon, Joseph C Bardin, Rami Barends, Sergio Boixo, Michael Broughton, Bob B Buckley, et al.  Hartree-Fock on a superconducting qubit quantum computer.  Science, 369(6507):1084–1089, 2020.
  * [64] Iordanis Kerenidis and Anupam Prakash.  Quantum machine learning with subspace states.  arXiv preprint arXiv:2202.00054, 2022.
  * [65] Alain Delgado, Juan Miguel Arrazola, Soran Jahangiri, Zeyue Niu, Josh Izaac, Chase Roberts, and Nathan Killoran.  Variational quantum algorithm for molecular geometry optimization.  Physical Review A, 104(5):052402, 2021.
  * [66] François Le Gall.  Powers of tensors and fast matrix multiplication.  In Proceedings of the 39th International Symposium on Symbolic and Algebraic Computation, pages 296–303. ACM, 2014.
  * [67] Scott Aaronson and Daniel Gottesman.  Improved simulation of stabilizer circuits.  Physical Review A, 70(5):052328, 2004.

##  Appendix A Some additional preliminaries

###  A.1 Cliffords and Paulis

The _single qubit Pauli matrices_ are as follows.

$\displaystyle\text{X}=\begin{bmatrix}0&1\\\
1&0\end{bmatrix}\qquad\text{Y}=\begin{bmatrix}0&-i\\\
i&0\end{bmatrix}\qquad\text{Z}=\begin{bmatrix}1&0\\\ 0&-1\end{bmatrix}$

Parenthesized subscripts are used to indicate qubits on which an operator acts
For example, $\text{X}_{(1)}=\text{X}\otimes\mathbb{I}^{\otimes(n-1)}$ implies
that Pauli X matrix acts on the first qubit and the remaining qubits are
unchanged.

The _$n$ -qubit Pauli operators_ are : $\mathcal{P}_{n}=\\{Q_{1}\otimes
Q_{2}\otimes\ldots\otimes
Q_{n}:Q_{i}\in\\{\mathbb{I},\text{X},\text{Y},\text{Z}\\}\\}.$

The _single-qubit Clifford group_ $\mathcal{C}_{1}$ is generated by the
Hadamard and phase gates : $\mathcal{C}_{1}=\braket{\text{H},\text{S}}$ where

$\displaystyle\text{H}=\frac{1}{\sqrt{2}}\begin{bmatrix}1&1\\\
1&-1\end{bmatrix}\qquad\text{S}=\begin{bmatrix}1&0\\\ 0&i\end{bmatrix}.$

When $n>1$ the _$n$ -qubit Clifford group_ $\mathcal{C}_{n}$ is generated by
these two gates (acting on any of the $n$ qubits) along with the two-qubit
$\text{CNOT}=\ket{0}\bra{0}\otimes\mathbb{I}+\ket{1}\bra{1}\otimes\text{X}$
gate (acting on any pair of qubits).

The Clifford group is special because of its relationship to the set of
$n$-qubit Pauli operators. Cliffords map Paulis to Paulis, up to a possible
phase of $-1$, i.e. for any $P\in\mathcal{P}_{n}$ and any
$C\in\mathcal{C}_{n}$ we have $CPC^{\dagger}=(-1)^{b}P^{\prime}$ for some
$b\in\\{0,1\\}$ and $P^{\prime}\in\mathcal{P}_{n}$. In fact, given two Paulis
(neither equal to the identity), it is always possible to efficiently find a
Clifford which maps one to the other.

######  Fact A.1 ([46]).

For any $P,P^{\prime}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}\\}$ there exists
a Clifford $C\in\mathcal{C}_{n}$ such that $CPC^{\dagger}=P^{\prime}$. A
circuit for $C$ over the gate set $\\{\text{H},\text{S},\text{CNOT}\\}$ can be
computed efficiently (as a function of $n$).

######  Fact A.2 ([15]).

Let $Q=\sum_{P\in\mathcal{P}_{n}}q_{P}P$ be the expansion of a matrix $Q$ in
the Pauli basis. Then

$q_{P}=\text{Tr}\left(QP\right)/N\qquad[N=2^{n}].$

Further, if $Q$ is a unitary then

$\sum_{P\in\mathcal{P}_{n}}\left|q_{P}\right|^{2}=1$

We observe the following when expanding a Clifford in the Pauli basis.

######  Fact A.3 ([59]).

If $C\in\mathcal{C}_{n}$ then for each $P\in\mathcal{P}_{n}$ $\exists
r_{P}\in\mathbb{C}$, such that $C=\sum_{P\in\mathcal{P}_{n}}r_{P}P$. Further,
if $r_{P},r_{P^{\prime}}\neq 0$ for any pair of $P,P^{\prime}$, then
$|r_{P}|=|r_{P^{\prime}}|=r$, for some $r\in\operatorname{\mathbb{R}}$.

######  Fact A.4.

If $[A,B]=0$ and $CAC^{\dagger}=A_{1}$, $CBC^{\dagger}=B_{1}$, then
$[A_{1},B_{1}]=0$.

###### Proof.

We have the following.

$\displaystyle A_{1}B_{1}$ $\displaystyle=$
$\displaystyle(CAC^{\dagger})(CBC^{\dagger})=C(AB)C^{\dagger}=C(BA)C^{\dagger}$
$\displaystyle=$ $\displaystyle(CBC^{\dagger})(CAC^{\dagger})=B_{1}A_{1}$

∎

##  Appendix B Generating set $\mathcal{G}_{Tof}$

######  Lemma B.1.

If $P_{1},P_{2},P_{3}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}\\}$ pair-wise
commutes, then

  1. 1.

$G_{P_{1},P_{2},P_{3}}$ remains same for any permutation of the Paulis.

  2. 2.

$G_{P_{1},-P_{1}P_{2},P_{3}}=G_{P_{1},P_{2},P_{3}}$.

  3. 3.

$G_{P_{1},-P_{2},P_{3}}=G_{P_{1},P_{2},P_{3}}C$, for some
$C\in\mathcal{C}_{n}$.

  4. 4.

$G_{P_{1},P_{2},P_{1}P_{2}}=\mathbb{I}$.

  5. 5.

$G_{P_{1},P_{2},P_{1}P_{2}P_{3}}=G_{P_{1},P_{2},P_{3}}$.

###### Proof.

(1) easily follows from the definition of $G_{P_{1},P_{2},P_{3}}$ in Equation
4.

(2) can be proved as follows.

$\displaystyle G_{P_{1},-P_{1}P_{2},P_{3}}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(P_{1}-P_{1}P_{2}+P_{3}+P_{1}P_{1}P_{2}+P_{1}P_{2}P_{3}-P_{1}P_{3}-P_{1}P_{1}P_{2}P_{3}\right)$
$\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(P_{1}+P_{2}+P_{3}-P_{1}P_{2}-P_{2}P_{3}-P_{3}P_{1}+P_{1}P_{2}P_{3}\right)=G_{P_{1},P_{2},P_{3}}$

Now we prove (3). We observe that $\text{TOF}_{i,j;k}$ can be alternatively
written as follows.

$\displaystyle\text{TOF}_{i,j;k}$ $\displaystyle=$
$\displaystyle\exp\left(\frac{i\pi}{8}(\mathbb{I}-\text{X}_{(i)})(\mathbb{I}-\text{Z}_{(j)})(\mathbb{I}-\text{Z}_{(k)})\right)$
$\displaystyle\implies CU_{tof}C^{\dagger}$ $\displaystyle=$
$\displaystyle\exp\left(\frac{i\pi}{8}(\mathbb{I}-Q)\right)=G_{P_{1},P_{2},P_{3}}$

and so

$\displaystyle G_{P_{1},P_{2},P_{3}}$ $\displaystyle=$ $\displaystyle
C(\text{TOF}_{i,j;k})C^{\dagger}$ $\displaystyle=$
$\displaystyle\exp\left(\frac{i\pi}{8}(\mathbb{I}-P_{1}-P_{2}-P_{3}+P_{1}P_{2}+P_{2}P_{3}+P_{3}P_{1}-P_{1}P_{2}P_{3})\right).$

And,

$\displaystyle G_{P_{1},-P_{2},P_{3}}$ $\displaystyle=$
$\displaystyle\exp\left(\frac{i\pi}{8}(\mathbb{I}-P_{1}+P_{2}-P_{3}-P_{1}P_{2}-P_{2}P_{3}+P_{1}P_{3}+P_{1}P_{2}P_{3})\right)$
$\displaystyle=$ $\displaystyle
G_{P_{1},P_{2},P_{3}}\exp\left(\frac{i\pi}{8}(2P_{2}-2P_{1}P_{2}-2P_{2}P_{3}+2P_{1}P_{2}P_{3})\right)$
$\displaystyle=$ $\displaystyle
G_{P_{1},P_{2},P_{3}}\exp\left(\frac{i\pi}{4}P_{2}\right)\exp\left(-\frac{i\pi}{4}P_{1}P_{2}\right)\exp\left(-\frac{i\pi}{4}P_{2}P_{3}\right)\exp\left(\frac{i\pi}{4}P_{1}P_{2}P_{3}\right)$

where the last line follows since each term commutes with the others. We prove
that $\exp\left(\frac{i\pi}{4}P_{2}\right)$,
$\exp\left(-\frac{i\pi}{4}P_{1}P_{2}\right)$,
$\exp\left(-\frac{i\pi}{4}P_{2}P_{3}\right)$ and
$\exp\left(\frac{i\pi}{4}P_{1}P_{2}P_{3}\right)$ are Cliffords and this will
prove (3). Let $B\in\pm\mathcal{P}_{n}$. Then,

$\displaystyle\exp(\frac{i\pi}{4}P_{2})B\exp(\frac{i\pi}{4}P_{2})^{\dagger}$
$\displaystyle=$
$\displaystyle\frac{\mathbb{I}+iP_{2}}{\sqrt{2}}B\frac{\mathbb{I}-iP_{2}}{\sqrt{2}}$
$\displaystyle=$ $\displaystyle
P_{2}\left(\frac{P_{2}B+BP_{2}}{2}\right)+i\left(\frac{P_{2}B-BP_{2}}{2}\right)$
$\displaystyle=$ $\displaystyle B\quad\text{if }P_{2}B=BP_{2};\quad
iP_{2}B\quad\text{if }P_{2}B=-BP_{2}$

Since $B$ is a Pauli, so $iP_{2}B$ is also a Paulis This proves that
$\exp\left(\frac{i\pi}{4}P_{2}\right)\in\mathcal{C}_{n}$.

Similarly,
$\exp\left(-\frac{i\pi}{4}P_{1}P_{2}\right),\exp\left(-\frac{i\pi}{4}P_{2}P_{3}\right),\exp\left(\frac{i\pi}{4}P_{1}P_{2}P_{3}\right)\in\mathcal{C}_{n}$.
This proves (3).

(4) can be proved as follows.

$\displaystyle G_{P_{1},P_{2},P_{1}P_{2}}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(P_{1}+P_{2}+P_{1}P_{2}-P_{1}P_{2}-P_{2}P_{1}P_{2}-P_{1}P_{2}P_{1}+P_{1}P_{2}P_{1}P_{2}\right)$
$\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(P_{1}+P_{2}-P_{1}-P_{2}+\mathbb{I}\right)=\mathbb{I}$

(5) can be proved as follows.

$\displaystyle G_{P_{1},P_{2},P_{1}P_{2}P_{3}}$ $\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(P_{1}+P_{2}+P_{1}P_{2}P_{3}-P_{1}P_{2}-P_{2}P_{1}P_{2}P_{3}-P_{1}P_{2}P_{3}P_{1}+P_{1}P_{2}P_{1}P_{2}P_{3}\right)$
$\displaystyle=$
$\displaystyle\frac{3}{4}\mathbb{I}+\frac{1}{4}\left(P_{1}+P_{2}+P_{3}-P_{1}P_{2}-P_{2}P_{3}-P_{3}P_{1}+P_{1}P_{2}P_{3}\right)=G_{P_{1},P_{2},P_{3}}$

∎

We know that $\mathcal{G}_{Tof}$ is the set of unitaries

$\displaystyle
G_{P_{1},P_{2},P_{3}}=\frac{3}{4}\mathbb{I}_{n}+\frac{1}{4}(P_{1}+P_{2}+P_{3}-P_{1}P_{2}-P_{2}P_{3}-P_{3}P_{1}+P_{1}P_{2}P_{3}),$

such that if $\pi_{3}$ is a permutation of 3 elements, then

  1. 1.

$P_{1},P_{2},P_{3}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}\\}$, $P_{1}\neq
P_{2}\neq P_{3}$ and $[P_{1},P_{2}]=[P_{2},P_{3}]=[P_{3},P_{1}]=0$.

  2. 2.

$(P_{1},P_{2},P_{3})\equiv(P_{1},\pm
P_{1}P_{2},P_{3})\equiv(P_{1},P_{2},P_{1}P_{2}P_{3})\equiv\pi_{3}(P_{1},P_{2},P_{3})$.
The equivalence implies that we include only one of them in the set.

We call this a _generating set_ (modulo Clifford), because any unitary exactly
implementable by the Clifford+Toffoli gate set can be written (up to a global
phase) as a product of unitaries from this set and a Clifford.

Before proceeding to prove a bound on $|\mathcal{G}_{Tof}|$, let us consider a
set $\mathcal{G}$ of pairs of Paulis $(P_{1},P_{2})$ such that

  1. 1.

$P_{1},P_{2}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}\\}$, $P_{1}\neq P_{2}$
and $[P_{1},P_{2}]=0$.

  2. 2.

$(P_{1},P_{2})\equiv(P_{2},P_{1})\equiv(P_{1},\pm P_{1}P_{2})$. The
equivalence implies that we include only one of them in the set.

This set is relevant for constructing a generating set ($\mathcal{G}_{CS}$)
for the unitaries exactly implementable by the Clifford+CS gate set. In fact,
a bound on $|\mathcal{G}|$ provides a bound on $|\mathcal{G}_{CS}|$, which we
have proved in a separate paper. We are reproducing the proof, for
completeness.

######  Theorem B.2.

$\displaystyle|\mathcal{G}_{CS}|\leq\frac{1}{8}(16^{n}-13^{n}-4^{n}+1)+\frac{1}{12}(12^{n}-2\cdot
6^{n})$

###### Proof.

Let
$P=\bigotimes_{j=1}^{n}P_{j}\in\mathcal{P}_{n}\setminus\\{\mathbb{I}_{n}\\}$
such that it is the tensor product of $m$ non-identity single-qubit Paulis.
Without loss of generality, let us assume that $P_{j}\neq\mathbb{I}$ when
$1\leq j\leq m$. Let $S_{1Pm}$ is the set of non-identity Paulis that commute
with $P$ and have $\mathbb{I}$ on the first $m$ 1-qubit subspaces i.e.
$Q=\bigotimes_{j=1}^{n}Q_{j}\in S_{1Pm}$ if $Q_{j}=\mathbb{I}$ for $1\leq
j\leq m$. This implies $[P,Q]=0$. When $m+1\leq j\leq n$, then
$Q_{j}\in\mathcal{P}_{1}$ and so there can be 4 possible options. Excluding
the condition when $Q_{j}=\mathbb{I}$ for all $j$, we have
$|S_{1Pm}|=4^{n-m}-1$.

Let $S_{2Pm}$ is the set of non-identity Paulis that commute with $P$ and have
$\mathbb{I}$ on the last $n-m$ 1-qubit subspaces i.e.
$Q^{\prime}=\bigotimes_{j=1}^{n}Q_{j}^{\prime}\in S_{2Pm}$ if
$Q_{j}^{\prime}=\mathbb{I}$ for $m+1\leq j\leq n$. Here, we want to enforce
the constraint that $(P,P^{\prime})\equiv(P,\pm PP^{\prime})$ i.e. $\nexists
R\in S_{2Pm}$ such that $Q^{\prime}=\pm PR$. Now, $[P,Q^{\prime}]=0$ if and
only if $k=\left|\\{j:P_{j}\neq Q_{j}^{\prime};1\leq j\leq m\\}\right|$ is
even. In the remaining $m-k$ subspaces $Q_{j}^{\prime}=P_{j}$ or $\mathbb{I}$.
There can be at most $2^{m-k}$ combinations. In the $k$ unequal (but neither
is identity) subspaces there can be at most $2^{k}$ combinations, but out of
these half can be obtained by myltiplying $P$ with the other half. For
example, consider $P=\text{X}\text{X}$, then
$Q^{\prime}\in\\{\text{Y}\text{Y},\text{Y}\text{Z}\\}$. Though
$\\{\text{Z}\text{Z},\text{Z}\text{Y}\\}$ also commute with $P$ but these can
be obtained by multiplying $P$ with YY and YZ. Thus, there can be
$2^{k-1}\cdot 2^{m-k}=2^{m-1}$ possibilities for $Q^{\prime}$. Now we can
select $k$ subspaces in $\binom{m}{k}$ ways. When $m<n$ then $k$ can vary from
$0$ to $m$ or $m-1$, whichever is even. When $m=n$ then $k$ varies from $2$ to
$m$ or $m-1$, because we want to avoid the all identity case. Hence,

$\displaystyle|S_{2Pm}|$ $\displaystyle\leq$ $\displaystyle
2^{m-1}\sum_{k=0}^{m^{\prime}}\binom{m}{k}:=2^{m-1}\cdot
A\quad[m<n,m^{\prime}=m\text{ or }m-1]$ $\displaystyle S_{2Pm}|$
$\displaystyle\leq$ $\displaystyle
2^{m-1}\sum_{k=2}^{m^{\prime}}\binom{m}{k}=2^{m-1}\left(\sum_{\ell=0}^{k}\binom{m}{\ell}-1\right):=2^{m-1}\cdot(A-1)\quad[m=n,m^{\prime}=m\text{
or }m-1]$

When $m$ is even then,

$\displaystyle A$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{m}\binom{m}{k}=\frac{1}{2}\left(\sum_{i=0}^{m}\binom{m}{i}+\sum_{i=0}^{m}(-1)^{i}\binom{m}{i}\right)=\frac{1}{2}((1+1)^{m}+(1-1)^{m})=2^{m-1};$

while when $m-1$ is even then,

$\displaystyle A$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{m-1}\binom{m}{k}=\frac{1}{2}\left(\sum_{i=0}^{m}\binom{m}{i}+\sum_{i=0}^{m}(-1)^{i}\binom{m}{i}\right)=2^{m-1}.$

So,

$\displaystyle|S_{2Pm}|$ $\displaystyle\leq$ $\displaystyle
2^{m-1}2^{m-1}=2^{2m-2}\qquad[1\leq m<n]$ (35) $\displaystyle|S_{2Pm}|$
$\displaystyle\leq$ $\displaystyle
2^{m-1}\left(2^{m-1}-1\right)=2^{2m-2}-2^{m-1}\qquad[m=n]$

Let $d_{Pm}$ is the number of Paulis that commute with $P$ and we maintain the
equivalence $(P,Q)\equiv(P,\pm PQ)$ i.e. we count either one of $P$ and $PQ$.
Then,

$\displaystyle d_{P}$ $\displaystyle\leq$ $\displaystyle|S_{1Pm}|\leq
4^{n-1}-1\qquad[m=1]$ $\displaystyle\leq$ $\displaystyle|S_{2Pm}|\leq
2^{2n-2}-2^{n-1}\qquad[m=n]$ $\displaystyle\leq$
$\displaystyle|S_{1Pm}|\cdot|S_{2Pm}|\leq\left(2^{2n-2m}-1\right)2^{2m-2}=2^{2n-2}-2^{2m-2}=\frac{1}{4}\left(4^{n}-4^{m}\right)\qquad[1<m<n]$

We count the number of Paulis $P^{\prime}=\bigotimes_{j=1}^{n}P_{j}^{\prime}$,
which are tensor product of $m$ non-identity Paulis. We denote it by $d_{m}$.
There can be $\binom{n}{m}3^{m}$ of them. When $m=n$, among them there are few
which are products of the other. Specifically we can take 2 possible values
for $P_{1}^{\prime}$ and 3 possible values for the remaining $P_{j}^{\prime}$,
where $1<j\leq n$. The remaining $n$-qubit Paulis that have $m$ non-identity
tensored components, can be obtained as product of these. And we do not count
these, due to the equivalence constraint. When $m=1$ there can be 3
possibilities for the only non-identity component. So number of possible
Paulis with $m$ non-identity tensored components is at most
$3\cdot\binom{n}{1}$ if $m=1$; $\binom{n}{n}2\cdot 3^{n-1}$ if $m=n$ and
$\binom{n}{m}3^{m}$ if $1<m<n$. Hence,

$\displaystyle d_{m}$ $\displaystyle\leq$
$\displaystyle\binom{n}{1}3(4^{n-1}-1)=3n(4^{n-1}-1)\qquad\text{when }m=1;$
$\displaystyle d_{m}$ $\displaystyle\leq$
$\displaystyle\binom{n}{m}\frac{3^{m}}{4}\left(2^{2n}-2^{2m}\right)=\binom{n}{m}\left(4^{n-1}\cdot
3^{m}-\frac{1}{4}\cdot 12^{m}\right)\qquad\text{when }n>m\geq 2;$
$\displaystyle d_{m}$ $\displaystyle\leq$ $\displaystyle\binom{n}{n}2\cdot
3^{n-1}\left(2^{2n-2}-2^{n-1}\right)=2\left(12^{n-1}-6^{n-1}\right)\qquad\text{when
}m=n.$

Total number of Pauli pairs is,

$\displaystyle\sum_{m}d_{m}$ $\displaystyle\leq$ $\displaystyle
3n(4^{n-1}-1)+2(12^{n-1}-6^{n-1})+\sum_{m=2}^{n-1}\binom{n}{m}\left(4^{n-1}\cdot
3^{m}-\frac{1}{4}\cdot 12^{m}\right)$ $\displaystyle:=$
$\displaystyle\frac{3n}{4}4^{n}-3n+\frac{12^{n}}{6}-\frac{6^{n}}{3}+B;$

where

$\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{4^{n}}{4}\left[\sum_{m=0}^{n}\binom{n}{m}3^{m}-\binom{n}{0}-\binom{n}{1}3-\binom{n}{n}3^{n}\right]-\frac{1}{4}\left[\sum_{m=0}^{n}\binom{n}{m}12^{m}-\binom{n}{0}-\binom{n}{1}12-\binom{n}{n}12^{n}\right]$
$\displaystyle=$
$\displaystyle\frac{4^{n}}{4}\left[4^{n}-1-3n-3^{n}\right]-\frac{1}{4}\left[13^{n}-1-12n-12^{n}\right]=\frac{16^{n}}{4}-\frac{13^{n}}{4}-4^{n}\frac{3n+1}{4}+3n+\frac{1}{4};$

and so,

$\displaystyle\sum_{m}d_{m}$ $\displaystyle\leq$
$\displaystyle\frac{16^{n}}{4}-\frac{13^{n}}{4}+\frac{12^{n}}{6}-\frac{6^{n}}{3}-\frac{4^{n}}{4}+\frac{1}{4}.$

Since the number of pairs are permutation invariant i.e.
$(P_{1},P_{2})\equiv(P_{2},P_{1})$, so,

$\displaystyle|\mathcal{G}_{CS}|=\frac{1}{2}\sum_{m}d_{m}\leq\frac{16^{n}}{8}-\frac{13^{n}}{8}+\frac{12^{n}}{12}-\frac{6^{n}}{6}-\frac{4^{n}}{8}+\frac{1}{8}.$
(36)

This proves the theorem.

∎

Now we prove a bound on $|\mathcal{G}_{Tof}|$.

######  Theorem B.3.

$\displaystyle|\mathcal{G}_{Tof}|$ $\displaystyle\leq$
$\displaystyle\frac{1}{384}[64^{n}-52^{n}-17(16^{n}-4^{n})]+\frac{1}{576}[48^{n}-2\cdot
24^{n}]-\frac{1}{36}[12^{n}-2\cdot 6^{n}]+\frac{1}{24}(13^{n}-1)$
$\displaystyle\in$ $\displaystyle O\left(n^{2}64^{n-2}\right)$

###### Proof.

First we count the number of commuting Pauli pairs $(P,Q)$ such that
$(P,Q)\equiv(P,\pm PQ)$. Let $\mathcal{G}$ denote the set of all such pairs.
From Theorem B.2 we know that

$\displaystyle|\mathcal{G}|\leq\frac{1}{8}\left(16^{n}-13^{n}-4^{n}+1\right)+\frac{1}{12}\left(12^{n}-2\cdot
6^{n}\right).$

Now, let us fix a particular $(P,Q)$ and count the number of commuting Paulis
$R$ that satisfy the constraints in Lemma B.1. Now number of Paulis that
commute with both $P$ and $Q$ but is not equal to $\mathbb{I},P,Q,PQ$ is
$\frac{4^{n}}{4}-4=4^{n-1}-4$. The constraints imply that
$(P,Q,R)\equiv(P,Q,PR)\equiv(P,Q,QR)\equiv(P,Q,PQR)$. Let $R^{\prime}$ is a
Pauli that commutes with both $P,Q$. If $R^{\prime}=PQ$ then
$PR^{\prime}=PPQ=Q$, $QR^{\prime}=PQR$, $PQR^{\prime}=PQPR=QR$. Thus we get a
cycle. Similary, for $R^{\prime}=QR,PQR$. For each commuting Pauli $R$ we get
such a cycle of length 4. Thus if $S_{PQ}$ is the number of commuting Paulis
$R$ that satisfy the constraints, then

$\displaystyle|S_{PQ}|\leq\frac{4^{n-1}-1}{4}=4^{n-2}-1.$

So total number of triples which we can have is

$\displaystyle|\mathcal{G}|\cdot|S_{PQ}|$ $\displaystyle\leq$
$\displaystyle(4^{n-2}-1)\left[\frac{1}{8}\left(16^{n}-13^{n}-4^{n}+1\right)+\frac{1}{12}\left(12^{n}-2\cdot
6^{n}\right)\right]$ $\displaystyle=$
$\displaystyle\frac{64^{n}}{128}-\frac{52^{n}}{128}+\frac{48^{n}}{192}-\frac{24^{n}}{96}-\frac{16^{n}}{128}+\frac{4^{n}}{128}-\frac{16^{n}}{8}+\frac{13^{n}}{8}-\frac{12^{n}}{12}+\frac{6^{n}}{6}+\frac{4^{n}}{8}-\frac{1}{8}$
$\displaystyle=$
$\displaystyle\frac{1}{128}[64^{n}-52^{n}-17(16^{n}-4^{n})]+\frac{1}{192}[48^{n}-2\cdot
24^{n}]-\frac{1}{12}[12^{n}-2\cdot 6^{n}]+\frac{1}{8}(13^{n}-1)$

Now there are $\binom{3}{2}=3$ ways in which we can first place the pair
$(P,Q)$ and then the position of $R$ is uniquely determined. Thus,

$\displaystyle|\mathcal{G}_{tof}|$ $\displaystyle\leq$
$\displaystyle\frac{1}{3}\left(|\mathcal{G}|\cdot|S_{PQ}|\right)$
$\displaystyle\leq$
$\displaystyle\frac{1}{384}[64^{n}-52^{n}-17(16^{n}-4^{n})]+\frac{1}{576}[48^{n}-2\cdot
24^{n}]-\frac{1}{36}[12^{n}-2\cdot 6^{n}]+\frac{1}{24}(13^{n}-1).$

We observe the fact that for any $a,b\in\operatorname{\mathbb{R}}$ such that
$a\geq b>1$ and $m\in\mathbb{Z}_{+}$ we have

$\displaystyle
a^{m}-b^{m}=(a-b)(a^{m-1}+a^{m-2}b+a^{m-3}b^{2}+\ldots+b^{m-1}),$

and so

$\displaystyle(a-b)mb^{m-1}\leq a^{m}-b^{m}\leq(a-b)ma^{m-1}.$

Thus,

$\displaystyle|\mathcal{G}_{Tof}|$ $\displaystyle\leq$
$\displaystyle\frac{1}{384}\left[(64-52)n64^{n-1}-17(16-4)n4^{n-1}\right]+\frac{2(48-24)}{576}n48^{n-1}+\frac{(13-1)}{24}n13^{n-1}$
$\displaystyle\leq$
$\displaystyle\frac{12n}{384}\left[64^{n-1}-4^{n-1}\right]+\frac{n48^{n-1}}{12}+\frac{n13^{n-1}}{2}$
$\displaystyle\leq$
$\displaystyle\frac{n}{32}(n-1)(64-4)64^{n-2}+\frac{n48^{n-2}48}{12}+\frac{n13^{n-2}13}{2}$
$\displaystyle=$
$\displaystyle\frac{15n(n-1)}{8}64^{n-2}+4n48^{n-2}+\frac{13n}{2}13^{n-2}\in
O\left(n^{2}64^{n-2}\right),$

thus proving the theorem. ∎

##  Appendix C Bound on Toffoli-count of arbitrary unitaries

In this section we derive some results required to prove a lower bound on the
Toffoli-count of arbitrary $n$-qubit unitaries, as discussed in Section 3.2.

######  Lemma C.1.

If $\widetilde{U}=\prod_{j=1}^{m}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}$, where
$G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}=a\mathbb{I}+bQ_{j}$,
$Q_{j}=P_{1_{j}}+P_{2_{j}}+P_{3_{j}}-P_{1_{j}}P_{2_{j}}-P_{2_{j}}P_{3_{j}}-P_{3_{j}}P_{1_{j}}+P_{1_{j}}P_{2_{j}}P_{3_{j}}$,
$a=\frac{3}{4}$ and $b=\frac{1}{4}$, then

$\displaystyle\widetilde{U}$ $\displaystyle=$ $\displaystyle
a^{m}\mathbb{I}+a^{m-1}b\left(\sum_{j}Q_{j}\right)+a^{m-2}b^{2}\left(\sum_{j_{1}<j_{2}}Q_{j_{1}}Q_{j_{2}}\right)+a^{m-3}b^{3}\left(\sum_{j_{1}<j_{2}<j_{3}}Q_{j_{1}}Q_{j_{2}}Q_{j_{3}}\right)+\cdots$
$\displaystyle\cdots+ab^{m-1}\left(\sum_{j_{1}<j_{2}<\cdots<j_{m-1}}Q_{j_{1}}Q_{j_{2}}\cdots
Q_{j_{m-1}}\right)+b^{m}\prod_{j=1}^{m}Q_{j}$

###### Proof.

We prove this by induction. The statement holds trivially for $m=1$. For $m=2$
we have

$\displaystyle\prod_{j=1}^{2}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}=(a\mathbb{I}+bQ_{1})(a\mathbb{I}+bQ_{2})=a^{2}\mathbb{I}+ab(Q_{1}+Q_{2})+b^{2}Q_{1}Q_{2}.$

We assume that the statement holds for $m-1$, that is

$\displaystyle\prod_{j=1}^{m-1}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}$
$\displaystyle=$ $\displaystyle
a^{m-1}\mathbb{I}+a^{m-2}b\left(\sum_{j=1}^{m-1}Q_{j}\right)+a^{m-3}b^{2}\left(\sum_{\begin{subarray}{c}j_{1},j_{2}=1\\\
j_{1}<j_{2}\end{subarray}}^{m-1}Q_{j_{1}}Q_{j_{2}}\right)+\cdots$
$\displaystyle\cdots+ab^{m-2}\left(\sum_{\begin{subarray}{c}j_{1},j_{2},\ldots,j_{m-2}=1\\\
j_{1}<j_{2}<\cdots<j_{m-2}\end{subarray}}^{m-1}Q_{j_{1}}Q_{j_{2}}\cdots
Q_{j_{m-2}}\right)+b^{m-1}\prod_{j=1}^{m-1}Q_{j}.$

Therefore,

$\displaystyle\widetilde{U}$ $\displaystyle=$
$\displaystyle\left(\prod_{j=1}^{m-1}G_{P_{1_{j}},P_{2_{j}},P_{3_{j}}}\right)(a\mathbb{I}+bQ_{m})$
$\displaystyle=$
$\displaystyle\left(a^{m-1}\mathbb{I}+a^{m-2}b\left(\sum_{j=1}^{m-1}Q_{j}\right)+a^{m-3}b^{2}\left(\sum_{\begin{subarray}{c}j_{1},j_{2}=1\\\
j_{1}<j_{2}\end{subarray}}^{m-1}Q_{j_{1}}Q_{j_{2}}\right)+\cdots\right.$