\big\lvert \bigwedge_{i=1}^{d+1} P'_{\alpha_i} \big\rvert &\equiv \sum_{k_{d+1}=0}^{1} \ldots \sum_{k_1=0}^{1} \big\lvert \left[ \bigwedge_{i=1}^{d+1} (k_i P_{\alpha_i} \vee (1 - k_i) z_{\alpha_i} \bs y) \right] \big\rvert \pmod{2} \\
&\equiv \big\lvert \bigwedge_{i=1}^{d+1} P_{\alpha_i} \big\rvert \pmod{2}.
\end{aligned}
\end{equation}

Based on Theorem <ref>, we can prove the following subadditivity theorem which is a generalization of Theorem <ref>.

Let $d$ be a non-negative integer, let $U_1 \in \mathcal{D}^C_{d+1}$, and let $U_2 \in \mathcal{D}_{d+1}$ where $\mathcal{D}_{d+1}$ is the diagonal subgroup of the $(d+1)$th level of the Clifford hierarchy, and $\mathcal{D}_{d+1}^C$ is the subgroup of $\mathcal{D}_{d+1}$ which can be implemented using only $C^{\otimes d}Z$ gates.
If $\tau[U_1], \tau[U_2] > 0$, then $\tau[U_1 \otimes U_2] < \tau[U_1] + \tau[U_2]$ where $\tau[U]$ is the optimal number of $R_Z(\pi/2^d)$ gates required to implement $U$ without ancillary qubits over the $\{\mathrm{CNOT}$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ gate set.

Let $W = \begin{bmatrix} P & Q \end{bmatrix}$ be a parity table such that $P$ and $Q$ are the parity tables associated with the implementation of $U_1$ and $U_2$ and which have $\tau[U_1]$ and $\tau[U_2]$ columns respectively.
Then, because $U_1 \in \mathcal{D}^C_{d+1}$, $P$ satisfies

\begin{equation}
\big\lvert \bigwedge_{i=1}^d P_{\alpha_i} \big\rvert \equiv 0 \pmod{2}
\end{equation}

for all $\alpha_1, \ldots, \alpha_d$ satisfying $0 \leq \alpha_1 < \ldots \leq \alpha_d < n$ where $n$ is the number of qubits on which $U_1$ is acting.
Let $\bs z = P_{:,i} \oplus Q_{:,j}$ for any $i$ and $j$ satisfying $0\leq i < \tau[U_1]$, $0\leq j < \tau[U_2]$, and let $P'$ be a parity table such that

P' =
\begin{cases}
P \oplus \bs z \bs 1^T &\text{if $\tau[U_1] \equiv 0 \pmod{2}$},\\
\begin{bmatrix} P \oplus \bs z \bs 1^T & \bs z \end{bmatrix} &\text{otherwise}.\\
\end{cases}

Then, as stated by Theorem <ref>, the parity table $P'$ satisfies

\begin{equation}
\big\lvert \bigwedge_{i=1}^{d+1} P'_{\alpha_i} \big\rvert \equiv \big\lvert \bigwedge_{i=1}^{d+1} P_{\alpha_i} \big\rvert \pmod{2}
\end{equation}

for all $0 \leq \alpha_1 \leq \ldots \leq \alpha_{d+1} < n$.
And so the parity table $W' = \begin{bmatrix} P' & Q \end{bmatrix}$, which has at most one more column than $W$, also satisfies

\begin{equation}\label{eq:d_w_equality}
\big\lvert \bigwedge_{i=1}^{d+1} W'_{\alpha_i} \big\rvert \equiv \big\lvert \bigwedge_{i=1}^{d+1} W_{\alpha_i} \big\rvert \pmod{2}
\end{equation}

for all $0 \leq \alpha_1 \leq \ldots \leq \alpha_{d+1} < n$.
However, we can notice that $P'_{:,i} = Q_{:,j}$.
Therefore, by removing these two columns from $W'$ Equation <ref> still holds and $W'$ has at least one less column than $W$.
The parity table $W'$ implements the unitary $U_1 \otimes U_2$ up to an operator implementable over the $\{\mathrm{CNOT},$ $R_Z(2\pi/2^d)\}$ gate set and has at most $\tau[U_1] + \tau[U_2] - 1$ columns, thus we have $\tau[U_1 \otimes U_2] \leq \tau[U_1] + \tau[U_2] - 1 < \tau[U_1] + \tau[U_2]$.

Based on this subadditivity theorem and on Theorem <ref>, we can derive the following upper on the number of $R_Z(\pi/2^d)$ gates in a $\{\mathrm{CNOT}$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ circuit.

Let $d$ be a non-negative integer.
The number of $R_Z(\pi/2^d)$ gates in an $n$-qubits $\{\mathrm{CNOT}$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ circuit can be upper bounded by

\begin{equation}
2 \left\lfloor \sum_{i=1}^{d} \frac{{n \choose i}}{2} \right\rfloor + 1 \leq \sum_{i=0}^{d} {n \choose i}
\end{equation}

in polynomial time.

Let $U$ be a unitary gate implementable over the $\{\mathrm{CNOT}$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ gate set, let $P$ be a parity table of size $n \times m$ which implements $U$ up to an operator implementable over the $\{\mathrm{CNOT}$, $R_Z(2\pi/2^d)\}$ gate set, and let $L$ be a matrix with rows labelled by $(\alpha_1\ldots\alpha_k)$ such that

\begin{align}
L_{\alpha_1\ldots\alpha_k} &= \bigwedge_{i=1}^k P_{\alpha_i}
\end{align}

for all $\alpha_1, \ldots, \alpha_k$ satisfying $0 \leq \alpha_i < \ldots < \alpha_k < n$ and $k$ satisfying $1 \leq k \leq d$.
If $P$ has strictly more than $\sum_{i=0}^{d} {n \choose i}$ columns then we can necessarily find a non-zero vector $\bs y$ satisfying $L \bs y = \bs 0$ and $\bs y \neq \bs 1$ because $L$ has $\sum_{i=1}^{d} {n \choose i}$ rows.
Note that such vector $\bs y$ necessarily satisfies Equations <ref> of Theorem <ref>.
We can then divide $P$ into two non-empty parity tables $P^{(1)}$ and $P^{(2)}$ where the column $P_{:,i}$ belongs to $P^{(1)}$ if and only if $y_i = 1$ and to $P^{(2)}$ otherwise.
The parity tables $P^{(1)}$ and $P^{(2)}$ are implementations of some unitary gates $U_1 \in \mathcal{D}^C_{d+1}$ and $U_2 \in \mathcal{D}_{d+1}$ respectively.
The subadditivity theorem (Theorem <ref>) can then be exploited to reduce the number of columns of $P$.
Let $\bs z = P_{:,i} \oplus P_{:,j}$ where $i$ and $j$ are satisfying $y_i = 1$ and $y_j = 0$, and let

P' =
\begin{cases}
P \oplus \bs z \bs y^T &\text{if $\lvert \bs y \rvert \equiv 0 \pmod{2}$},\\
\begin{bmatrix} P \oplus \bs z \bs y^T & \bs z \end{bmatrix} &\text{otherwise}.\\
\end{cases}

By Theorem <ref>, $P'$ implements the same unitary gate as $P$ up to an operator implementable over the $\{\mathrm{CNOT}$, $R_Z(2\pi/2^d)\}$ gate set.
The parity table $P'$ has at most one more column than $P$ and we have $P'_{:,i} = P'_{:,j}$.
Therefore, the columns $i$ and $j$ can be removed from $P'$, which entails that $P'$ has at least one less column than $P$.
We showed that if the number of columns of $P$ is strictly greater than $\sum_{i=0}^{d} {n \choose i}$, then the number of columns of $P$ can be reduced by at least one in polynomial time.
Moreover, if the number of columns of $P$ is equal to $\sum_{i=0}^{d} {n \choose i}$ and is even, then we can necessarily find a non-zero vector $\bs y$ satisfying $L \bs y = \bs 0$ because $L$ has $\sum_{i=1}^{d} {n \choose i}$ rows.
If there exist $i$ such that $y_i = 0$ then we can reduced the number of columns of $P$ as described above.
Otherwise we must have $\lvert \bs y \rvert \equiv 0 \pmod{2}$, and so $\bs y$ satisfies the Equations of Theorem <ref>.
Therefore, if $P' = P \oplus \bs z \bs y^T$ where $\bs z$ is equal to the $i$th column of $P$ for any $i$, then $P'$ implements the same unitary gate as $P$ up to an operator implementable over the $\{\mathrm{CNOT}$, $R_Z(2\pi/2^d)\}$ gate set.
The parity table $P'$ has the same number of columns as $P$ but its $i$th column is equal to the null vector and can therefore be removed, which leads to a parity table containing $\sum_{i=1}^{d} {n \choose i}$ columns.

The polynomial-time procedure described above to reduce the number of columns of $P$ can be repeated until $P$ has a number of columns lower or equal to

\begin{equation}
2 \left\lfloor \sum_{i=1}^{d} \frac{{n \choose i}}{2} \right\rfloor + 1 \leq \sum_{i=0}^{d} {n \choose i}.
\end{equation}

The overall complexity of the algorithm described in the proof of Theorem <ref> is $\mathcal{O}(n^{2d + 1} m)$ where $m$ is the number of $R_Z(\pi/2^d)$ gates in the initial circuit.
For $d=0$, the algorithm yields a parity table which contains at most one column, which is optimal.
For $d=1$, the algorithm corresponds to Lempel's matrix factorization algorithm <cit.>, which is optimal.
For $d=2$, the algorithm is similar to the  algorithm (Algorithm <ref>).
Note that the  algorithm has a complexity of $\mathcal{O}(n^2 m^3)$ instead of $\mathcal{O}(n^5 m)$ because of the heuristic used by the algorithm to better optimize the number of $T$ gates.
The same heuristic can be used in the case where $d>2$, which would lead to an algorithm having a complexity of $\mathcal{O}(n^d m^3)$.

§.§ Generalization of the  algorithm for optimizing $R_Z(\pi/2^d)$ gates

In this subsection we generalize the  algorithm (Algorithm <ref>) for the optimization of the number of $R_Z(\pi/2^d)$ gates in a $\{\mathrm{CNOT}$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ circuit.

Let $d$ be a non-negative integer, let $P$ be a parity table of size $n \times m$ and $P' = P \oplus \bs z \bs y^T$ where $\bs z$ and $\bs y$ are vectors of size $n$ and $m$ respectively and such that $\lvert \bs y \rvert \equiv 0 \pmod{2}$.
And let $L$ and $X$ be matrices and $\bs v$ be a vector, all with rows labelled by $(\alpha_1 \ldots \alpha_k)$ such that

\begin{align}
L_{\alpha_1 \ldots \alpha_k} &= \bigwedge_{i=1}^{k} P_i \\
X_{\alpha_1 \ldots \alpha_k,\beta_1 \ldots \beta_{k-1}} &= \bigvee_{i=1}^k z_{\alpha_i} \delta_{\alpha_1 \ldots \alpha_{i-1} \alpha_{i+1} \ldots \alpha_k, \beta_1 \ldots \beta_{k-1}} \\
v_{\alpha_1 \ldots \alpha_k} &= \bigwedge_{i=1}^{k} z_{\alpha_i}
\end{align}

for all $\alpha_1, \ldots, \alpha_k$ and $\beta_1, \ldots, \beta_{k-1}$ satisfying $0 \leq \alpha_1 < \ldots < \alpha_k < n$ and $0 \leq \beta_1 < \ldots < \beta_{k-1} < n$ where $k$ satisfies $1 \leq k \leq d$, and where $\delta$ is defined as follows:

\begin{equation}
\delta_{\alpha_1 \ldots \alpha_{i-1} \alpha_{i+1} \ldots \alpha_k, \beta_1 \ldots \beta_{k-1}} =
\begin{cases}
1 &\text{if $\alpha_1 = \beta_1, \ldots, \alpha_{i-1} = \beta_{i-1}, \alpha_{i+1} = \beta_{i}, \ldots, \alpha_k = \beta_{k-1}$},\\
0 &\text{otherwise}.
\end{cases}
\end{equation}

All entries of $L, X$ and $\bs v$ are set to zero in the case where $d=0$ for $L$ and in the case where $d\leq 1$ for $X$ and $\bs v$.
Then, the equality

\begin{equation}\label{eq:iff_higer_orders_1}
\big\lvert \bigwedge_{i=1}^{k} P'_{\alpha_i} \big\rvert \equiv \big\lvert \bigwedge_{i=1}^{k} P_{\alpha_i} \big\rvert \pmod{2}
\end{equation}

holds for all $0 \leq \alpha_1 < \ldots < \alpha_k < n$ and $1 \leq k \leq d+1$ if and only if there exists a vector $\bs y'$ and a Boolean $b$ such that $L\bs y \oplus X\bs y' \oplus b\bs v = \bs 0$.

In the case where $k=1$, we have

\begin{equation}
\begin{aligned}
\lvert P'_{\alpha_1} \rvert &= \lvert P_{\alpha_1} \oplus z_{\alpha_1}\bs y \rvert \\
&\equiv \lvert P_{\alpha_1} \rvert + \lvert z_{\alpha_1} \bs y \rvert \pmod{2} \\
&\equiv \lvert P_{\alpha_1} \rvert \pmod{2} \\
\end{aligned}
\end{equation}

because $\lvert \bs y \rvert \equiv 0 \pmod{2}$.
This proves Theorem <ref> in the case where $d=0$, because in such case the equation $L\bs y \oplus X \bs y' \oplus b \bs v = \bs 0$ is always satisfied.
We now prove Theorem <ref> in the case where $d > 0$ and $\bs z$ is satisfying $\lvert \bs z \rvert = 1$.
Without loss of generality we will assume that $z_{\alpha_1} = 1$ and $z_{\alpha_i} = 0$ for all $i$ such that $i \neq 1$, then

\begin{equation}
\begin{aligned}
\big\lvert \bigwedge_{i=1}^{k} P'_{\alpha_i} \big\rvert &= \big\lvert (P_{\alpha_1} \oplus z_{\alpha_1}\bs y) \bigwedge_{i=2}^{k} P_{\alpha_i} \big\rvert \\
&= \big\lvert (P_{\alpha_1} \oplus \bs y) \bigwedge_{i=2}^{k} P_{\alpha_i} \big\rvert \\
&= \big\lvert \bigwedge_{i=1}^{k} P_{\alpha_i} \oplus \bigwedge_{i=2}^{k} P_{\alpha_i} \wedge \bs y \big\rvert \\
&\equiv \big\lvert \bigwedge_{i=1}^{k} P_{\alpha_i} \big\rvert + \big\lvert \bigwedge_{i=2}^{k} P_{\alpha_i} \wedge \bs y \big\rvert \pmod{2} \\
\end{aligned}
\end{equation}

where $\alpha_1, \ldots, \alpha_k$ are satisfying $0 \leq \alpha_1 < \ldots < \alpha_k < n$ and $k$ is satisfying $1 \leq k \leq d+1$.
Therefore, proving Theorem <ref> in the case where $d > 0$ and $\lvert \bs z \rvert = 1$ can be done by showing that there exists a vector $\bs y'$ and a Boolean $b$ such that $L\bs y \oplus X\bs y' \oplus b \bs v = \bs 0$ if and only if the equation

\begin{equation}
\begin{aligned}\label{eq:iff_higer_orders_0}
\big\lvert \bigwedge_{i=2}^{k} P_{\alpha_i} \wedge \bs y \big\rvert &\equiv 0 \pmod{2} \\
\end{aligned}
\end{equation}

holds for all $\alpha_2, \ldots, \alpha_k$ satisfying $\alpha_1 < \alpha_2 < \ldots < \alpha_{k} < n$ and $k$ satisfying $2 \leq k \leq d+1$.
By definition we have $X_{\alpha_2\ldots\alpha_{k}} = \bs 0$ and $v_{\alpha_2\ldots\alpha_{k}} = 0$ for all $\alpha_2, \ldots, \alpha_k$ satisfying $0 \leq \alpha_1 < \ldots < \alpha_{k} < n$ and $k$ satisfying $2 \leq k \leq d+1$.

\begin{equation}
L_{\alpha_2\ldots\alpha_{k}} \bs y \oplus X_{\alpha_2\ldots\alpha_{k}} \bs y' \oplus b v_{\alpha_2\ldots\alpha_{k}} = L_{\alpha_2\ldots\alpha_{k}} \bs y = \big\lvert \bigwedge_{i=2}^{k} P_{\alpha_i} \wedge \bs y \big\rvert \pmod{2}
\end{equation}

and so if $L\bs y \oplus X\bs y' \oplus b \bs v = \bs 0$ then $L_{\alpha_2\ldots\alpha_{k}}\bs y = 0$ which imply that Equation <ref> is satisfied.
Conversely, if Equation <ref> is satisfied then we must have

\begin{equation}
L_{\alpha_2\ldots\alpha_{k}}\bs y \oplus X_{\alpha_2\ldots\alpha_{k}}\bs y' \oplus b \bs v_{\alpha_2\ldots\alpha_{k}} = L_{\alpha_2\ldots\alpha_{k}} = \bs 0
\end{equation}

for all $\alpha_2, \ldots, \alpha_k$ satisfying $\alpha_1 < \alpha_2 < \ldots < \alpha_{k} < n$ and $k$ satisfying $2 \leq k \leq d+1$.
Moreover we have

\begin{equation}
L_{\alpha_1}\bs y \oplus X_{\alpha_1}\bs y' \oplus b v_{\alpha_1} = L_{\alpha_1}\bs y \oplus b
\end{equation}

because $X_{\alpha_1} = \bs 0$ and $v_{\alpha_1} = 1$.

\begin{equation}
L_{\alpha_1\ldots\alpha_{k}}\bs y \oplus X_{\alpha_1\ldots\alpha_{k}}\bs y' \oplus b v_{\alpha_1\ldots\alpha_{k}} = L_{\alpha_1\ldots\alpha_{k}}\bs y \oplus y'_{\alpha_2 \ldots \alpha_{k}}
\end{equation}

for all $\alpha_2, \ldots, \alpha_k$ satisfying $\alpha_1 < \alpha_2 < \ldots < \alpha_{k} < n$ and $k$ satisfying $2 \leq k \leq d$, because $v_{\alpha_1 \ldots \alpha_{k}} = 0$ and

\begin{equation}\label{eq:iff_higer_orders_2}
X_{\alpha_1 \ldots \alpha_k}\bs y' = \bigoplus_{i=1}^k z_{\alpha_i} y'_{\alpha_1 \ldots \alpha_{i-1} \alpha_{i+1} \ldots \alpha_{k}}
\end{equation}

for all $\alpha_1, \ldots, \alpha_k$ satisfying $0 \leq \alpha_1 < \ldots < \alpha_{k} < n$ and $k$ satisfying $2 \leq k \leq d$ and where the rows of $\bs y'$ are labelled the same way as the columns of $X$.
Let $\bs y'$ be such that $y'_{\alpha_2 \ldots \alpha_k} = L_{\alpha_1\ldots\alpha_{k}}\bs y$ for all $\alpha_2, \ldots, \alpha_k$ satisfying $\alpha_1 < \alpha_2 < \ldots < \alpha_{k} < n$ and $k$ satisfying $2 \leq k \leq d$ and $b$ be such that $b = L_{\alpha_1}\bs y$, then we have $L\bs y \oplus X\bs y' \oplus b\bs v = \bs 0$.
Thus, we proved that

\begin{equation}
L\bs y \oplus X\bs y' \oplus b\bs v = \bs 0 \iff \big\lvert \bigwedge_{i=2}^{k} P_{\alpha_i} \wedge \bs y \big \rvert \equiv 0 \pmod{2}
\end{equation}

for all $\alpha_2, \ldots, \alpha_k$ satisfying $\alpha_1 < \alpha_2 < \ldots < \alpha_{k} < n$ and $k$ satisfying $2 \leq k \leq d+1$, and so Theorem <ref> is true in the case where $\lvert \bs z \rvert = 1$.

Let $B$ be a full rank binary matrix of size $n\times n$ which represents a change of basis, and let $\tilde{L}, \tilde{X}$ and $\tilde{\bs v}$ be constructed in the same way as $L, X$ and $\bs v$ but with respect to $BP$ and $B\bs z$.
Notice that if Equation <ref> is satisfied for $P$ and $\bs z$ then it is also satisfied for $BP$ and $B\bs z$.
The proof of Theorem <ref> can then be completed by proving that if there exists a vector $\bs y'$ and a Boolean $b$ such that $L\bs y \oplus X\bs y' \oplus b \bs v = \bs 0$ then there exists a vector $\tilde{\bs y}'$ and a Boolean $\tilde{b}$ such that $\tilde{L}\bs y \oplus \tilde{X}\tilde{\bs y}' \oplus \tilde{b}\tilde{\bs v} = \bs 0$.
Suppose this proposition to be true and let $B$ be such that $\lvert B\bs z \rvert = 1$.
Then, as previously demonstrated, Equation <ref> holds for $BP$ and $B\bs z$ if and only if there exists a vector $\tilde{\bs y}'$ and a Boolean $\tilde{b}$ such that $\tilde{L}\bs y \oplus \tilde{X}\tilde{\bs y}' \oplus \tilde{b}\tilde{\bs v} = \bs 0$, which would imply that there exists a vector $\bs y'$ and a Boolean $b$ such that $L\bs y \oplus X\bs y' \oplus b\bs v = \bs 0$ if and only if Equation <ref> holds for $P$ and $\bs z$.

Let $\tilde{\bs z}$ and $\tilde{P}$ be such that $\tilde{z}_{\alpha_1} = z_{\alpha_1} \oplus z_{\alpha_2}$ and $\tilde{P}_{\alpha_1} = P_{\alpha_1} \oplus P_{\alpha_2}$ for some fixed $\alpha_1, \alpha_2$ satisfying $\alpha_1 \neq \alpha_2$ and $\tilde{z}_\beta = z_\beta$, $\tilde{P}_\beta = P_\beta$ for all $\beta$ such that $\beta \neq \alpha_1$.
Without loss of generality, we will assume that $\alpha_1 = 0$ and $\alpha_2 = 1$.
Indeed, all the other possibilities can be reduced to this case simply by permuting the rows of the parity table $P$ and the vector $\bs z$.
Let $\tilde{L}$, $\tilde{X}$ and $\tilde{\bs v}$ be constructed in the same way as $L$, $X$ and $\bs v$ but with respect to $\tilde{P}$ and $\tilde{\bs z}$.
We will now prove that if there exists a vector $\bs y'$ and a Boolean $b$ such that $L\bs y \oplus X\bs y' \oplus b\bs v = \bs 0$ then $\tilde{L}\bs y \oplus \tilde{X} \tilde{\bs y}' \oplus b \tilde{\bs v} = \bs 0$, where $\tilde{\bs y}'$ satisfies $\tilde{y}'_{\alpha_1\alpha_3\ldots\alpha_k} = y'_{\alpha_1\alpha_3\ldots\alpha_k} \oplus y'_{\alpha_2\ldots\alpha_k}$ and $\tilde{y}'_{\alpha_3\ldots\alpha_k} = y'_{\alpha_3\ldots\alpha_k}$ such that $\alpha_i \neq \alpha_1$ for all $i$.
By using Equation <ref> we can deduce that

\begin{equation}\label{eq:iff_higer_orders_3}
\tilde{X}_{\alpha_3\ldots \alpha_k}\tilde{\bs y}' = \bigoplus_{i=3}^{k} \tilde{z}_{\alpha_i} \tilde{y}'_{\alpha_3\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} = \bigoplus_{i=3}^{k} z_{\alpha_i} y'_{\alpha_3\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} = X_{\alpha_3\ldots \alpha_k}\bs y'
\end{equation}

for all $\alpha_3, \ldots, \alpha_k$ satisfying $\alpha_2 < \alpha_3 < \ldots < \alpha_k < n$ and $k$ satisfying $3 \leq k \leq d+2$, and where the rows of $\tilde{y}'$ and $y'$ are labelled in the same as the columns of $\tilde{X}$ and $X$.
And from the definitions of the vectors $\tilde{\bs v}$ and $\bs v$, we can deduce that

\begin{equation}\label{eq:iff_higer_orders_4}
\tilde{v}_{\alpha_3\ldots\alpha_k} = \bigwedge_{i=3}^k \tilde{z}_{\alpha_i} = \bigwedge_{i=3}^k z_{\alpha_i} = v_{\alpha_3\ldots\alpha_k}
\end{equation}

for all $\alpha_3, \ldots, \alpha_k$ satisfying $\alpha_2 < \alpha_3 < \ldots < \alpha_k < n$ and $k$ satisfying $3 \leq k \leq d+2$.
Equations <ref> and <ref> imply that

\begin{equation}
\tilde{L}_{\alpha_3\ldots \alpha_k}\bs y \oplus \tilde{X}_{\alpha_3\ldots \alpha_k} \tilde{\bs y}' \oplus b\tilde{v}_{\alpha_3\ldots \alpha_k} = L_{\alpha_3\ldots \alpha_k}\bs y \oplus X_{\alpha_3\ldots \alpha_k}\bs y' \oplus b v_{\alpha_3\ldots \alpha_k} = 0
\end{equation}

for all $\alpha_3, \ldots, \alpha_k$ satisfying $\alpha_2 < \alpha_3 < \ldots < \alpha_k < n$ and $k$ satisfying $3 \leq k \leq d+2$.
Furthermore, we have

\begin{equation}
\begin{aligned}
\tilde{L}_{\alpha_1\ldots\alpha_k}\bs y &\equiv \big\lvert \bigwedge_{i=1}^{k} \tilde{P}_{\alpha_i} \wedge \bs y \big\rvert \pmod{2} \\
&\equiv \big\lvert (P_{\alpha_1} \oplus P_{\alpha_2}) \wedge \bigwedge_{i=2}^k P_{\alpha_k} \wedge \bs y \big\rvert \pmod{2} \\
&\equiv \big\lvert \bigwedge_{i=1}^k P_{\alpha_k} \wedge \bs y \big\rvert + \big\lvert \bigwedge_{i=2}^k P_{\alpha_k} \wedge \bs y \big\rvert\pmod{2} \\
&\equiv L_{\alpha_1\ldots\alpha_k}\bs y + L_{\alpha_2\ldots\alpha_k}\bs y \pmod{2} \\
&\equiv X_{\alpha_1\ldots\alpha_k}\bs y' + b v_{\alpha_1\ldots\alpha_k} + X_{\alpha_2\ldots\alpha_k}\bs y' + b v_{\alpha_2\ldots\alpha_k} \pmod{2} \\
\end{aligned}
\end{equation}

for all $\alpha_3, \ldots, \alpha_k$ satisfying $\alpha_2 < \alpha_3 < \ldots < \alpha_k < n$ and $k$ satisfying $3 \leq k \leq d$, which entails

\begin{equation}
\begin{aligned}
&\tilde{L}_{\alpha_1\ldots\alpha_k}\bs y \oplus \tilde{X}_{\alpha_1\ldots\alpha_k}\tilde{\bs y}' \oplus b \tilde{v}_{\alpha_1\ldots\alpha_k} \\
&= X_{\alpha_1\ldots\alpha_k}\bs y' \oplus X_{\alpha_2\ldots\alpha_k}\bs y' \oplus \tilde{X}_{\alpha_1\ldots\alpha_k}\tilde{\bs y}'\oplus b v_{\alpha_1\ldots\alpha_k} \oplus b v_{\alpha_2\ldots\alpha_k} \oplus b \tilde{v}_{\alpha_1\ldots\alpha_k} \\
&= \bigoplus_{i=1}^{k} z_{\alpha_i} y'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus \bigoplus_{i=2}^{k} z_{\alpha_i} y'_{\alpha_2\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus \bigoplus_{i=1}^{k} \tilde{z}_{\alpha_i} \tilde{y}'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \\
&\quad \oplus b \bigwedge_{i=1}^k z_{\alpha_i} \oplus b \bigwedge_{i=2}^k z_{\alpha_i} \oplus b \bigwedge_{i=1}^k \tilde{z}_{\alpha_i}\\
&= \bigoplus_{i=1}^{k} z_{\alpha_i} y'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus \bigoplus_{i=2}^{k} z_{\alpha_i} y'_{\alpha_2\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus (z_{\alpha_1} \oplus z_{\alpha_2})y'_{\alpha_2\ldots\alpha_k} \\
&\quad \oplus \bigoplus_{i=2}^{k} z_{\alpha_i} \tilde{y}'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus b \bigwedge_{i=1}^k z_{\alpha_i} \oplus b \bigwedge_{i=2}^k z_{\alpha_i} \oplus b \bigwedge_{i=2}^k z_{\alpha_i} \wedge (z_{\alpha_1} \oplus z_{\alpha_2}) \\
&= \bigoplus_{i=2}^{k} z_{\alpha_i} (y'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus y'_{\alpha_2\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k}) \oplus z_{\alpha_2} y'_{\alpha_2\ldots\alpha_k} \\
& \quad \oplus \bigoplus_{i=2}^{k} z_{\alpha_i} \tilde{y}'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus b \bigwedge_{i=1}^k z_{\alpha_i} \oplus b \bigwedge_{i=2}^k z_{\alpha_i} \oplus b \bigwedge_{i=1}^k z_{\alpha_i} \oplus b \bigwedge_{i=2}^k z_{\alpha_i} \\
&= \bigoplus_{i=2}^{k} z_{\alpha_i} (y'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus y'_{\alpha_2\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k}) \oplus z_{\alpha_2} y'_{\alpha_2\ldots\alpha_k} \oplus z_{\alpha_2} \tilde{y}'_{\alpha_1\alpha_{3}\ldots\alpha_k} \\
&\quad \oplus \bigoplus_{i=3}^{k} z_{\alpha_i} \tilde{y}'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \\
&= \bigoplus_{i=2}^{k} z_{\alpha_i} (y'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus y'_{\alpha_2\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k}) \oplus z_{\alpha_2} y'_{\alpha_2\ldots\alpha_k} \oplus z_{\alpha_2} (y'_{\alpha_1\alpha_{3}\ldots\alpha_k} \oplus y'_{\alpha_2\ldots\alpha_k} \\
&\quad \oplus y'_{\alpha_{3}\ldots\alpha_k}) \oplus \bigoplus_{i=3}^{k} z_{\alpha_i} (y'_{\alpha_1\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k} \oplus y'_{\alpha_2\ldots\alpha_{i-1}\alpha_{i+1}\ldots\alpha_k})\\
&= 0.
\end{aligned}
\end{equation}

Finally, we have

\begin{equation}
\begin{aligned}
\tilde{L}_{\alpha_1}\bs y \oplus \tilde{X}_{\alpha_1}\tilde{\bs y}' \oplus b\tilde{v}_{\alpha_1}
&= \tilde{L}_{\alpha_1}\bs y \oplus b\tilde{z}_{\alpha_1} \\
&= \tilde{L}_{\alpha_1}\bs y \oplus b(z_{\alpha_1} \oplus z_{\alpha_2}) \\
&\equiv \lvert \tilde{P}_{\alpha_1} \wedge \bs y \rvert + bv_{\alpha_1} + bv_{\alpha_2} \pmod{2} \\
&\equiv \lvert (P_{\alpha_1} \oplus P_{\alpha_2}) \wedge \bs y \rvert + bv_{\alpha_1} + bv_{\alpha_2} \pmod{2} \\
&\equiv \lvert P_{\alpha_1} \wedge \bs y \rvert + \lvert P_{\alpha_2} \wedge \bs y \rvert + bv_{\alpha_1} + bv_{\alpha_2} \pmod{2} \\
&= L_{\alpha_1} \oplus bv_{\alpha_1} \oplus L_{\alpha_2} \oplus bv_{\alpha_2} \\
&= 0.
\end{aligned}
\end{equation}

Thus, we proved that

\begin{equation}
L\bs y \oplus X\bs y' \oplus b \bs v = \bs 0 \implies \tilde{L}\bs y \oplus \tilde{X}\tilde{\bs y}' \oplus b \tilde{v} = \bs 0
\end{equation}

which concludes the proof of Theorem <ref>.

§ $R_Z(\PI/2^D)$-COUNT UPPER BOUND IN UNIVERSAL GATE SETS

In this section we demonstrate an upper bound achievable in polynomial time for the $R_Z(\pi/2^d)$-count within a $\{\mathrm{CNOT},$ $H$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ circuit where $d$ is a non-negative integer.
Note that it corresponds to an upper bound for the number of $T$ gates in a Clifford$+T$ circuit in the case where $d=4$.
Our proof for the upper bound will partially rest on the following lemma, which has already been proven with a different approach in Reference <cit.> for the case where $d=2$.

Let $d$ be a non-negative integer and let $U \in \mathcal{D}_{d+1}$ act on $n$ qubits where $\mathcal{D}_{d+1}$ is the diagonal subgroup of the $(d+1)$th level of the Clifford hierarchy.
Then $U$ can be decomposed into two unitary gates $U = \tilde{U} \tilde{U}'$ which can be found in polynomial time and such that $\tilde{U} \in \mathcal{D}_{d+1}$ is acting on $n-1$ qubits.
Furthermore, we can find a $\{\mathrm{CNOT},$ $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ circuit implementing $\tilde{U}' \in \mathcal{D}_{d+1}$ with no more than

\begin{equation}
\sum_{i=0}^{d-1} {n-1 \choose i} + 1
\end{equation}

$R_Z(\pi/2^d)$ gates in polynomial time.

Let $P$ be a parity table associated with an implementation of $U$, let $\beta$ be the qubit on which $\tilde{U}$ is not acting, and let $\tilde{P}'$ be a parity table such that

\begin{equation}\label{eq:eq_1_d_upper_bound}
\big\lvert \bigwedge_{i=1}^d \tilde{P}'_{\alpha_i} \big\rvert \equiv \big\lvert \bigwedge_{i=1}^d P_{\alpha_i} \wedge P_\beta \big\rvert \pmod{2}
\end{equation}

for all $\alpha_1, \ldots, \alpha_d$ satisfying $0 \leq \alpha_1 \leq \ldots \leq \alpha_d < n$ and $\alpha_i \neq \beta$ for all $i$.
The row $\beta$ of $\tilde{P}'$ can be ignored, and, as stated by Theorem <ref>, we can optimize the number of columns of $\tilde{P}'$ such that  it is lower or equal to

\begin{equation}
\sum_{i=0}^{d-1} {n-1 \choose i}
\end{equation}

and still satisfies Equation <ref>.
If we add a null column to $\tilde{P}'$ in the case where the number of columns of $\tilde{P}'$ is not equal to $\lvert P_\beta \rvert \pmod{2}$, and then set $\tilde{P}'_\beta = \bs 1$, then we have

\begin{equation}\label{eq:eq_2_d_upper_bound}
\big\lvert \bigwedge_{i=1}^d \tilde{P}'_{\alpha_i} \wedge \tilde{P}'_\beta \big\rvert = \big\lvert \bigwedge_{i=1}^d \tilde{P}'_{\alpha_i} \big\rvert \equiv \big\lvert \bigwedge_{i=1}^d P_{\alpha_i} \wedge P_\beta \big\rvert \pmod{2}
\end{equation}

for all $0 \leq \alpha_1 \leq \ldots \leq \alpha_d < n$ and the number of columns of $\tilde{P}$ is at most

\begin{equation}
\sum_{i=0}^{d-1} {n-1 \choose i} + 1.
\end{equation}

Finally, we can easily find a parity table $\tilde{P}$ such that $\tilde{P}_\beta = \bs 0$ and

\begin{equation}\label{eq:eq_3_d_upper_bound}
\big\lvert \bigwedge_{i=1}^{d+1} \tilde{P}_{\alpha_i} \big\rvert \equiv \big\lvert \bigwedge_{i=1}^{d+1} \tilde{P}'_{\alpha_i} \big\rvert + \big\lvert \bigwedge_{i=1}^{d+1} P_{\alpha_i} \big\rvert \pmod{2}
\end{equation}

for all $0 \leq \alpha_1 \leq \ldots \leq \alpha_d < n$, which implies its associated unitary gate $\tilde{U}$ doesn't act on qubit $\beta$.
Thus, $\tilde{U}$ and $\tilde{U}'$ are forming a decomposition of $U$ up to an operator $V$ implementable over $\{\mathrm{CNOT}$, $R_Z(2\pi/2^d)\}$ gate set: $U = \tilde{U}\tilde{U}'V$, where $\tilde{U}'$ is the unitary gate associated with the parity table $\tilde{P}'$.
The operator $V$ can then be merged with the unitary $\tilde{U}'$ to comply with the decomposition of $U$ established in Lemma <ref>.

Based on Lemma <ref> and Theorem <ref>, we can prove the following theorem which provides an upper bound for the number of $R_Z(\pi/2^d)$ gates in a $\{\mathrm{CNOT}$, $H$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ circuit.
Note that in the case where $d = 2$, we get an upper bound of $(n+1)(n+2h)/2 + 1$ for the number of $T$ gates in a Clifford$+T$ circuit, which is significantly better than the previously best-known ancilla-free upper bound of $\mathcal{O}(n^2h)$ <cit.>, where $n$ is the number of qubits and $h$ is the number of Hadamard gates in the circuit.

Let $U$ be a unitary gate acting on $n$ qubits and implementable over the $\{\mathrm{CNOT}$, $H$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ gate set where $d$ is a non-negative integer and with $h$ internal Hadamard gates such that

\begin{equation}\label{eq:d_u_decomposition}
U = C_f U_0 \prod_{i=1}^h \left[ H_{\alpha_i} U_i \right] C_{init}
\end{equation}

where $C_f, C_{init}$ are implementable over the $\{\mathrm{CNOT}$, $H$, $R_Z(2\pi/2^d)\}$ gate set, $H_{\alpha_i}$ denotes the Hadamard gate applied on some qubit $\alpha_i$, and $U_0,\ldots, U_h$ are implementable over the $\{\mathrm{CNOT}$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ gate set.
Then we can find a implementation of $U$ over the $\{\mathrm{CNOT}$, $H$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ gate set, in polynomial time, and such that the number of $R_Z(\pi/2^d)$ gates is lower or equal to

\begin{equation}
\sum_{i=0}^{d} {n \choose i} + h\left(\sum_{i=0}^{d-1} {n-1 \choose i} + 1\right).
\end{equation}

As stated by Lemma <ref>, the unitary $U_h$ can be decomposed into two unitary gates: $U_h = \tilde{U}_h \tilde{U}'_h$, such that $\tilde{U}_h$ is not acting on qubit $\alpha_h$, which imply that $\tilde{U}_h$ commutes with the $h$th internal Hadamard gate.
Let $U_i \tilde{U}_{i+1} = \tilde{U}_i \tilde{U}'_i$ for all $i$ satisfying $0 < i < h$ where $\tilde{U}_i \tilde{U}'_i$ is a decomposition of $U_i \tilde{U}'_{i+1}$ such as given by Lemma <ref> where $\tilde{U}'_i$ is not acting on qubit $\alpha_i$.
Then Equation <ref> is equivalent to

\begin{equation}
U = C_f U_0\tilde{U}_1 \prod_{i=1}^h \left[ H_{\alpha_i} \tilde{U}'_i \right] C_{init}
\end{equation}

and, as stated by Lemma <ref>, we can find an implementation of $\tilde{U}'_i$ in polynomial time and with no more than

\begin{equation}
\sum_{k=0}^{d-1} {n-1 \choose k} + 1
\end{equation}

$R_Z(\pi/2^d)$ gates, for all $i$ satisfying $0 < i \leq h$.
Moreover, as stated by Theorem <ref>, we can find an implementation of $U_0 \tilde{U}_1$ in polynomial time and with no more than

\begin{equation}
\sum_{i=0}^{d} {n \choose i}
\end{equation}

$R_Z(\pi/2^d)$ gates.
Thus, $U$ can be implemented in polynomial time over the $\{\mathrm{CNOT}$, $H$, $R_Z(\pi/2^d)$, $R_Z(2\pi/2^d)\}$ gate set and with a number of $R_Z(\pi/2^d)$ gates that is lower or equal to

\begin{equation}
\sum_{i=0}^{d} {n \choose i} + h\left(\sum_{i=0}^{d-1} {n-1 \choose i} + 1\right).
\end{equation}

§ CONCLUSION

We presented several polynomial-time algorithms for reducing the number of $T$ gates in a Clifford$+T$ circuit.
Benchmarks show that our algorithms consistently achieve the lowest execution times and provide the best $T$-count reduction on almost all the quantum circuits evaluated when compared to state-of-the-art $T$-count optimizers.
As such, our algorithms not only achieve state-of-the-art $T$-count reduction but also offer much greater scalability, thereby allowing efficient $T$-count optimization on larger quantum circuits.

We proved that our algorithms are producing a circuit in which the number of $T$ gates is upper bounded by $(n^2 + n)/2 + 1$ when they are executed on a Hadamard-free circuit, where $n$ is the number of qubits.
It has been shown that there exists an asymptotic upper bound of $n^2/2 - 1$ <cit.>.
The question of whether or not there exists a polynomial-time algorithm satisfying this asymptotic upper bound remains open.

We also demonstrated how the number of $T$ gates in a Clifford$+T$ circuit can be optimized such that it is lower or equal to $(n+1)(n+2h)/2 + 1$ where $n$ is the number of qubits and $h$ is the number of internal Hadamard gates in the circuit, without using any ancillary qubits and in polynomial time.
This reinforces the interdependence between the problems of optimizing the number of internal Hadamard gates and optimizing the number of $T$ gates.
As part of future work, it would be beneficial to more clearly determine the roles that internal Hadamard gates and ancillary qubits have in the $T$-count optimization problem.

§ ACKNOWLEDGMENTS

We acknowledge funding from the Plan France 2030 through the projects NISQ2LSQ ANR-22-PETQ-0006 and EPIQ ANR-22-PETQ-007.

§ PROOF OF THEOREM <REF>

Equation <ref> holds if and only if the two following equations are true:

\begin{align}
\!g(\bs a, \bs b, \bs c) &= g(\bs a_L \oplus \bs a_R, \bs b_L \oplus \bs b_R, \bs c_R) \oplus g(\bs a_R, \bs b_R, \bs c_R) \oplus g(\bs a_L, \bs b_L, \bs c_L \oplus \bs c_R) \oplus h(\bs a_L, \bs b_L, \bs c_R)  \label{eq:gf_2} \\
\!h(\bs a, \bs b, \bs c') &= h(\bs a_L \oplus \bs a_R, \bs b_L \oplus \bs b_R, \bs c'_L) \oplus g(\bs a_R, \bs b_R, \bs c'_L) \oplus h(\bs a_R, \bs b_R, \bs c'_L \oplus \bs c'_R) \oplus h(\bs a_L, \bs b_L, \bs c'_L) \label{eq:gf_1}
\end{align}

We first prove Equation <ref>:

\begin{equation}
\begin{aligned}
&h(\bs a_L \oplus \bs a_R, \bs b_L \oplus \bs b_R, \bs c'_L) \oplus g(\bs a_R, \bs b_R, \bs c'_L) \oplus h(\bs a_R, \bs b_R, \bs c'_L \oplus \bs c'_R) \oplus h(\bs a_L, \bs b_L, \bs c'_L) \\
&= \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} (a_j + a_{\frac{n}{2} + j}) (b_{\frac{n}{2} + i - j} + b_{n + i - j}) c'_i + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_{\frac{n}{2} + j} b_{\frac{n}{2} + i - j} c'_i \\
&+ \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_{\frac{n}{2} + j} b_{n + i - j} (c'_i + c'_{\frac{n}{2} + i}) + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_j b_{\frac{n}{2} + i - j} c'_i \pmod{2} \\
&= \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_j b_{n + i - j} c'_i  + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_{\frac{n}{2} + j} b_{\frac{n}{2} + i - j} c'_i + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_{\frac{n}{2} + j} b_{\frac{n}{2} + i - j} c'_i \\
&+ \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_{\frac{n}{2} + j} b_{n + i - j} c'_{\frac{n}{2} + i} \pmod{2} \\
&= \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_j b_{n + i - j} c'_i  + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=\frac{n}{2}+i+1}^{n-1} a_{j} b_{n + i - j} c'_i  + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=\frac{n}{2}}^{\frac{n}{2} + i} a_{j} b_{n + i - j} c'_i \\
&+ \sum_{i=\frac{n}{2}}^{n-2} \sum_{j=i+1}^{n-1} a_{j} b_{n + i - j} c'_{i} \pmod{2} \\
&= \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_j b_{n + i - j} c'_i  + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=\frac{n}{2}+i+1}^{n-1} a_{j} b_{n + i - j} c'_i + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=\frac{n}{2}}^{\frac{n}{2} + i} a_{j} b_{n + i - j} c'_i \\
&+  \sum_{i=\frac{n}{2}-1}^{\frac{n}{2}-1} \sum_{j=i+1}^{n-1} a_{j} b_{n + i - j} c'_i + \sum_{i=\frac{n}{2}}^{n-2} \sum_{j=i+1}^{n-1} a_{j} b_{n + i - j} c'_{i} \pmod{2} \\
&= \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_j b_{n + i - j} c'_i  + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=\frac{n}{2}+i+1}^{n-1} a_{j} b_{n + i - j} c'_i + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=\frac{n}{2}}^{\frac{n}{2} + i} a_{j} b_{n + i - j} c'_i \\
&+ \sum_{i=\frac{n}{2}-1}^{n-2} \sum_{j=i+1}^{n-1} a_{j} b_{n + i - j} c'_{i} \pmod{2} \\
&= \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}+i} a_j b_{n + i - j} c'_i  + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=\frac{n}{2}+i+1}^{n-1} a_{j} b_{n + i - j} c'_i + \sum_{i=\frac{n}{2}-1}^{n-2} \sum_{j=i+1}^{n-1} a_{j} b_{n + i - j} c'_{i} \pmod{2} \\
&= \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{n-1} a_{j} b_{n + i - j} c'_i + \sum_{i=\frac{n}{2}-1}^{n-2} \sum_{j=i+1}^{n-1} a_{j} b_{n + i - j} c'_{i} \pmod{2} \\
&= \sum_{i=0}^{n-2} \sum_{j=i+1}^{n-1} a_{j} b_{n + i - j} c'_{i} \pmod{2} \\
& = h(\bs a, \bs b, \bs c').
\end{aligned}
\end{equation}

We now prove Equation <ref>:

\begin{equation}
\begin{aligned}
&g(\bs a_L \oplus \bs a_R, \bs b_L \oplus \bs b_R, \bs c_R) \oplus g(\bs a_R, \bs b_R, \bs c_R) \oplus g(\bs a_L, \bs b_L, \bs c_L \oplus \bs c_R) \oplus h(\bs a_L, \bs b_L, \bs c_R) \\
&= \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} (a_j + a_{\frac{n}{2} + j}) (b_{i - j} + b_{\frac{n}{2} + i - j}) c_{\frac{n}{2} + i} + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_{\frac{n}{2} + j} b_{\frac{n}{2} + i - j} c_{\frac{n}{2} + i} \\
&+ \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_j b_{i - j} (c_i + c_{\frac{n}{2} + i}) + \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_j b_{\frac{n}{2} + i - j} c_{\frac{n}{2} + i} \pmod{2} \\
&= \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_j b_{\frac{n}{2} + i - j} c_{\frac{n}{2} + i} + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_{\frac{n}{2} + j} b_{i - j} c_{\frac{n}{2} + i} + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_j b_{i - j} c_i \\
&+ \sum_{i=0}^{\frac{n}{2}-2} \sum_{j=i+1}^{\frac{n}{2}-1} a_j b_{\frac{n}{2} + i - j} c_{\frac{n}{2} + i} \pmod{2} \\
&= \sum_{i=\frac{n}{2}}^{n-1} \sum_{j=0}^{i - \frac{n}{2}} a_j b_{i - j} c_{i} + \sum_{i=\frac{n}{2}}^{n-1} \sum_{j=\frac{n}{2}}^{i} a_{j} b_{i - j} c_{i} + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_j b_{i - j} c_i \\
&+ \sum_{i=\frac{n}{2}}^{n-2} \sum_{j=i-\frac{n}{2}+1}^{\frac{n}{2}-1} a_j b_{i - j} c_{i} \pmod{2} \\
&= \sum_{i=\frac{n}{2}}^{n-1} \sum_{j=0}^{\frac{n}{2} - 1} a_j b_{i - j} c_{i} + \sum_{i=\frac{n}{2}}^{n-1} \sum_{j=\frac{n}{2}}^{i} a_{j} b_{i - j} c_{i} + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_j b_{i - j} c_i \pmod{2} \\
&= \sum_{i=\frac{n}{2}}^{n-1} \sum_{j=0}^{i} a_j b_{i - j} c_{i} + \sum_{i=0}^{\frac{n}{2}-1} \sum_{j=0}^{i} a_j b_{i - j} c_i \pmod{2} \\
&= \sum_{i=0}^{n-1} \sum_{j=0}^{i} a_j b_{i - j} c_i \pmod{2} \\
&= g(\bs a, \bs b, \bs c).
\end{aligned}
\end{equation}


# Statistical analysis of inter coding in VVC Test Model (VTM)

Yiqun Liu12, Mohsen Abdoli1, Thomas Guionnet1, Christine Guillemot2, and Aline
Roumy2  1Ateme, Rennes, France  2INRIA, Rennes, France

###### Abstract

The promising improvement in compression efficiency of Versatile Video Coding
(VVC) compared to High Efficiency Video Coding (HEVC) [1] comes at the cost of
a non-negligible encoder side complexity. The largely increased complexity
overhead is a possible obstacle towards its industrial implementation. Many
papers have proposed acceleration methods for VVC. Still, a better
understanding of VVC complexity, especially related to new partitions and
coding tools, is desirable to help the design of new and better acceleration
methods. For this purpose, statistical analyses have been conducted, with a
focus on Coding Unit (CU) sizes and inter coding modes.

###### Index Terms:

Versatile Video Coding, Inter Coding, Rate Distortion Optimization, Complexity
Analysis

##  I Introduction

Standardization of VVC in 2020 has brought significant improvement to the
capacity of video compression in terms of bitrate saving. It offers the 50%
compression efficiency [2] compared to one of the most efficient video
compression standards HEVC. This improvement in VVC is mainly due to newly
adopted coding techniques. Most decoder devices could afford the additional
complexity brought by these novel coding techniques taking into account
current hardware capacities of these devices. Specifically, studies in [3]
have shown that the relative complexity of the decoder of VVC is from 150% to
200% compared to HEVC in different configurations.

On the contrary, it is far from affordable for real time application for VVC
encoder side since industrial encoding applications have strict limitations in
terms of resources and execution time. Tests in [4] in VVC reference software
VVC Test Model (VTM) 7 show that the encoding time of VVC is 5x, 7x, and 37x
times of the encoding of HEVC in the configurations Low-Delay (LD), Random-
Access (RA), and All-Intra (AI), respectively. Hence, it is vital to develop
acceleration algorithms or methods to largely reduce the encoding complexity
while preserving the majority of encoding efficiency. Complexity analysis
papers could help researchers to have a clear understanding of what is
happening inside a VVC encoder (e.g. VTM), and what potentially interests them
for their design process of the acceleration method.

Various studies have contributed to the complexity analysis of VVC. In [5], a
detailed complexity analysis based on VVC intra prediction tools has been
performed. Pakdaman et al. in [3] have broken down the encoding process into
encoding modules such as motion estimation, intra prediction, entropy coding,
etc. and then analyzed the complexity partition of modules in multiple
encoding configurations. [6] reviews complexity aspects of the different
modules of the VVC standard and provide a complexity breakdown of these
modules in a more precise way. In [4], VVC and HEVC are compared in terms of
rate-distortion and complexity analysis. These aforementioned papers present
complexity analysis at the level of encoding modules for inter coding. Our
paper is the first to provide an analysis from CU sizes and coding modes
perspective for inter coding in VVC.

Figure 1: High-level view of the RDO process involving partitioning, test
modes and possible VTM shortcuts.

In this paper, a statistical analysis of the Rate-Distortion Optimization
(RDO) process in inter coding of VTM-15.0 is presented. The main focus is put
on the statistics of two factors: CU sizes and inter coding modes. The goal is
to provide useful information for related works aimed at speeding up inter
coding in VVC. The rest of this paper is organized as follows. Section II
presents a summary of VVC specification in terms of CU sizes and available
coding modes. In Section III, all statistical observations are presented,
which are later analyzed and concluded in Section IV.

##  II RDO of inter coding in VTM

To represent the RDO process in VTM, there exist numerous coding parameters
such as Intra Prediction Mode (IPM) of intra coding, Motion Vector (MV)
representation mode, choice of transform etc.. However, if we ignore these
trivial parameters, the RDO process could be described as the search for the
best trade-off between bit rate and distortion. More precisely, this search is
executed on different coding modes of different CU sizes. Therefore, in this
paper, the statistics of CU sizes and coding modes are jointly considered.

As presented in Fig.1, various CU sizes are the result of partitioning in the
RDO process. The partitioning consists of splitting the CU of size
128$\times$128 recursively by five split modes, namely Quaternary Tree (QT)
split, Horizontal Binary Tree (HBT) split, Vertical Binary Tree (VBT) split,
Horizontal Tenary Tree (HTT) split, Vertical Tenary Tree (VTT) split. Compared
to codec HEVC in which only QT is available for partitioning, the added
directional splits give rise to a larger variety of CU sizes. In VVC CU size
is authorised if its widths and heights are any power of two between 4 and
128, except for sizes 128$\times$4, 128$\times$8, 128$\times$16 and
128$\times$32\. It is worth noting that the same CU size could be obtained by
different series of split modes.

For each CU, 12 coding modes are available. Two of these modes, namely hash-
based inter prediction and palette modes are not enabled in JVET Common Test
Condition (CTC). Hence, they are not included in the analysis of this paper.
In the inter prediction of VVC, motion compensation is executed after motion
estimation. Subsequently, the residuals and MV information need to be
transmitted. Depending on the transmission of MV and residuals, three coding
modes are available: Advanced Motion Vector Prediction (AMVP) [7] mode,
_Merge_ mode, and _Skip_ mode. For the purpose of simplification, we refer to
AMVP as _Reg_ mode for the remainder of the paper. Before transmitting MV s, a
candidate list of MV s is constructed based on the spatial and temporal
neighboring CU s by exploiting the correlations of MV s between them. Four
types of inter prediction data can be signaled, including the index of
reference frame (i.e. Ref Frame Idx), the index of the best MV candidate (i.e.
MV Cand Idx), the difference between the best candidate and MV determined by
motion estimation (i.e. Motion Vector Difference (MVD)), and residuals. Tab.I
presents the signaled data types for _Reg_ , _Merge_ and _Skip_.

TABLE I: Data types to transmit for motion data coding | Ref Frame Idx | MV Cand Idx | MVD | Residuals
---|---|---|---|---
_Reg_ | ✓ | ✓ | ✓ | ✓
_Merge_ | X | ✓ | X | ✓
_Skip_ | X | ✓ | X | X

In addition to the _Affine_ mode and the _Intra_ mode, two novel coding modes
are available in VVC. For CU s coded in merge mode, Combined Intra-Inter
Prediction (CIIP) [8] combines the inter prediction and the intra prediction
to form a final prediction. The Geometric Partitioning Mode, denoted as _Geo_
, is designed to better predict moving objects in video. Additionally, _Geo_
is conventionally coded with _Merge_ mode. For inter coding configuration,
coding modes of _Reg_ , _Merge_ and _Skip_ could combine with _Affine_ ,
_CIIP_ , and _Geo_ , which results in a total of 10 coding modes: _Intra_ ,
_Reg_ , _Merge_. _Skip_ , _Affine_ , _AffineMerge_ , _AffineSkip_ , _GeoMerge_
, _GeoSkip_ , _CiipMerge_. Statistics of these coding modes are further
collected and analyzed in the following part.

Although VVC is computationally expensive, the Joint Video Exploration Team
(JVET) group has already adopted various shortcuts or conditional early exits
as presented in [9] for the VTM. We have deactivated existing shortcuts in
VTM-15.0 and evaluated its performance. As a result, the complexity increases
by 138%. Furthermore, the performance of the tested encoder (i.e. without
shortcuts) is 0.76% better than the reference encoder (i.e. with shortcuts),
in terms of BD-rate. This trade-off might be interpreted as an indicator that
the shortcuts in VTM are efficient in terms of identifying useless tests and
partitioning depths. Many of shortcuts are based on history of the tested
split modes. However, aspects of CU sizes and coding modes are overlooked.

Figure 2: Encoding time in different QPs comparing to QP 22

##  III Statistics

Our main purpose in the following analysis is to find CU sizes and/or coding
modes with a relatively high complexity occupation and a low selection rate in
the RDO process. From the perspective of encoder acceleration, CU sizes or
coding modes with a higher complexity portion and a significantly lower
selection rate are more favorable to the design of acceleration rules based on
CU size/coding mode. The size or mode with larger complexity portion has more
potential in accelerating. Lower selection rate indicates it is less likely to
make wrong decisions when skipping the RDO of current CU size or mode.

All our experiments and analyses are performed on the first 64 frames of the
CTC sequences in the RandomAccess Group Of Picture 32 (RAGOP32) configuration
in which intra frames are excluded. Exceptionally, Fig.2 is based on sequences
in Class A, B, E of CTC.

Figure 3: Complexity distribution for CU sizes in QP22 and QP37

The encoding complexity for Chrominance channel only accounts for a small part
comparing to Luminance channel. Thus we focus on Luminance channel in the
remaining of the paper. From a high-level perspective, the encoding complexity
of VTM significantly depends on the selected Quantization Parameter (QP).
Particularly, larger QP values tend to have faster encoding with the VTM. Fig.
2 is obtained by measuring the encoding times of sequences of resolution
2160p, 1080p and 720p in QP 22, 27, 32 and 37. Then the average ratio is
calculated between the encoding time of each QP and that of QP 22 is
calculated. It shows that the encoding time at QP 22 could be five times as
much as at QP 37.

Fig.3 shows the percentage distribution of the encoding time spent on
different CU sizes in QP 22 and QP 37\. In addition to the fact that the
overall encoding time is higher for QP 22, it can be observed that a
relatively higher portion of the time in QP 22 is passed on smaller CU sizes.
This could be partly explained by the existing shortcuts in VTM disallowing
excessively small CU s in QP 37\. We could declare that larger CU sizes are in
general more crucial to speeding up the partitioning process, especially CU
64$\times$64 and 128$\times$128 which take in total from 20% complexity in QP
22 to 30% in QP 37.

In another test, the selection percentages of different CU sizes are
calculated. This metric is defined as the ratio between the total number of
times it is selected and the total number of times a CU size is tested. Fig. 4
shows the values of this metric in QP 22 and 37. As we can see from this
figure, larger CU sizes correspond to larger selection rates compared with
smaller CU s. Another phenomenon worth noting is that the selection rate of
128$\times$128 increase dramatically from 14% in QP 22 to 37% in QP 37.

Figure 4: Selection rate for different CU sizes

Combining the above figures, it is observed that CU sizes such as 16$\times$8,
8$\times$16 and 16$\times$16 are sizes with low selection rate and high
complexity. For example, 16$\times$16 CU s have the same level of complexity,
while its selection rate is half of 32$\times$32 CU s in QP 22.

Figure 5: Pie chart of complexities of inter coding modes

To take one step further in the statistical analysis of inter coding, we
present how different inter coding modes are involved in RDO search. The first
experiment in Fig.5 presents the distribution of the encoding time at inter
coding mode level in QP 22 and QP 37.

In general, _Intra_ , _AffineMerge_ , _AffineSkip_ , _Merge_ , and _Skip_ are
main contributors to encoding time. Fig.6 provides selection rates as the
ratio between the number of selected inter coding modes and number of tested
modes. We could observe that the three modes, namely _AffineMerge_ ,
_AffineSkip_ , and _Merge_ have relatively low selection rates, although they
collectively account for nearly half of the complexity.

Figure 6: Selection rate of inter coding modes Figure 7: Stacked chart of
selected inter modes in different CU sizes

Fig.7 shows the distribution of inter modes for encoded CU s of different
sizes. The fact that the aforementioned three coding modes are less chosen
could also be proved by this figure. We find that the number of _Skip_ is
dominant for most CU sizes and that the number of these three modes is
relatively small, which is consistent with Fig.6. From Fig.7 we also observe
that smaller CU s tend to be encoded with intra mode. Another remark is that
the skip modes (i.e. _Skip_ , _AffineSkip_ , and _GeoSkip_) are more
frequently selected for larger CU s. It is probably because the residual of
larger CU s is more expensive to be encoded. In addition, _Merge_ mode has a
higher chance to be selected in smaller QP which is in contrast to
_AffineSkip_ and _GeoSkip_. For some CU sizes, we could make shortcuts or
conditions for early termination for RDO of inter modes which are rarely
selected to speed up the encoding process, such as the _AffineMerge_ mode with
merely 1.9% selected for CU 128$\times$128.

##  IV Analysis and Conclusion

In this study, complexity analysis of CU sizes and inter coding modes has been
combined with selection rate analysis. From the perspective of CU size, CU
sizes with high complexity generally correspond to a high selection rate. CU
sizes 128$\times$128 and 64$\times$64 are responsible for one-third of the
complexity. In addition, CU sizes such as 16$\times$8, 8$\times$16,
16$\times$16 exhibit relatively low selection rate while requiring a
significant share of the overall complexity. Therefore, these CU sizes are
relevant targets for acceleration algorithms. From a coding mode perspective,
_AffineMerge_ , _AffineSkip_ , and _Merge_ tend to be less likely to be
selected. Thus, dedicated shortcuts to adaptively skip these coding modes
might be promising. Shortcuts on coding modes and partitioning acceleration
method are in different scopes. The former focus on reducing number of CU for
RDO. The latter speeds up RDO for CU of certain sizes. The combination of
these two could lead to a larger speed-up of encoding.

## References

  * [1] Gary J Sullivan, Jens-Rainer Ohm, Woo-Jin Han, and Thomas Wiegand,  “Overview of the high efficiency video coding (hevc) standard,”  IEEE Transactions on circuits and systems for video technology, vol. 22, no. 12, pp. 1649–1668, 2012.
  * [2] B. Bross et al.,  “Overview of the versatile video coding (VVC) standard and its applications,”  IEEE Trans. on Circuits and Systems for Video Technology, vol. 31, no. 10, pp. 3736–3764, 2021.
  * [3] Farhad Pakdaman, Mohammad Ali Adelimanesh, Moncef Gabbouj, and Mahmoud Reza Hashemi,  “Complexity analysis of next-generation vvc encoding and decoding,”  in 2020 IEEE International Conference on Image Processing (ICIP). IEEE, 2020, pp. 3134–3138.
  * [4] Alexandre Mercat, Arttu Mäkinen, Joose Sainio, Ari Lemmetti, Marko Viitanen, and Jarno Vanne,  “Comparative rate-distortion-complexity analysis of vvc and hevc video codecs,”  IEEE Access, vol. 9, pp. 67813–67828, 2021.
  * [5] Mário Saldanha, Gustavo Sanchez, César Marcon, and Luciano Agostini,  “Complexity analysis of vvc intra coding,”  in 2020 IEEE International Conference on Image Processing (ICIP). IEEE, 2020, pp. 3119–3123.
  * [6] Frank Bossen, Karsten Sühring, Adam Wieckowski, and Shan Liu,  “Vvc complexity and software implementation analysis,”  IEEE Transactions on Circuits and Systems for Video Technology, vol. 31, no. 10, pp. 3765–3778, 2021.
  * [7] Wei-Jung Chien, Li Zhang, Martin Winken, Xiang Li, Ru-Ling Liao, Han Gao, Chih-Wei Hsu, Hongbin Liu, and Chun-Chi Chen,  “Motion vector coding and block merging in the versatile video coding standard,”  IEEE Transactions on Circuits and Systems for Video Technology, vol. 31, no. 10, pp. 3848–3861, 2021.
  * [8] Xin Jin, King Ngan, and Guangxi Zhu,  “Combined inter-intra prediction for high definition video coding,”  01 2007.
  * [9] Adam Wieckowski, Jackie Ma, Heiko Schwarz, Detlev Marpe, and Thomas Wiegand,  “Fast partitioning decision strategies for the upcoming versatile video coding (vvc) standard,”  in 2019 IEEE International Conference on Image Processing (ICIP). IEEE, 2019, pp. 4130–4134.