# Fast quantum integer multiplication with zero ancillas

Gregory D. Kahanamoku-Meyer  Corresponding author. Email<EMAIL_ADDRESS>Department of Physics, University of California at Berkeley, Berkeley, CA
94720  Research Laboratory of Electronics, Massachusetts Institute of
Technology, Cambridge, MA 02139  Norman Y. Yao  Department of Physics,
University of California at Berkeley, Berkeley, CA 94720  Department of
Physics, Harvard University, Cambridge, MA 02139

###### Abstract

The multiplication of superpositions of numbers is a core operation in many
quantum algorithms. The standard method for multiplication (both classical and
quantum) has a runtime quadratic in the size of the inputs. Quantum circuits
with asymptotically fewer gates have been developed, but generally exhibit
large overheads, especially in the number of ancilla qubits. In this work, we
introduce a new paradigm for sub-quadratic-time quantum multiplication with
zero ancilla qubits—the only qubits involved are the input and output
registers themselves. Our algorithm achieves an asymptotic gate count of
$\mathcal{O}(n^{1+\epsilon})$ for any $\epsilon>0$; with practical choices of
parameters, we expect scalings as low as $\mathcal{O}(n^{1.3})$. Used as a
subroutine in Shor’s algorithm, our technique immediately yields a factoring
circuit with $\mathcal{O}(n^{2+\epsilon})$ gates and only $2n+\mathcal{O}(\log
n)$ qubits; to our knowledge, this is by far the best qubit count of any
factoring circuit with a sub-cubic number of gates. Used in Regev’s recent
factoring algorithm, the gate count is $\mathcal{O}(n^{1.5+\epsilon})$.
Finally, we demonstrate that our algorithm has the potential to outperform
previous proposals at problem sizes relevant in practice, including yielding
the smallest circuits we know of for classically-verifiable quantum advantage.

##  I Introduction

Quantum circuits that perform arithmetic on superpositions of numbers are the
fundamental building block of many quantum algorithms, from factoring to
protocols for certifiable random number generation and efficiently-verifiable
quantum computational advantage [1, 2, 3, 4, 5, 6, 7, 8]. The feasibility of
realizing these applications depends on the efficiency with which this
arithmetic, and in particular multiplication, can be implemented. The standard
way of performing multiplication, in both the classical and quantum setting,
is via the “schoolbook” algorithm using $\mathcal{O}(n^{2})$ gates, where $n$
is the size of the input. Asymptotically sub-quadratic-time algorithms have
been known for over half a century, but have overheads that make them useful
only for the multiplication of very large values.111For a classical example,
the GNU multiple-precision arithmetic library uses a threshold of 2176 bit
inputs to switch away from the schoolbook method. In the quantum setting, the
reversibility constraint imposed by unitarity generally makes these overheads
even worse.

Despite this challenge, several works have explored the implementation of
quantum circuits for sub-quadratic-time multiplication, largely focusing on
the Karatsuba algorithm which exhibits an asymptotic run-time of roughly
$\mathcal{O}(n^{1.58})$ [9, 10, 11, 12, 13]. A significant obstacle to this
effort is the recursive structure of fast multiplication algorithms—making
them reversible involves storing intermediate data, which ultimately requires
utilizing a large number of ancilla qubits. Notable recent work has, for the
first time, reduced the number of ancillas required for quantum Karatsuba
multiplication to linear in the size of the inputs [11]. However, in practice,
the number of ancillas still dominates the total qubit cost. Moreover, it has
remained an open question whether sub-quadratic-time quantum multiplication
with fewer than $\mathcal{O}(n)$ ancillas is even possible.

Figure 1: Recursive structure of the quantum multiplication circuit, for $k=3$
classical-quantum multiplication. The green boxes represent phase rotations
proportional to the product of registers (the PhaseProduct operation described
in Section III), which are recursively decomposed; the unlabeled tall, narrow
unitaries represented by white boxes are in-place quantum addition circuits.
The qubits of the input and output registers are interleaved; the top pair of
qubits corresponds to the least significant bit of the input and output
registers, the next pair the second least significant bits of each, and so on.
At the last level of recursion, phase rotations on individual pairs of qubits
can be implemented directly as two-qubit controlled phase rotation gates
(represented here by green circles). Note that the locality of the gates
improves as the recursion proceeds.

In this work, we introduce a new paradigm for sub-quadratic-time
multiplication on a quantum computer, that requires zero ancilla qubits (Fig.
1): the only qubits needed are those used to store the input and output data.
In terms of gate counts, our algorithms achieve asymptotic scaling as low as
$\mathcal{O}(n^{1+\epsilon})$ for any $\epsilon>0$ (see Table 1). Our approach
combines ideas from classical fast multiplication with an inherently quantum
technique where arithmetic is performed in the phase of a quantum state. Using
our algorithm as a subroutine in Shor’s algorithm for factoring yields
circuits that require only $2n+\mathcal{O}(\log n)$ total qubits, where $n$ is
the length of the integer to be factored, while obtaining an asymptotic gate
count of $\mathcal{O}(n^{2+\epsilon})$ for arbitrarily small $\epsilon$; to
the best of our knowledge, this represents by far the fewest qubits required
to perform Shor’s algorithm in fewer than $\mathcal{O}(n^{3})$ gates.222See
Table 1 of [14] for a recent comparison of proposals for the implementation of
Shor’s algorithm. Moreover, used as a subroutine in a space-optimized version
of Regev’s recent factoring algorithm, our multiplication algorithm yields
quantum circuits for factoring with as few as $\mathcal{O}(n^{1.5+\epsilon})$
gates and $\mathcal{O}(n)$ qubits [8, 15].

One might naturally wonder whether our algorithm’s elimination of ancilla
qubits comes at the cost of significant overhead in constant factors. This is
not the case. Indeed, for multiplication of a quantum integer register by a
classical value (the key operation for implementing Shor’s algorithm), we find
it is possible to construct circuits for which the number of qubits required
is dramatically reduced, while the gate count remains small (see Table 2).
Intuitively, this gate efficiency comes from a fundamental feature of our
algorithm’s structure: that a substantial portion of the overhead implicit in
fast multiplication can be performed in classical precomputation. Highlighting
the versatility of our algorithm, we also investigate its application to an
efficiently-verifiable proof of quantumness [6] (see Table 3).

##  II Background and Framework

We focus on the implementation of two related unitaries, which are defined by
their action on product states (and can be extended by linearity to
superpositions of inputs):

$\displaystyle\mathcal{U}_{c\times q}(a)\ket{x}\ket{w}$
$\displaystyle=\ket{x}\ket{w+ax}$ (1) $\displaystyle\mathcal{U}_{q\times
q}\ket{x}\ket{y}\ket{w}$ $\displaystyle=\ket{x}\ket{y}\ket{w+xy}.$ (2)

The _classical-quantum multiply_ , $\mathcal{U}_{c\times q}(a)$, corresponds
to the multiplication of an integer stored in a quantum register, $\ket{x}$,
by a classical number, $a$. The result is added to an output register that
initially contains the value $\ket{w}$. Similarly, the _quantum-quantum
multiply_ , $\mathcal{U}_{q\times q}$, corresponds to the multiplication of a
quantum integer $\ket{x}$ by another quantum integer $\ket{y}$. The value is
again added to an output register $\ket{w}$.

In order to perform both of the above unitaries in sub-quadratic-time without
ancillas, we draw inspiration from two key ideas: (i) the classical fast (sub-
quadratic) multiplication algorithm of Toom and Cook and (ii) the quantum
Fourier transform based multiplication algorithm of Draper. At first glance,
these two seem entirely incompatible. Toom-Cook gains its speed from the reuse
of intermediate values that have already been computed; however, Draper
arithmetic is performed via quantum phase rotations, which, once applied,
cannot be reused. The crux of our result uses an inexpensive classical
precomputation to obviate the reuse of intermediate values; this
simultaneously removes the need to store intermediate values, while also
reducing the size of the quantum circuit.

###  II.1 Classical Toom-Cook multiplication

The most straightforward algorithm for multiplication is known as the
schoolbook algorithm and has a runtime of $\mathcal{O}(n^{2})$. It decomposes
the product as $xy=\sum_{i,j}b^{i+j}x_{i}y_{j}$ for a base $b$, where $x_{i}$
is the $i^{\mathrm{th}}$ base-$b$ digit of $x$ (and $y_{j}$ is defined
similarly). With $b=2$, the circuit is quite simple, because the binary
product $x_{i}y_{j}$ is an AND gate.

Starting in the 1960s, it was realized that multiplication could be performed
faster than $\mathcal{O}(n^{2})$ by breaking the inputs into pieces, and
carefully combining linear combinations of those pieces. To establish the
notation that we will use throughout this work, we begin by expressing the
classical Toom-Cook multiplication algorithm in the language of linear algebra
[16, 17].

Consider an $n$-bit value $x$ written in base $b=2^{n/k}$ for some integer
$k$, such that $x=\sum_{i=0}^{k-1}x_{i}b^{i}$. In this representation, each of
the $k$ digits, $x_{i}$, corresponds to $n/k$ bits of the binary
representation of $x$. The above sum can be expressed as the inner product of
two length-$k$ vectors: $\mathbf{x}=(x_{k-1},\cdots,x_{1},x_{0})$ and
$\mathbf{e}_{b}=(b^{k-1},\cdots,b,1)$,

$x=\sum_{i=0}^{k-1}x_{i}b^{i}=\mathbf{x}^{\intercal}\mathbf{e}_{b}.$ (3)

Allowing $b$ to vary, Eq. 3 defines a polynomial $x(b)$ whose coefficients are
the elements of $\mathbf{x}$. Letting $y(b)$ be the polynomial corresponding
to another integer $y$, one can define a third polynomial $p(b)=x(b)y(b)$ for
which $p(2^{n/k})$ equals the numerical value of the integer product $p=xy$.
Thus, the problem of integer multiplication can be expressed as polynomial
multiplication.

The key insight of the Toom-Cook algorithm is that the polynomial $p(b)$ is
uniquely determined by its value at $q=2k-1$ points $w_{\ell}$, where
$\ell\in\\{0,\cdots,q-1\\}$, and that each value, $p(w_{\ell})$, can be
computed as the pointwise multiplication $x(w_{\ell})y(w_{\ell})$. With an
appropriate choice of $w_{\ell}$, the values $x(w_{\ell})$ and $y(w_{\ell})$
are integers of roughly $n/k$ bits. Thus $p(b)$ can be computed via only
$q=2k-1$ multiplications (compare to the $k^{2}$ required by the schoolbook
algorithm) of size $\sim n/k$. By recursively applying this construction to
the smaller products of size $n/k$, and then the resulting products of size
$n/k^{2}$ (and so on), the Toom-Cook algorithm achieves an asymptotic runtime
of $\mathcal{O}(n^{\log_{k}q})$, which is sub-quadratic for all $k\geq 2$.

The evaluation of the polynomial $x(b)$ at the points $w_{\ell}$ can be
expressed via the matrix $A$ whose $\ell^{\mathrm{th}}$ row is
$\mathbf{e}_{w_{\ell}}$ (which is known as a Vandermonde matrix):

$A=\begin{pmatrix}w_{0}^{q-1}&\cdots&w_{0}^{2}&w_{0}&1\\\
w_{1}^{q-1}&\cdots&w_{1}^{2}&w_{1}&1\\\ \vdots&\ddots&\vdots&\vdots&\vdots\\\
w_{q-1}^{q-1}&\cdots&w_{q-1}^{2}&w_{q-1}&1\\\ \end{pmatrix}.$ (4)

Then, the vector $\tilde{\mathbf{x}}$ containing the values $x(w_{\ell})$ is
given by $\tilde{\mathbf{x}}=A\mathbf{x}$, and the vector
$\tilde{\mathbf{p}}$, containing the values $p(w_{\ell})$, can be expressed as
$\tilde{\mathbf{p}}=\tilde{\mathbf{x}}\circ\tilde{\mathbf{y}}=A\mathbf{x}\circ
A\mathbf{y}$, where $\circ$ denotes the pointwise product of vectors. (We
implicitly extend the vectors $\mathbf{x}$ and $\mathbf{y}$ to length $q$
elements by inserting zeros for all indices larger than $k-1$).

As long as the points $w_{\ell}$ are chosen such that $A$ is invertible, the
vector of polynomial coefficients $\mathbf{p}$ can be computed as
$A^{-1}\tilde{\mathbf{p}}$. This polynomial can then be evaluated at
$b=2^{n/k}$ by taking the inner product of $\mathbf{p}$ and
$\mathbf{e}_{2^{n/k}}$, yielding the integer product $p$. Putting it all
together, the integer value $p=xy$ can be computed as

$xy=\mathbf{e}_{2^{n/k}}^{\intercal}A^{-1}(A\mathbf{x}\circ A\mathbf{y}).$ (5)

For clarity of intuition, we provide an explicit instantiation of this
expression for $k=2$—corresponding to the Karatsuba algorithm—in Appendix A.

###  II.2 Quantum multiplication via the QFT

Many quantum circuits for multiplication mirror classical ones, with clever
optimizations to account for unitarity [18, 19, 20, 21, 22, 23, 11]. Here, we
describe an early result by Draper, which stands out as an example of a
quantum arithmetic circuit that has no classical analogue [24]. Consider the
quantum Fourier transform (QFT) of the output register in Eq. 1, before and
after the application of the classical-quantum multiply unitary
$\mathcal{U}_{c\times q}(a)$,

$\displaystyle\ket{x}\ket{w}$
$\displaystyle\xrightarrow{(\mathbb{I}\otimes\mathsf{QFT})}\sum_{z}\ket{x}\exp\left(2\pi
iwz/2^{2n}\right)\ket{z}$ $\displaystyle\ket{x}\ket{w+ax}$
$\displaystyle\xrightarrow{(\mathbb{I}\otimes\mathsf{QFT})}\sum_{z}\ket{x}\exp\left(2\pi
i(w+ax)z/2^{2n}\right)\ket{z}.$

The two only differ by a phase, $2\pi axz/2^{2n}$, on each element of the
superposition over $\ket{z}$. This suggests the following strategy for
implementing $\mathcal{U}_{c\times q}(a)$: (i) apply a QFT to the output
register, (ii) apply the diagonal unitary

$\tilde{\mathcal{U}}_{c\times q}(a)\ket{x}\ket{z}=\exp(2\pi
iaxz/2^{2n})\ket{x}\ket{z}$ (6)

and then (iii) apply an inverse QFT. In effect, this decomposes the original
unitary as: $\mathcal{U}_{c\times
q}(a)=(\mathbb{I}\otimes\mathsf{IQFT})\tilde{\mathcal{U}}_{c\times
q}(a)(\mathbb{I}\otimes\mathsf{QFT})$. For the quantum-quantum multiplication
unitary $\mathcal{U}_{q\times q}$ [Eq. 2], the strategy is identical, with the
diagonal unitary being,

$\tilde{\mathcal{U}}_{q\times q}\ket{x}\ket{y}\ket{z}=\exp(2\pi
ixyz/2^{2n})\ket{x}\ket{y}\ket{z}$ (7)

To date, the phase rotation unitaries, $\tilde{\mathcal{U}}_{c\times q}(a)$
and $\tilde{\mathcal{U}}_{q\times q}$, have been constructed via the binary
schoolbook decomposition of the products $xz$ and $xyz$:

$\displaystyle\exp(2\pi iaxz/2^{2n})$ $\displaystyle=\prod_{ik}\exp(2\pi
ia2^{i+k}x_{i}z_{k}/2^{2n})$ (8) $\displaystyle\exp(2\pi ixyz/2^{2n})$
$\displaystyle=\prod_{ijk}\exp(2\pi i2^{i+j+k}x_{i}y_{j}z_{k}/2^{2n}).$ (9)

Since the products $x_{i}z_{k}$ and $x_{i}y_{j}z_{k}$ are over binary values,
they are 0 if any of the bits are zero (corresponding to no phase shift) and 1
if all the bits are 1. Thus, the phases can be implemented as a series of
either singly-controlled or doubly-controlled phase rotations, for
$\tilde{\mathcal{U}}_{c\times q}(a)$ and $\tilde{\mathcal{U}}_{q\times q}$,
respectively. For the classical-quantum case, the number of controlled phase
rotations is $\mathcal{O}(n^{2})$, while for the quantum-quantum case, one has
to perform $\mathcal{O}(n^{3})$ phase rotations; the difference comes from the
fact that one can scale the entire phase rotation by $a$ since it is a
classical value.333In both cases, the quantum Fourier transforms (and their
inverses) can be performed in $\mathcal{O}(n^{2})$ time with no ancilla
qubits. Thus, they do not affect the asymptotic complexity. In both cases, the
number of ancilla qubits required is zero.

##  III Sub-quadratic quantum multiplication without ancillas

Our main result is that it is possible to implement the phases of Eqs. 6 and 7
via a construction based on the Toom-Cook algorithm, without ancilla qubits.
As a corollary, we also find an algorithm for the exact quantum Fourier
transform using a sub-quadratic number of gates and no ancillas. Taken
together, these results yield ancilla-free algorithms for fast quantum integer
multiplication, whose explicit asymptotic scalings for are shown in Table 1.

Algorithm | Gate count $\mathcal{O}(n^{\log_{k}(2k-1)})$
---|---
Eq. 8 | $\mathcal{O}(n^{2})$
$k=2$ | $\mathcal{O}(n^{1.58\cdots})$
$k=5$ | $\mathcal{O}(n^{1.37\cdots})$
$k=8$ | $\mathcal{O}(n^{1.30\cdots})$
(a) Classical-quantum multiplication

Algorithm | Gate count $\mathcal{O}(n^{\log_{k}(3k-1)})$
---|---
Eq. 9 | $\mathcal{O}(n^{3})$
$k=2$ | $\mathcal{O}(n^{2})$
$k=3$ | $\mathcal{O}(n^{1.77\cdots})$
$k=6$ | $\mathcal{O}(n^{1.55\cdots})$
$k=9$ | $\mathcal{O}(n^{1.46\cdots})$
(b) Quantum-quantum multiplication

Table 1:  Asymptotic scaling of gate counts for various selected $k$. Note
that the constant factors (which are not visible in big-$\mathcal{O}$
notation) become worse as $k$ increases, resulting in a tradeoff. The values
of $k$ shown here are practical for input integers up to a few thousand bits
(see Sec. V).

Phase rotation algorithms—Consider a generalization of the phase rotation of
Eq. 6, which we denote the PhaseProduct and define as follows:

$\mathsf{PhaseProduct}\left(\phi\right)\ket{x}\ket{z}=\exp(i\phi
xz)\ket{x}\ket{z}.$ (10)

Note that $\mathsf{PhaseProduct}(2\pi a/2^{2n})=\tilde{\mathcal{U}}_{c\times
q}(a)$. Our goal is to decompose the phase $\phi xz$ into a sum of many
phases, that are each easier to implement. At first glance, the Toom-Cook
construction does not seem helpful: it involves both recursion and a
complicated linear algebra expression (Eq. 5) that does not obviously
decompose into a single sum. Our key observation is that
$\mathbf{e}_{2^{n/k}}^{\intercal}$ and $A^{-1}$ are classically-known values,
and thus a vector $\vec{\phi}=\phi\mathbf{e}_{2^{n/k}}^{\intercal}A^{-1}$ can
be pre-computed. Using this in Eq. 5 for the product $xz$ yields

$\phi
xz=\sum_{\ell=0}^{q-1}\phi_{\ell}(A\mathbf{x})_{\ell}(A\mathbf{z})_{\ell},$
(11)

where $\phi_{\ell}$ are the elements of $\vec{\phi}$. To implement this
decomposition, the only linear algebra that needs to be done by the quantum
circuit is the computation of $(A\mathbf{x})_{\ell}$ and
$(A\mathbf{z})_{\ell}$. As in the classical Toom-Cook algorithm, by choosing
$w_{\ell}$ appropriately, these values can be computed using only a small
number of additions. Furthermore, each term of the sum in Eq. 11 has the form
of a PhaseProduct itself, but applied to values of size roughly $n/k$—thus,
the decomposition can be applied recursively (Fig. 1). Doing so yields a
decomposition of $\phi xz$ into $\mathcal{O}(n^{\log_{k}q})$ phase rotations
that can each be implemented in a constant number of gates, plus an
asymptotically negligible number of gates to compute (and uncompute) each
$(A\mathbf{x})_{\ell}$ and $(A\mathbf{z})_{\ell}$. Thus, the total number of
gates is $\mathcal{O}(n^{\log_{k}q})$ matching the asymptotic complexity of
the classical Toom-Cook algorithm.

Naively, ancilla registers would be needed to store the values
$(A\mathbf{x})_{\ell}$ and $(A\mathbf{z})_{\ell}$. We use the fact that
addition is reversible to instead overwrite parts of the input registers
$\ket{x}$ and $\ket{z}$ with those values—for example, to compute
$\ket{x_{0}+x_{1}}$ we may simply sum $x_{1}$ into the part of the $\ket{x}$
register already containing $x_{0}$.444The idea of using the reversibility of
addition to avoid allocating extra ancilla registers was proposed in [11] in
the context of the Karatsuba algorithm, however the overflow bits were not
considered in that work because $\mathcal{O}(n)$ ancilla qubits were already
otherwise being used. Because $(A\mathbf{x})_{\ell}$ and
$(A\mathbf{z})_{\ell}$ may take on slightly larger values than the registers
they overwrite, there will be some “overflow” carry bits. We avoid having to
store those in ancillas by devising a way to directly implement the part of
the PhaseProduct that involves the overflow bits, without ever explicitly
storing their value in ancillas (see Appendix C). Thus, no ancilla qubits are
needed for any part of the algorithm.

Next, we discuss quantum-quantum multiplication, which requires only a small
modification. We define an operation analogous to PhaseProduct, but which
operates on three registers instead of two:

$\displaystyle\mathsf{PhaseTripleProduct}\left(\phi\right)$
$\displaystyle\ket{x}\ket{y}\ket{z}$ $\displaystyle=\exp(i\phi xyz)$
$\displaystyle\ket{x}\ket{y}\ket{z}.$ (12)

In a less constrained setting, no special algorithm is required for the
product of three integers—one would simply take the product of two and then
multiply the third by the result. However, we cannot make use of intermediate
values like this in the phase; instead, the entire product must be decomposed
into a single sum. To this end, we introduce a modified version of the Toom-
Cook construction: we define three polynomials $x(b)$, $y(b)$, $z(b)$, and
compute their product $p(b)=x(b)y(b)z(b)$ by evaluating the polynomials at a
set of points $\\{w_{\ell}\\}$, finding their pointwise product, and then
interpolating. The only difference from the conventional Toom-Cook
construction is that $p(b)$ will now have degree $3(k-1)$ and is uniquely
determined by $q=3k-2$ points. Aside from the slightly larger value of $q$,
the structure of the quantum algorithm is the same. We may decompose $\phi
xyz$ as

$\phi
xyz=\sum_{\ell=0}^{q-1}\phi_{\ell}(A\mathbf{x})_{\ell}(A\mathbf{y})_{\ell}(A\mathbf{z})_{\ell}$
(13)

and recursively apply $\mathsf{PhaseTripleProduct}$ to each of the terms on
the right hand side. The larger value of $q$ leads to a slightly worse
asymptotic scaling, but still outperforms $\mathcal{O}(n^{2})$ for $k\geq 3$
[Table 1]. Again, overwriting parts of the input register to compute the
linear combinations, and performing the phase rotations corresponding to the
overflow bits directly, our algorithm requires no ancilla qubits.

Before moving on, we note two generalizations to $\mathcal{U}_{c\times q}(a)$
and $\mathcal{U}_{q\times q}$. Firstly, we observe that a classical-quantum-
quantum product $\mathcal{U}_{c\times q\times q}(a)$ may be implemented via
our quantum-quantum multiplication algorithm with no extra cost, by simply
scaling $\phi$ by $a$. Secondly, we note that the classical value $a$ does not
need to be an integer—it can be an arbitrary real number. Of course, for an
exact answer, the product must be representable by the output qubits, but even
if it is not, the final state will be a superposition of values near the
correct value (see Section 5.2.1 of [25]).

Fast exact quantum Fourier transform—Our algorithm for PhaseProduct can be
directly used to implement an exact quantum Fourier transform with the same
sub-quadratic gate count as classical-quantum multiplication, and again no
ancilla qubits. It is known that the exact quantum Fourier transform with
modulus $2^{n}$, $\mathsf{QFT}_{2^{n}}$, can be implemented via the following
three steps for any positive integer $m<n$ [26]:

  1. 1.

Apply $\mathsf{QFT}_{2^{m}}$ to the first $m$ qubits

  2. 2.

Apply $\mathsf{PhaseProduct}(2\pi/2^{n-m})$ to $\ket{x}\ket{y}$, where $x$ is
the value of the first $m$ qubits and $y$ is the value of the remaining $n-m$
qubits

  3. 3.

Apply $\mathsf{QFT}_{2^{n-m}}$ to the final $n-m$ qubits.

Setting $m=1$ corresponds to the standard construction for the QFT, using
$\mathcal{O}(n^{2})$ gates. By setting $m=n/2$ and using recursion to perform
steps 1 and 3, the asymptotic runtime can be improved. Prior to this work, the
proposed implementation for step 2 was to directly compute the product $xy$
into an ancilla register, via a classical fast multiplication algorithm that
has been compiled into a quantum circuit (at the cost of many ancillas).
Single-qubit phase rotations are then applied to the bits of the product,
after which it must be uncomputed. Instead, using our algorithm for
PhaseProduct in step 2 yields an algorithm for the exact QFT whose asymptotic
runtime matches that of our implementation of PhaseProduct without any
ancillas.

##  IV Space-time trade-offs

In the previous sections, we have shown that sub-quadratic quantum
multiplication can be performed with zero ancillas. However, we find that the
use of just a few ancillas is helpful for improving the practical performance
of our algorithm. As an example, explicitly storing the overflow bits
eliminates the extra work required to implement that portion of the phase
rotation (see Appendix C), and requires at most a few dozen extra qubits for
inputs of even several thousand bits. In what follows, we discuss three other
space-time trade-offs: (i) improving the base case of the recursion, (ii)
implementing arbitrary phase rotations with low overhead, and (iii) achieving
sub-linear circuit depth.

Base case optimization—One instance where a few extra ancillas can improve
performance considerably is in the base case of PhaseTripleProduct. In
particular, when the recursion has reached a sufficiently small size, denoted
$n_{\mathrm{base}}$, the phase rotation is performed directly. For
PhaseTripleProduct, implementing the base case via Eq. 9 requires
$n_{\mathrm{base}}^{3}$ doubly-controlled phase gates. Instead, we propose a
“semi-digital” implementation which requires only $\sim n_{\mathrm{base}}^{2}$
gates and is structured as follows. To implement the base-case phase rotation
of $\phi^{\prime}x^{\prime}y^{\prime}z^{\prime}$, start by explicitly
computing $\ket{x^{\prime}y^{\prime}}$ into an ancilla register using a
standard multiplication circuit. Next, use $2n_{\mathrm{base}}^{2}$ controlled
phase rotations (off of the ancilla register) to apply a phase of
$\phi(x^{\prime}y^{\prime})z^{\prime}$. Finally, uncompute the ancilla
register. This procedure requires at most $2n_{\mathrm{base}}$ ancilla
qubits.555We may even reduce the qubit count further by reusing a smaller
number of qubits to compute and uncompute small chunks of
$x^{\prime}y^{\prime}$ in sequence. We note that for PhaseProduct, the above
optimization is less important, since its base case can be implemented
directly with Eq. 8 using only $n_{\mathrm{base}}^{2}$ singly controlled phase
rotations.

Implementing $CR_{\phi}$ gates with low overhead via a phase gradient—In
settings where only a discrete gate set is native (such as in the logical
qubits of quantum error correcting codes), arbitrary phase rotations cannot be
directly implemented [27, 28, 29]. However, the overall cost of synthesizing
them can be dramatically reduced by using a few ancilla qubits to store a so-
called phase gradient state: $\ket{\Phi}=\sum_{\omega=0}^{2^{m}-1}e^{-2\pi
i\omega/2^{m}}\ket{\omega}$ [30, 31, 32, 33].666For a pedagogical exposition
of the use of phase gradients to implement the quantum Fourier transform, see
[31]. Here, $m=\lceil\log_{2}1/\eta\rceil$ where $\eta$ is the desired
precision and $\omega$ are integers. To apply an arbitrary phase $\phi$, one
can simply use a quantum addition circuit to increment the $\omega$ register
by $a=\lceil 2^{m}\phi/2\pi\rfloor$, where $\lceil\cdot\rfloor$ denotes
rounding to the nearest integer. This yields a phase shift of $\phi$, up to
precision $\eta=2^{-m}$:

$\sum_{\omega=0}^{2^{m}-1}e^{-2\pi i\omega/2^{m}}\ket{\omega+a}=e^{2\pi
ia/2^{m}}\ket{\Phi}\approx e^{i\phi}\ket{\Phi}.$ (14)

Thus, the $CR_{\phi}$ gates of our algorithm can be implemented via doubly-
controlled addition circuits. Because $\ket{\Phi}$ is an eigenstate of the
addition circuit, it is not destroyed by this process and can be reused an
arbitrary number of times.

As discussed earlier, the base case of both PhaseProduct and
PhaseTripleProduct can be implemented via $\sim n_{\mathrm{base}}^{2}$
$CR_{\phi}$ gates. A naive use of the phase gradient state to implement these
$\sim n_{\mathrm{base}}^{2}$ $CR_{\phi}$ gates yields a total cost of $\sim
n_{\mathrm{base}}^{2}\log(1/\eta)$ Toffoli gates. Instead, by again applying
the semi-digital optimization introduced above, one can convert the $\sim
n_{\mathrm{base}}^{2}$ $CR_{\phi}$ gates to $\sim n_{\mathrm{base}}^{2}$
Toffoli gates and $\sim n_{\mathrm{base}}$ $R_{\phi}$ gates. Putting
everything together, this technique implements the required rotations in
PhaseProduct and PhaseTripleProduct using a total of only
$\sim(n_{\mathrm{base}}^{2}+n_{\mathrm{base}}\log(1/\eta))$ Toffoli gates
(plus some Clifford gates). Since one expects
$\log_{2}(1/\eta)<n_{\mathrm{base}}$, this represents a considerable
improvement in the overall cost.

Circuit depth—The recursive tree structure of our algorithms lends itself well
to parallelization, as multiple branches of the tree can be performed
simultaneously (Fig. 1). Since each register is divided into $k$ parts,
ostensibly it is possible to perform $k$ recursive calls in parallel without
the need for ancillas. The challenge is that the $k$ linear combinations
$(A\mathbf{x})_{\ell}$ must be simultaneously computed in-place. We do not
know of any generic way to compute many specific linear combinations
reversibly and in-place like this; instead we have designed sequences of
operations by hand that are optimal or nearly optimal, which we record in
Appendix F. By parallelizing in this way (and using sub-linear-depth circuits
for the additions [21]), PhaseProduct can be implemented in depth
$\mathcal{O}(n^{\log_{k}2})$ and PhaseTripleProduct in
$\mathcal{O}(n^{\log_{k}3})$. These scalings are sublinear for $k>2$ and $k>3$
respectively. This sublinear-depth construction uses $\mathcal{O}(n)$ ancilla
qubits because the sublinear depth additions require $\mathcal{O}(n)$
ancillas. We also note that although the PhaseProduct and PhaseTripleProduct
can be implemented in sublinear depth, surprisingly, for the full
multiplication circuit we seem to be limited by the depth of the quantum
Fourier transforms; indeed, as far as we are aware, achieving a sub-linear
depth QFT requires a super-linear number of ancilla qubits.

##  V Applications

We now explore the use of our fast quantum multiplication algorithm as a
subroutine in two specific applications: Shor’s algorithm for factoring, and a
cryptographic proof of quantum computational advantage. While we leave the
explicit construction and optimization of quantum circuits to future work,
here, we perform careful estimates of the resources required and find
promising gate counts along with a dramatic reduction in the number of ancilla
qubits (see Tables 2 and 3).777The code used to estimate these gates counts is
available online: https://zenodo.org/doi/10.5281/zenodo.10871109 The
constructions we have presented so far implement multiplication over all
integers. Many applications, including both that we discuss in this section,
instead require multiplication over the integers modulo some $n$-bit integer
$N$.

Modular multiplication—It is straightforward to use generic multiplication to
implement modular multiplication [34], but doing so increases the size of the
circuit considerably and also requires the allocation of a $2n$-bit register
to store the full product.888In some cases, Zalka’s coset representation
allows approximate quantum modular multiplication to be performed with
essentially no overhead compared to standard multiplication [35]; fast
multipliers like the one presented in this paper and the Karatsuba circuits of
[36, 11] seem to have a structure that is fundamentally incompatible with
Zalka’s trick. Here, we propose to utilize an alternate strategy that
leverages the cyclic nature of the quantum phase to perform the modulo
operation automatically [6]. That is, we replace Eqs. 6 and 7 with

$\displaystyle\tilde{\mathcal{U}}_{c\times q}^{\prime}(a)\ket{x}\ket{z}$
$\displaystyle=\exp(2\pi iaxz/N)\ket{x}\ket{z}$ (15)
$\displaystyle\tilde{\mathcal{U}}_{q\times q}^{\prime}\ket{x}\ket{y}\ket{z}$
$\displaystyle=\exp(2\pi ixyz/N)\ket{x}\ket{y}\ket{z}$ (16)

such that multiples of $N$ in the product become multiples of $2\pi$ in the
phase. The only challenge with this strategy is that performing (inverse)
quantum Fourier transforms modulo arbitrary $N$ is too expensive. Instead, we
propose to still perform the QFT modulo $2^{n}$, which will yield a state
heavily weighted on the binary fraction that most closely approximates the
rational value $ax/N$ (classical-quantum multiplication) or $xy/N$ (quantum-
quantum multiplication); this binary fraction then uniquely identifies the
integer value. From quantum phase estimation, it is known that an output
register of $n+\mathcal{O}(\log(1/\eta))$ qubits is sufficient to ensure that
the final state is within $\eta$ of $ax/N$ or $xy/N$, respectively [25].

###  V.1 Classical-quantum multiplication: Shor’s algorithm

Perhaps the most obvious application of our result is to Shor’s algorithm for
integer factorization, which can be implemented via $\mathcal{O}(n)$
controlled, in-place classical-quantum modular multiplications. In particular,
it is possible to implement these multiplications to precision $\eta$ using a
total of $2n+\mathcal{O}(\log(1/\eta))$ qubits (see Appendix D for details).
Thus, asymptotically, our technique yields a circuit for factoring with
$\mathcal{O}(n^{2+\epsilon})$ gates (for arbitrarily small $\epsilon$) and
only $2n+\mathcal{O}(\log n)$ qubits. To our knowledge, this is by far the
best qubit count of any factoring circuit with a sub-cubic number of gates.

Algorithm | Asymptotic | Gate count (millions) | Ancillas
---|---|---|---
scaling | Toffoli | $CR_{\phi}$ |  $H$,$X$,CNOT
This work (standard QFT) | $\bm{\mathcal{O}(n^{1.29})}^{*}$ | 0.6 | 0.3 | 1.9 | 79
This work (phase gradient QFT) | $\bm{\mathcal{O}(n^{1.29})}^{*}$ | 0.9 | 0.1 | 3.2 | 80
Karatsuba [36] | $\mathcal{O}(n^{1.58})$ | 5.6 | — | 34 | 12730
Windowed [36] | $\mathcal{O}(n^{2}/\log^{2}n)$ | 1.8 | — | 2.5 | 4106
Schoolbook [36] | $\mathcal{O}(n^{2})$ | 6.4 | — | 38 | $1^{**}$
Table 2:  Circuit size estimates for one classical-quantum multiplication of
2048 bit numbers. Results from this work in bold, previous works in non-bold.
All estimates are in the “abstract circuit model” (no error correction or
routing costs included). Costs of previous works were computed using the Q#
code from [36]. Note that that code implements a regular multiplication rather
than a multiplication $\bmod N$; using Zalka’s coset representation of
integers [35] the cost of modular multiplication is expected to be roughly the
same as standard multiplication (except in the case of Karatsuba, for which it
does not seem possible to apply Zalka’s optimization and thus modular
multiplication may be considerably more expensive). Results for “this work”
are reported for modular multiplication, and include the quantum Fourier
transforms before and after the PhaseProduct, performed to a precision of
$\eta=10^{-12}$ per qubit. Note also that the “phase gradient QFT” estimates
here do not perform the arbitrary phase rotations in the base case of the
PhaseProduct via phase gradient state, only the QFTs. ∗We vary the parameter
$k$ throughout the recursion; this is the asymptotic scaling for $k=9$ which
is the value at the top level of recursion. ∗∗An ancilla count of 2048 was
reported in [36] and by the associated code; via qubit reuse it should be
possible to reduce this to 1. We did not explore how qubit reuse could be
applied to the other two constructions from [36].

In order to get a sense of circuit sizes in practice, we perform a detailed
analysis of the gate and qubit counts for the out-of-place multiplication of
2048-bit inputs (which is used to construct the in-place multiplications
described above). Specifically, we build the recursive tree (Fig. 1) and
carefully tabulate the number and type of gates in each layer. Via a careful
choice of the points $w_{\ell}$, we are able to devise efficient addition
sequences that minimize the overhead of forming linear combinations such as
$(A\mathbf{x})_{\ell}$. Moreover, at each layer, we explicitly optimize over
$k$, the number of pieces into which we divide the inputs (see e.g. Sec. II.1
and Table 1). Interestingly, we find that it is sometimes optimal to use
different values of $k$ even at the same level of recursion (see Appendix E).

Our results are depicted in Table 2. For context, we compare to recent work
which implements quantum circuits for classical-quantum multiplication,
including via the sub-quadratic Karatsuba algorithm [36].999Note that the
estimates from [36] are for non-modular arithmetic. As discussed in the
previous section, we have used a small number of extra ancillas to reduce the
constant factors on the gate count. Even with this optimization, the ancilla
counts are orders of magnitude less than comparable implementations! Moreover,
the gate counts are very promising, although impossible to compare directly to
alternate strategies without compiling to a common native gate set. We note
that the majority of the $CR_{\phi}$ gates occur not in the PhaseProduct but
in the QFTs performed on either side of it. As shown in Table 2, implementing
these QFTs using a phase gradient state reduces the total number of
$CR_{\phi}$ gates considerably.

###  V.2 $x^{2}\bmod N$: cryptographic proofs of quantum computational
advantage

Recently, much excitement has centered on the experimental realization of
quantum computational advantage—specifically, the performance of random
sampling tasks that are infeasible for even the world’s fastest classical
supercomputers [37, 38, 39, 40, 41]. So far, such experiments have had the
subtle characteristic that checking the results classically is as hard as, if
not harder than, generating the solution. However, recent theoretical work has
explored the use of cryptography to design “test of quantum computational
advantage” protocols which are classically efficiently verifiable, yet remain
hard to spoof [4, 5, 6]. When instantiated with quantum-secure cryptographic
assumptions, this class of protocols is also potentially useful for other
cryptographic tasks such as certifiable random number generation, remote state
preparation, and classical delegation of quantum computations to untrusted
devices [4, 3, 2, 7, 42]. Because much of the underlying cryptography involves
arithmetic, our algorithm for multiplication is directly relevant to the
protocols’ implementation.

Here, we focus specifically on a recent proposal for a test of quantum
computational advantage in which the quantum computer must evaluate the
function $f(x)=x^{2}\bmod N$ on a superposition of inputs $x$ [6]. This
function, and the protocol it is used in, permit a number of helpful
optimizations. For our multiplication algorithm, computing a square instead of
a general quantum-quantum multiplication is less costly because the linear
combinations like $(A\mathbf{x})_{\ell}$ need only be computed on two
registers ($x$ and $z$) instead of three. The protocol also has a built-in
measurement based uncomputation scheme, which allows ancilla qubits to be
discarded by measuring them in the Hadamard basis and recording the
measurement results, avoiding the need for explicit uncomputation of garbage
bits [6].

Algorithm | Asymptotic | Gate count (millions) | Total
---|---|---|---
scaling | Toffoli | $CR_{\phi}$ | Measmt. |  $H$,$X$,CNOT | qubits
This work, “fast” | $\bm{\mathcal{O}(n^{1.49})}^{*}$ | 0.6 | 0.7 | 0.4 | 0.8 | 2937
This work, “balanced” | $\bm{\mathcal{O}(n^{1.49})}^{*}$ | 0.7 | 0.7 | 0.4 | 1.1 | 2140
This work, “narrow” | $\bm{\mathcal{O}(n^{1.55})}^{*}$ | 1.9 | 1.4 | 1.0 | 2.7 | 1583
Prev. Fourier 1 [6] | $\mathcal{O}(n^{3})$ | — | $539^{**}$ | — | — | 1025
Prev. Fourier 2 [6] | $\mathcal{O}(n^{2}\log n)$ | — | 35 | — | — | 2062
“Digital” Karatsuba [6] | $\mathcal{O}(n^{1.58})$ | 1.6 | — | 1.1 | 1.6 | 6801
“Digital” Schoolbook [6] | $\mathcal{O}(n^{2})$ | 3.5 | — | 2.2 | 2.9 | 4097
Table 3:  Circuit size estimates for $x^{2}\bmod N$ in the context of the
proof of quantumness protocol, for 1024-bit $N$. Results from this work in
bold, previous results in non-bold. All estimates are in the “abstract circuit
model” (no error correction or routing costs included). ∗We vary the parameter
$k$ throughout the recursion; these are the asymptotic scalings for the values
of $k$ used at the top level of the recursions ($k=6$ for “narrow” and $k=8$
for “balanced” and “fast”). ∗∗These gates are doubly, rather than singly,
controlled $R_{\phi}$ gates.

We perform analogous estimates of the gate and qubit counts for computing
$f(x)$ with 1024-bit $x$ and $N$.101010This is a problem size we expect to be
infeasible for modern supercomputers: the protocol’s hardness is based on the
hardness of factoring $N$, and the largest publicly-known factorization of an
integer without a special form is of length 829 bits [43]. We explore three
constructions (“fast”, “balanced”, and “narrow”), with varying trade-offs
between qubit and gate count, and our results are shown in Table 3. The “fast”
version of our circuits aggressively uses the measurement-based uncomputation
scheme to reduce gate counts at the expense of a moderate number of ancillas.
The “balanced” version computes and uncomputes the sums in-place, requiring a
few more gates but many fewer ancillas. Finally, the “narrow” version never
stores the entire output and instead reuses a register of $n/2$ qubits twice,
further reducing qubit counts at the expense of a few more gates.

We observe that all three constructions reduce the total qubit counts
substantially when compared to previous “digital” implementations of the
Karatsuba and Schoolbook algorithms. When compared with two other algorithms
optimized for qubit count (at the expense of gate count), our algorithms
achieve similarly low qubit counts while reducing gate counts by orders of
magnitude.

##  VI Outlook

The optimization of quantum circuits for performing arithmetic on
superpositions of inputs has been the subject of study for decades. In this
work, we have introduced a technique for performing quantum multiplication in
a sub-quadratic number of gates with zero ancillas.

Our results open the door to a number of intriguing open questions. First, can
factoring be performed in under $\mathcal{O}(n^{3})$ gates with even fewer
than the $2n+\mathcal{O}(\log n)$ qubits reported here? One obvious path
towards achieving this would be to devise a new way to implement modular
arithmetic in our algorithm without requiring $\mathcal{O}(\log n)$ extra
qubits [35]. Another path is to apply our multiplication algorithm to a space-
efficient version of Regev’s recent fast factoring algorithm [8, 15]; however,
it seems that further space optimization would be required. Second, we note
that there are likely more clever constructions for implementing the
controlled arbitrary phase rotations of the base case, which avoid the need to
compute and uncompute an intermediate ancilla register. Finally, it is worth
exploring whether our strategy is compatible with even faster classical
multiplication algorithms, such as the $\mathcal{O}(n\log n\log\log n)$
Schonhage-Strassen algorithm [44, 17].

The authors would like acknowledge the insights of and discussions with Tanuj
Khattar, Craig Gidney, Isaac Chuang, Seyoon Ragavan, Katherine van Kirk, and
John Blue. This work was supported by the NSF QLCI Award OMA-2016245 and the
NSF STAQ II program. N.Y.Y. acknowledges support from a Simons Investigator
Award.

## References

  * [1] P. W. Shor.  “Polynomial-Time Algorithms for Prime Factorization and Discrete Logarithms on a Quantum Computer.”  _SIAM Journal on Computing_ , 26(5), 1484–1509, October 1997.  ISSN 0097-5397.  doi:10.1137/S0097539795293172.
  * [2] U. Mahadev.  “Classical Verification of Quantum Computations.”  In “2018 IEEE 59th Annual Symposium on Foundations of Computer Science (FOCS),” pp. 259–267. October 2018\.  ISSN 2575-8454.  doi:10.1109/FOCS.2018.00033.
  * [3] A. Gheorghiu and T. Vidick.  “Computationally-Secure and Composable Remote State Preparation.”  In “2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS),” pp. 1024–1033. November 2019\.  ISSN 2575-8454.  doi:10.1109/FOCS.2019.00066.
  * [4] Z. Brakerski, P. Christiano, U. Mahadev, et al.  “A Cryptographic Test of Quantumness and Certifiable Randomness from a Single Quantum Device.”  _Journal of the ACM (JACM)_ , August 2021.  doi:10.1145/3441309.
  * [5] Z. Brakerski, V. Koppula, U. Vazirani, and T. Vidick.  “Simpler Proofs of Quantumness.”  In S. T. Flammia (Editor), “15th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2020),” volume 158 of _Leibniz International Proceedings in Informatics (LIPIcs)_ , pp. 8:1–8:14. Schloss Dagstuhl–Leibniz-Zentrum für Informatik, Dagstuhl, Germany, 2020.  ISBN 978-3-95977-146-7.  ISSN 1868-8969.  doi:10.4230/LIPIcs.TQC.2020.8.
  * [6] G. D. Kahanamoku-Meyer, S. Choi, U. V. Vazirani, and N. Y. Yao.  “Classically-Verifiable Quantum Advantage from a Computational Bell Test.”  _Nature Physics_ , 18(8), 918–924, August 2022.  ISSN 1745-2473, 1745-2481.  doi:10.1038/s41567-022-01643-7.
  * [7] Z. Brakerski, A. Gheorghiu, G. D. Kahanamoku-Meyer, et al.  “Simple Tests of Quantumness Also Certify Qubits.”  In H. Handschuh and A. Lysyanskaya (Editors), “Advances in Cryptology – CRYPTO 2023,” Lecture Notes in Computer Science, pp. 162–191. Springer Nature Switzerland, Cham, 2023.  ISBN 978-3-031-38554-4.  doi:10.1007/978-3-031-38554-4˙6.
  * [8] O. Regev.  “An Efficient Quantum Factoring Algorithm.”, August 2023.  doi:10.48550/arXiv.2308.06572.
  * [9] A. Parent, M. Roetteler, and M. Mosca.  “Improved reversible and quantum circuits for Karatsuba-based integer multiplication.”  In M. M. Wilde (Editor), “12th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2017),” volume 73 of _Leibniz International Proceedings in Informatics (LIPIcs)_ , pp. 7:1–7:15. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl, Germany, 2018.  ISBN 978-3-95977-034-7.  ISSN 1868-8969.  doi:10.4230/LIPIcs.TQC.2017.7.
  * [10] S. Dutta, D. Bhattacharjee, and A. Chattopadhyay.  “Quantum circuits for Toom-Cook multiplication.”  _Physical Review A_ , 98(1), 012311, July 2018.  doi:10.1103/PhysRevA.98.012311.
  * [11] C. Gidney.  “Asymptotically Efficient Quantum Karatsuba Multiplication.”  _arXiv:1904.07356 [quant-ph]_ , April 2019.
  * [12] H. T. Larasati, A. M. Awaludin, J. Ji, and H. Kim.  “Quantum Circuit Design of Toom 3-Way Multiplication.”  _Applied Sciences_ , 11(9), 3752, January 2021.  ISSN 2076-3417.  doi:10.3390/app11093752.
  * [13] L. Kowada, R. Portugal, and C. Figueiredo.  “Reversible Karatsuba’s Algorithm.”  _J. UCS_ , 12, 499–511, January 2006.
  * [14] C. Gidney and M. Ekerå.  “How to factor 2048 bit RSA integers in 8 hours using 20 million noisy qubits.”  _Quantum_ , 5, 433, April 2021.  doi:10.22331/q-2021-04-15-433.
  * [15] S. Ragavan and V. Vaikuntanathan.  “Space-Efficient and Noise-Robust Quantum Factoring.”, February 2024.  doi:10.48550/arXiv.2310.00899.
  * [16] M. Bodrato.  “Towards Optimal Toom-Cook Multiplication for Univariate and Multivariate Polynomials in Characteristic 2 and 0.”  In C. Carlet and B. Sunar (Editors), “Arithmetic of Finite Fields,” Lecture Notes in Computer Science, pp. 116–133. Springer, Berlin, Heidelberg, 2007.  ISBN 978-3-540-73074-3.  doi:10.1007/978-3-540-73074-3˙10.
  * [17] D. E. Knuth.  _The Art of Computer Programming, Volume II: Seminumerical Algorithms, 3rd Edition_.  Addison-Wesley, 1998.  ISBN 0-201-89684-2.
  * [18] V. Vedral, A. Barenco, and A. Ekert.  “Quantum networks for elementary arithmetic operations.”  _Physical Review A_ , 54(1), 147–153, July 1996.  doi:10.1103/PhysRevA.54.147.
  * [19] C. Zalka.  “Fast versions of Shor’s quantum factoring algorithm.”  _arXiv:quant-ph/9806084_ , June 1998.
  * [20] S. A. Cuccaro, T. G. Draper, S. A. Kutin, and D. P. Moulton.  “A new quantum ripple-carry addition circuit.”  _arXiv:quant-ph/0410184_ , October 2004.
  * [21] T. G. Draper, S. A. Kutin, E. M. Rains, and K. M. Svore.  “A logarithmic-depth quantum carry-lookahead adder.”  _Quantum Information & Computation_, 6(4), 351–369, July 2006.  ISSN 1533-7146.
  * [22] T. Häner, M. Roetteler, and K. M. Svore.  “Factoring using 2n + 2 qubits with Toffoli based modular multiplication.”  _Quantum Information & Computation_, 17(7-8), 673–684, June 2017\.  ISSN 1533-7146.
  * [23] C. Gidney.  “Factoring with n+2 clean qubits and n-1 dirty qubits.”  _arXiv:1706.07884 [quant-ph]_ , January 2018.
  * [24] T. G. Draper.  “Addition on a Quantum Computer.”  _arXiv:quant-ph/0008033_ , August 2000.
  * [25] M. A. Nielsen and I. L. Chuang.  _Quantum Computation and Quantum Information: 10th Anniversary Edition_.  Cambridge University Press, USA, 10th anniversary ed edition, 2011.  ISBN 978-1-107-00217-3.
  * [26] R. Cleve and J. Watrous.  “Fast parallel circuits for the quantum Fourier transform.”  In “Proceedings 41st Annual Symposium on Foundations of Computer Science,” pp. 526–536. November 2000.  ISSN 0272-5428.  doi:10.1109/SFCS.2000.892140.
  * [27] R. W. Heeres, P. Reinhold, N. Ofek, et al.  “Implementing a universal gate set on a logical qubit encoded in an oscillator.”  _Nature Communications_ , 8(1), 94, July 2017.  ISSN 2041-1723.  doi:10.1038/s41467-017-00045-1.
  * [28] R. Acharya, I. Aleiner, R. Allen, et al.  “Suppressing quantum errors by scaling a surface code logical qubit.”  _Nature_ , 614(7949), 676–681, February 2023.  ISSN 1476-4687.  doi:10.1038/s41586-022-05434-1.
  * [29] D. Bluvstein, H. Levine, G. Semeghini, et al.  “A quantum processor based on coherent transport of entangled atom arrays.”  _Nature_ , 604(7906), 451–456, April 2022.  ISSN 1476-4687.  doi:10.1038/s41586-022-04592-6.
  * [30] A. Kitaev, A. Shen, and M. Vyalyi.  _Classical and Quantum Computation_ , volume 47 of _Graduate Studies in Mathematics_.  American Mathematical Society, May 2002.  ISBN 978-0-8218-3229-5 978-1-4704-0927-2 978-1-4704-2011-6 978-1-4704-1800-7.  doi:10.1090/gsm/047.
  * [31] C. Gidney.  “Turning Gradients into Additions into QFTs.”  https://algassert.com/post/1620, July 2016.
  * [32] C. Gidney.  “Halving the cost of quantum addition.”  _Quantum_ , 2, 74, June 2018.  doi:10.22331/q-2018-06-18-74.
  * [33] Y. Nam, Y. Su, and D. Maslov.  “Approximate quantum Fourier transform with O ( n log( n )) T gates.”  _npj Quantum Information_ , 6(1), 1–6, March 2020.  ISSN 2056-6387.  doi:10.1038/s41534-020-0257-5.
  * [34] P. L. Montgomery.  “Modular multiplication without trial division.”  _Mathematics of Computation_ , 44(170), 519–521, 1985.  ISSN 0025-5718, 1088-6842.  doi:10.1090/S0025-5718-1985-0777282-X.
  * [35] C. Zalka.  “Shor’s algorithm with fewer (pure) qubits.”  _arXiv:quant-ph/0601097_ , January 2006.
  * [36] C. Gidney.  “Windowed quantum arithmetic.”, May 2019.  doi:10.48550/arXiv.1905.07682.
  * [37] F. Arute, K. Arya, R. Babbush, et al.  “Quantum supremacy using a programmable superconducting processor.”  _Nature_ , 574(7779), 505–510, October 2019.  ISSN 1476-4687.  doi:10.1038/s41586-019-1666-5.
  * [38] H.-S. Zhong, H. Wang, Y.-H. Deng, et al.  “Quantum computational advantage using photons.”  _Science_ , 370(6523), 1460–1463, December 2020.  ISSN 0036-8075, 1095-9203.  doi:10.1126/science.abe8770.
  * [39] Y. Wu, W.-S. Bao, S. Cao, et al.  “Strong Quantum Computational Advantage Using a Superconducting Quantum Processor.”  _Physical Review Letters_ , 127(18), 180501, October 2021.  doi:10.1103/PhysRevLett.127.180501.
  * [40] Q. Zhu, S. Cao, F. Chen, et al.  “Quantum computational advantage via 60-qubit 24-cycle random circuit sampling.”  _Science Bulletin_ , 67(3), 240–245, February 2022.  ISSN 2095-9273.  doi:10.1016/j.scib.2021.10.017.
  * [41] A. Morvan, B. Villalonga, X. Mi, et al.  “Phase transition in Random Circuit Sampling.”, December 2023\.  doi:10.48550/arXiv.2304.11119.
  * [42] A. Natarajan and T. Zhang.  “Bounding the Quantum Value of Compiled Nonlocal Games: From CHSH to BQP Verification.”  In “2023 IEEE 64th Annual Symposium on Foundations of Computer Science (FOCS),” pp. 1342–1348. November 2023\.  ISSN 2575-8454.  doi:10.1109/FOCS57990.2023.00081.
  * [43] F. Boudot, P. Gaudry, A. Guillevic, et al.  “The State of the Art in Integer Factoring and Breaking Public-Key Cryptography.”  _IEEE Security & Privacy_, 20(2), 80–86, March 2022.  ISSN 1558-4046.  doi:10.1109/MSEC.2022.3141918.
  * [44] A. Schönhage and V. Strassen.  “Schnelle Multiplikation großer Zahlen.”  _Computing_ , 7(3), 281–292, September 1971.  ISSN 1436-5057.  doi:10.1007/BF02242355.
  * [45] A. Barenco, C. H. Bennett, R. Cleve, et al.  “Elementary gates for quantum computation.”  _Physical Review A_ , 52(5), 3457–3467, November 1995.  doi:10.1103/PhysRevA.52.3457.
  * [46] S. Beauregard.  “Circuit for Shor’s algorithm using 2n+3 qubits.”  _Quantum Information & Computation_, 3(2), 175–185, March 2003\.  ISSN 1533-7146.

##  Appendix A The Karatsuba algorithm in the notation of Section II.1

Here for pedagogical purposes we briefly describe the Karatsuba algorithm, of
which Toom-Cook is a generalization, and show how it can be expressed in the
linear algebra notation introduced in Section II.1.

We begin by introducing the Karatsuba algorithm. It is an instance of Toom-
Cook with $k=2$, meaning that the input values are split into two pieces,
corresponding to the low and high halves of the bits: $x=2^{n/2}x_{1}+x_{0}$.
A product $xy$ of values expressed in this way can be expanded to

$xy=2^{n}x_{1}y_{1}+2^{n/2}(x_{1}y_{0}+x_{0}y_{1})+x_{0}y_{0}$ (17)

which replaces the product of $n$-bit values with four products of $n/2$-bit
values. The key to the Karatsuba algorithm is to observe that the quantity in
parentheses can be written as

$x_{1}y_{0}+x_{0}y_{1}=(x_{0}+x_{1})(y_{0}+y_{1})-x_{1}y_{1}-x_{0}y_{0}$ (18)

and since the second two products there can be reused from the other terms of
Eq. 17, the entire product can be computed using only three products of size
$n/2$ instead of four (plus a few additions). By recursively applying this
technique, one may compute the product $xy$ in only
$\mathcal{O}(n^{\log_{2}3})=\mathcal{O}(n^{1.58\cdots})$ operations.

In the notation of Section II.1, Karatsuba corresponds to setting $k=2$ and
using the evaluation points $w_{\ell}\in\\{0,1,\infty\\}$. For the sake of
gaining intuition, it may be helpful for the reader to work through Eq. 5
using these values—the result should be the Karatsuba decomposition of the
product $xy$:

$2^{n}x_{1}y_{1}+2^{n/2}\left[(x_{0}+x_{1})(y_{0}+y_{1})-x_{1}y_{1}-x_{0}y_{0}\right]+x_{0}y_{0}$
(19)

Note that arriving at the correct expression is easiest if one considers the
point $\infty$ as the unit fraction $1/0$; see Appendix B.

##  Appendix B Use of fractional $w_{\ell}$

In both classical Toom-Cook and in our algorithms, it is common to take some
of the evaluation points $w_{\ell}$ to be unit fractions $1/c$ for some $c$
(frequently one such point is $\infty$, which we may consider as a unit
fraction with $c=0$). In this section we describe how careful maneuvering
allows us to maintain a Vandermode matrix $A$ of integer values, even when
using such unit fractional $w_{\ell}$.

Consider a length-$q$ vector $\mathbf{e}_{1/c}$ (which is row $\ell$ of the
matrix $A$ for some $\ell$):

$\mathbf{e}_{1/c}=(1/c^{q-1},1/c^{q-2},\cdots,1/c^{2},1/c,1)$ (20)

Then the $\ell^{\mathrm{th}}$ element of $\tilde{\mathbf{x}}=A\mathbf{x}$ is
$\mathbf{e}_{1/c}\cdot\mathbf{x}$, where the length-$k$ vector $\mathbf{x}$
has been padded to length $q$ by inserting zero elements on the left.

The Toom-Cook algorithm is built on the fact that the $\ell^{\mathrm{th}}$
element of the product vector $\tilde{\mathbf{p}}$ is the pointwise product of
the $\ell^{\mathrm{th}}$ elements of $\tilde{\mathbf{x}}$ and
$\tilde{\mathbf{y}}$, that is,

$\mathbf{e}_{1/c}\cdot\mathbf{p}=(\mathbf{e}_{1/c}\cdot\mathbf{x})(\mathbf{e}_{1/c}\cdot\mathbf{y})$
(21)

Now, the key is to scale both sides of the expression by $c^{2k-2}$ and then
distribute that constant across the $\mathbf{e}$ vectors:

$c^{2k-2}\mathbf{e}_{1/c}\cdot\mathbf{p}=(c^{k-1}\mathbf{e}_{1/c}\cdot\mathbf{x})(c^{k-1}\mathbf{e}_{1/c}\cdot\mathbf{y})$
(22)

Now we define $\mathbf{e}^{\prime\prime}_{1/c}=c^{2k-2}\mathbf{e}_{1/c}$ and
(noting that $q=2k-1$) see that is has the form

$\mathbf{e}^{\prime\prime}_{1/c}=(1,c,c^{2},\cdots,c^{q-2},c^{q-1})$ (23)

all of which are integers. Similarly (dropping the elements that are
multiplied by padded zeros of $\mathbf{x}$) we have

$\mathbf{e}^{\prime}_{1/c}=(1,c,c^{2},\cdots,c^{k-2},c^{k-1})$ (24)

In summary, we may perform Toom-Cook using slightly adjusted matrices
$A^{\prime}$ and $A^{\prime\prime}$ which have only integer elements, as
follows. Let the first $m$ evaluation points $w_{\ell}$ be integers, and the
remaining $q-m$ points be unit fractions $1/c_{\ell}$. Then we construct the
$q\times k$ matrix

$A^{\prime}=\begin{pmatrix}w_{0}^{k-1}&w_{0}^{k-2}&\cdots&w_{0}&1\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
w_{m-1}^{k-1}&w_{m-1}^{k-2}&\cdots&w_{m-1}&1\\\
1&c_{m}&\cdots&c_{m}^{k-2}&c_{m}^{k-1}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&c_{q-1}&\cdots&c_{q-1}^{k-2}&c_{q-1}^{k-1}\\\ \end{pmatrix}$ (25)

and the $q\times q$ matrix

$A^{\prime\prime}=\begin{pmatrix}w_{0}^{q-1}&w_{0}^{q-2}&\cdots&w_{0}&1\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
w_{m-1}^{q-1}&w_{m-1}^{q-2}&\cdots&w_{m-1}&1\\\
1&c_{m}&\cdots&c_{m}^{q-2}&c_{m}^{q-1}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&c_{q-1}&\cdots&c_{q-1}^{q-2}&c_{q-1}^{q-1}\\\ \end{pmatrix}.$ (26)

Finally, we use them in a slightly modified version of Eq. 5:

$xy=\mathbf{e}_{2^{n/k}}^{\intercal}(A^{\prime\prime})^{-1}(A^{\prime}\mathbf{x}\circ
A^{\prime}\mathbf{y}).$ (27)

Now for a couple of remarks. In the previous paragraphs we have discussed only
the standard Toom-Cook decomposition, but this strategy applies equally well
to our modified Toom-Cook decomposition of $xyz$. In that case, the
appropriate rows of $A^{\prime}$ get scaled by $c_{\ell}^{k-1}$, and because
there are three factors, this matches a factor of $c^{q-1}=c^{3(k-1)}$ in the
corresponding rows of the $A^{\prime\prime}$ matrix.

Finally, we note that while $A^{\prime\prime}$ has integer entries, that does
not imply that $(A^{\prime\prime})^{-1}$ will. Indeed, it has been shown that
no choice of $w_{\ell}$ can avoid the need for division if $k>2$ [16]. For our
algorithm this is not a problem at all—the entries of
$(A^{\prime\prime})^{-1}$ get wrapped up into the $\phi_{\ell}$ in classical
precomputation, and thus the fact that they are not always integers does not
affect the quantum circuit. This is actually a big challenge for fast
multiplication circuits other than ours, however, because it is necessary to
perform division during the interpolation step, which can be expensive. This
issue was encountered by a recent work which explored a quantum implementation
of $k=3$ Toom-Cook, and found that the division operation led to a large
Toffoli count [12].

##  Appendix C Overflow bit phase application

As described in the main text, we avoid allocating extra registers to store
the linear combinations
$(A\mathbf{x})_{\ell}=\mathbf{x}^{\intercal}\mathbf{e}_{w_{\ell}}$ by
performing addition in-place, reusing existing qubits to store the summed
values.111111See footnote 4 However, because the value $(A\mathbf{x})_{\ell}$
may be a few bits larger than the portion of the input it is overwriting,
without additional tricks this technique would yield a multiplication
algorithm requiring $\mathcal{O}(1)$ ancillas at each level of the recursive
tree. Here we describe a construction by which these extra ancillas can be
avoided.

Figure 2: A circuit for implementing a phase rotation by
$\phi(x_{0}+x_{1})(z_{0}+z_{1})$, using only one dirty ancilla that is
borrowed from another part of the computation. As discussed in the main text,
this dirty ancilla can be borrowed from the inputs themselves, yielding an
implementation that requires zero ancilla qubits, dirty or clean. Here, the
unitaries with rounded corners represent phase rotations of the labeled factor
times the value of the inputs on which the unitary is being applied. The
operations labeled “Cuccaro 1/2” and “Cuccaro 2/2” represent the first and
second halves of a Cuccaro adder, respectively. The unitaries controlled off
of the top bit of a many-qubit register are controlled off of the qubit
temporarily holding the value of the carry bit of the Cuccaro adder in
progress. The operation “Cuccaro†” is the Cuccaro adder performed backwards,
uncomputing the sums.

###  C.1 Simple case

For pedagogical purposes, let us first consider a simple example: the
implementation of a PhaseProduct corresponding to a phase rotation of
$\phi(x_{0}+x_{1})(z_{0}+z_{1})$, where $\phi$ is a classically-known phase
factor, and we have as input the four registers
$\ket{x_{0}}\ket{x_{1}}\ket{z_{0}}\ket{z_{1}}$. For simplicity in this first
example, we will assume that all four registers are of the same length, which
we can denote as $m$ bits. The straightforward method for implementing this
phase rotation is to use a quantum adder to map
$\ket{x_{0}}\ket{x_{1}}\to\ket{x_{0}}\ket{x_{0}+x_{1}}$ (and similarly on the
$z$ registers), and then use a recursive call to our PhaseProduct algorithm on
the registers now holding $\ket{x_{0}+x_{1}}$ and $\ket{z_{0}+z_{1}}$. The
issue is that the value $x_{0}+x_{1}$ (resp. $z$) takes one more bit to
represent than $x_{1}$, so naively we need to use an ancilla qubit to hold the
“carry bit” produced at the top of the adder. This can be avoided as follows.

Observe that we can split up the value $x_{0}+x_{1}$ into the carry bit which
we denote $x_{c}$ and the remaining $m$ bits of the sum, such that
$x_{0}+x_{1}=2^{m}x_{c}+(x_{0}+x_{1}\bmod 2^{m})$ (and same for the $z$
values). Then with some rearrangement the product can be written as follows
(the reason for the precise form of which will become clear promptly):

$\phi(x_{0}+x_{1})(z_{0}+z_{1})=\phi 2^{m}x_{c}z_{0}+\phi 2^{m}x_{c}z_{1}+\phi
2^{m}z_{c}(x_{0}+x_{1}\bmod 2^{m})+\phi(x_{0}+x_{1}\bmod
2^{m})(z_{0}+z_{1}\bmod 2^{m})$ (28)

Thus we can implement our desired phase rotation by implementing separate
rotations for each of the terms on the right hand side of Eq. 28.

Cuccaro’s quantum adder [20] has a variant that takes two length $m$ quantum
integers, plus an “incoming carry bit” (of the same value as the least
significant bits of the input), and computes the sum of the integers plus the
incoming carry bit, writing the result into one of the input registers. It
only requires one extra qubit, which is the qubit in which the outgoing carry
(e.g. $x_{c}$) is ultimately stored. Conveniently, the value of that outgoing
carry is set simply with a CNOT from one of the other qubits, which
temporarily holds the value $x_{c}$. This means that if we pause the Cuccaro
adder mid-execution, we can forego the extra qubit needed to store the carry,
as well as the CNOT setting it, instead using the value directly. Crucially,
we have designed the first three terms of Eq. 28 such that they consist of the
product of a single bit value by an $m$ bit value—and thus their phase
rotation can be implemented directly via $m$ $CR_{\phi}$ gates, which is
asymptotically negligible compared to the superlinear cost of the recursive
call that is used to complete the last term in the sum.

But what about the “incoming carry bit” into the Cuccaro adder? We don’t have
any incoming carry value here, but that qubit is crucial to the adder, so we
cannot simply get rid of it (normally when there is no incoming carry, the
Cuccaro adder uses an ancilla initialized to zero as the incoming carry).
Instead we replace this incoming carry with a dirty qubit—a qubit of unknown
state, which we guarantee will be returned to its initial state after use.
Dirty qubits are useful because they can be borrowed from idle parts of the
computation, avoiding the allocation of an ancilla [45, 22, 23]. By linearity,
if we can show that the correct operation is implemented with the dirty qubit
in the $\ket{0}$ basis state and in the $\ket{1}$ basis state, the correct
operation will be implemented using a dirty qubit in any state (even some
complicated entangled one). The Cuccaro construction automatically returns the
incoming carry to its original state, so we simply need to ensure that we can
apply the correct phase rotation. The case in which the dirty qubit is in the
$\ket{0}$ basis state is trivial—if the incoming carry is zero, the sum is
unchanged. The tricky case is when the dirty qubit is in the $\ket{1}$ state,
because the resulting sum will be one too large and then the phase rotation
based on the resulting register will be slightly too large. We propose a
simple fix: before starting the adder, perform a negative phase rotation
controlled off of the dirty qubit, that exactly cancels the extra phase that
will be applied due to the increase in the sum. This way, the total
accumulated phase will be exactly the product that was originally desired. We
show an explicit circuit for the application of the entire phase
$\phi(x_{0}+x_{1})(z_{0}+z_{1})$ using only one, dirty ancilla in Figure 2.

Finally, we show how to avoid the need for even an external dirty ancilla
qubit, making the multiplier entirely self-contained with zero ancillas. The
idea is simple: stop the sums one bit earlier, not including the most
significant bit of the inputs, and simply use one of these now-untouched bits
as the dirty qubit input to the adder. Doing so will remove their value from
the sum, but as was done with the carry bit, we can easily perform the portion
of the PhaseProduct corresponding to those topmost bits by directly
implementing the extra phase rotation via $m$ $CR_{\phi}$ gates. In fact, we
are already performing $m$ $CR_{\phi}$ gates controlled off of the dirty
ancilla to account for the extra phase it will add into the sum, so we can
simply adjust those phase rotations to include the qubit’s value when taken as
the top bit of the inputs as well.

###  C.2 Full case

There are a few ways we must generalize the above to cover all situations that
may arise in our algorithm: 1) the inputs may not be exactly the same length,
for example if $n$ is not perfectly divisible by $k$; 2) there may be several
terms in the sum, rather than simply two; 3) the terms of the sum may be
multiplied by a constant; and 4) the PhaseTripleProduct has three factors
instead of just two. It is easy to see how to handle the first three of these
generalizations, by simply allowing more of the product to be done directly
via $CR_{\phi}$ gates. As long as we ensure that the number of bits for which
we perform the PhaseProduct directly is sufficiently small, it will not affect
the asymptotic scaling. As shown in Figure 2, we simply take the “overlapping”
bits of the input and apply the algorithm from the previous sub-section to
implement their PhaseProduct, directly using $\mathcal{O}(n)$ $CR_{\phi}$
gates to implement the phase rotation proportional to the other bits that do
not overlap (a number of bits which is constant in $n$). We may use one of the
non-overlapping bits as the dirty qubit input into the Cuccaro adder. For sums
with multiple terms, we can simply apply the Cuccaro adder several times,
making sure to rotate by an amount proportional to the carry bit of each sum
when we have it available. By choosing the $w_{\ell}$ to be powers of 2 we
ensure that any coefficients multiplying terms in our sums are also powers of
two, and thus can be implemented as logical bit shifts in the inputs to the
adders.

Figure 3: A diagram of the bits of two values, here a 12-bit value $x_{0}$ and
a 10-bit value $x_{1}$ being summed to form the linear combination
$2x_{0}+x_{1}$. The setting is the implementation of a phase rotation
proportional to, say, $(2x_{0}+x_{1})(2z_{0}+z_{1})$ ($z$ registers not
shown). The purple bits are the “main” bits, which are summed together using
Cuccaro’s adder, with the top bit of $x_{0}$ (marked by diagonal hashes) being
used as the “incoming carry” into Cuccaro’s adder. Separately, a phase
rotation is applied controlled off of that top bit of $x_{0}$, but adjusted
for the fact that its value was carried into the sum of the main bits. The
polka-dotted bits represent other extra “non-overlapping” bits whose phase
rotations are implemented directly. Crucially, the number of extra bits is
constant in $n$, so the cost of performing the portion of the PhaseProduct
involving them is asymptotically negligible.

The PhaseTripleProduct can be implemented without ancillas via essentially the
same strategy, with one small hiccup: because there are three inputs, the
portion of PhaseTripleProduct involving a constant number of qubits of one of
the inputs requires
$\mathcal{O}(1)\cdot\mathcal{O}(n)\cdot\mathcal{O}(n)=\mathcal{O}(n^{2})$
gates to implement via the schoolbook algorithm of Eq. 9, which is too many to
maintain the larger algorithm’s asymptotic scaling. Fortunately there is an
easy workaround: a PhaseTripleProduct in which one of the inputs is a single
qubit is equivalent to a PhaseProduct controlled off of that single bit. So,
we may simply implement that portion of the operation via a constant number of
controlled fast PhaseProduct operations, making sure to set $k$ in such a way
that the asymptotic scaling of each PhaseProduct is better than the asymptotic
scaling of the overall PhaseTripleProduct algorithm. Thus the contribution of
these controlled PhaseProduct operations is asymptotically negligible and the
overall asymptotic scaling of the fast PhaseTripleProduct is maintained.

##  Appendix D In-place modular multiplication

The core of Shor’s algorithm for factoring an integer of $n$ bits can be
implemented via a series of in-place multiplications by classical constants,
controlled off of a single qubit (see the “one controlling qubit” trick, Sec.
2.4 of [46]). Written out, a single one of these multiplications consists of
the in-place operation

$\ket{x}\to\ket{cx\bmod N}$ (29)

for a classical integer $c$. Here we apply the modular multiplication
introduced in Sec. V to implement this operation to within an error $\eta$
using $2n+\mathcal{O}(\log 1/\eta)$ qubits (while maintaining the same
subquadratic gate count of the multiplication algorithm).

To do so, we make the following observation: our classical-quantum
multiplication algorithm does not require that the classical value $c$ be an
integer. It can be a floating-point number, which we may classically compute
to whatever arbitrary precision we desire before using it to compute the
values of the phase rotations in our multiplication algorithm. We use this
fact to perform roughly the following trick: compute the fractional value
$w=(cx\bmod N)/N$ up to some precision $m$, and then multiply by $N$ to
convert the fractional value $w$ into an integer. The accuracy of this
operation will depend on the precision to which we compute $w$; however, we
will find that we only need $\mathcal{O}(\log(1/\eta))$ extra bits to achieve
an error of less than $\eta$.

Algorithm 1 applies this idea to approximately implement the unitary of Eq.
29. It is clear by inspection that it uses $2n+\mathcal{O}(\log 1/\eta)$ total
qubits and $\mathcal{O}(n^{\log_{k}q})$ gates. In the following theorem we
prove the error bound.

Input : Quantum state $\ket{x}$ (extended to superpositions by linearity)

Classical constant $c$

Error level $\eta$

Output : Quantum state $\ket{cx\bmod N}$ (up to error $\eta$)

Let $m=n+\lceil 2\log(2+1/2\eta)\rceil$

Allocate a register of $m$ ancillas initialized to $\ket{0}$

1 Compute $\ket{x}\ket{0}\to\ket{x}\ket{w}$ for $w=((c-1)x\mod N)/N$ via
classical-quantum multiplication

2 Adding one ancilla to the top of the first register, compute
$\ket{x}\ket{w}\to\ket{x+Nw}\ket{w}$. State is now (approximately)
$\ket{cx\bmod N}\ket{w}$ or $\ket{(cx\bmod N)+N}\ket{w}$.

3 Using an ancilla qubit, compute whether the left register is greater than
$N$; subtract $N$ controlled by the ancilla. State is now $\ket{cx\bmod
N}\ket{w}$

4 Uncompute the ancilla qubit by computing whether $cx\bmod N<Nw$ via a
comparison operator

5 Subtract the value $w=((1-c^{-1})(cx)\bmod N)/N$ from the second register,
where $c^{-1}$ is the multiplicative inverse of $c\pmod{N}$. State is now
$\ket{cx\bmod N}\ket{0}$.

Algorithm 1 In-place classical-quantum modular multiplication

######  Theorem 1.

The final state $\ket{\psi}$ produced by Algorithm 1 has
$|1-\braket{\psi}{cx\bmod N}|<\mathcal{O}(\eta)$, for arbitrary inputs
$\ket{x}$ and $c$.

###### Proof.

We enumerate each step of the algorithm, at various points replacing the state
at that step with another state that is $\eta$-close (measured via inner
product as above). By the triangle inequality and the fact that unitary
transformations preserve inner products, we thus prove the theorem.

We begin with step 1. Applying Hadamard gates to the second register to
generate $\ket{x}\sum_{z}\ket{z}$, then a phase rotation of $\phi=2\pi wz$ and
then an inverse $\mathsf{QFT}$ modulo $2^{m}$ yields a state which in full
generality can be written

$\ket{x}\sum_{\tilde{w}}\alpha_{\tilde{w}}\ket{\tilde{w}}$ (30)

where the sum is over all $m$-bit binary fractions $\tilde{w}$. Let $S$ be the
set of $m$-bit binary fractions for which $|w-\tilde{w}|<\eta/2^{n}$. The key
to Algorithm 1 is that from quantum phase estimation (see Sec. 5.2.1 of [25]),
we have that $\sum_{\tilde{w}\in
S}\alpha_{\tilde{w}}\ket{\tilde{w}}\approx\sum_{\tilde{w}}\alpha_{\tilde{w}}\ket{\tilde{w}}$;
formally, the inner product of the actual state and the truncated state in
which only terms with $|w-\tilde{w}|<\eta/2^{n}$ are kept is within
$\mathcal{O}(\eta)$ of 1. Therefore, we may consider for the analysis of step
2 the truncated state for which $|w-\tilde{w}|<\eta/2^{n}$ for all terms in
the superposition.

Applying a quantum Fourier transform now to the register containing $\ket{x}$,
and then applying $\mathsf{PhaseProduct}(2\pi 2^{2n-m}/N)$ to the two
registers (where the second register’s value is now considered as an integer)
yields the state

$\sum_{z}e^{2\pi i(x+((c-1)x\bmod
N)z/2^{n+1}}\ket{z}\sum_{\tilde{w}}\alpha_{\tilde{w}}\ket{\tilde{w}}$ (31)

up to a phase error of at most $\mathcal{O}(\eta)$. Replacing that state with
the state having precisely the correct phase, and applying an inverse quantum
Fourier transform to the first register yields

$\ket{x+((c-1)x\bmod N)}\sum_{\tilde{w}}\alpha_{\tilde{w}}\ket{\tilde{w}}.$
(32)

The conditional subtraction of steps 3 and 4 is exact, yielding the state

$\ket{cx\bmod N)}\sum_{\tilde{w}}\alpha_{\tilde{w}}\ket{\tilde{w}}.$ (33)

Now performing step 5 by precisely inverting the operations of step 1 that
produced the state of the second register yields the all zero state on the
second register, completing the computation. ∎

For Shor’s algorithm, this operation must be repeated $\mathcal{O}(n)$ times.
Thus for each product we set $\eta=\mathcal{O}(1/n)$ so that the error of the
overall algorithm does not grow with $n$. This yields a total qubit count of
$2n+\mathcal{O}(\log n/\eta^{\prime})$, where $\eta^{\prime}$ is the error
over the entire algorithm.

##  Appendix E Estimation of gate and qubit counts

In this section, we describe the details of the resource estimates presented
in Tables 2 and 3 of the main text. The resource counts for previous works
were computed using code presented with those works; in this section we
describe the estimation for the algorithms presented in this work. The code
used to perform these estimates is available online:
https://zenodo.org/doi/10.5281/zenodo.10871109.

###  E.1 Recursive structure and choice of $k$

One of the most critical choices for efficiency is the choice of $k$, the
number of pieces into which we divide the values being multiplied. Because
increasing $k$ decreases the asymptotic cost but increases costant factors,
the obvious strategy is to find some “cutoff” value for each $k$, and use the
largest $k$ for which the desired product is above the cutoff. However, the
situation is complicated by the fact that in practice, the cost for a
particular $k$ does not increase smoothly with $n$—there is a range of
“overlap” where neighboring $k$ alternately jump past each other in
efficiency. (This is due to the fact that $k$ in general will not evenly
divide $n$). It is complicated even further by consideration that the linear
combinations like $(Ax)_{\ell}$ may be slightly different bit lengths for
different $\ell$, so it may be optimal to use different $k$ for different
$\ell$, even at the same level of recursion.

Because of all of those complications, for each recursive call we simply
exhaustively search through all $k$, and choose the construction that yields
the lowest cost (for some cost function that we may define; we use the number
of $CR_{\phi}$ gates as our cost function). In practice, we can perform this
exhaustive search in a reasonable amount of time because the tree depth is
logarithmic in $n$, and we memoize the optimal result for each value of $n$
encountered in the recursion.

###  E.2 Evaluation points $w_{\ell}$

We choose the $q$ evaluation points $w_{\ell}$ as follows. We always include
$0$ and $\infty$ because both can be evaluated directly from the inputs,
without the need for any linear combinations. Next we include $-1$ and $1$ if
needed, prioritizing $-1$ because the absolute value of the resulting linear
combination is bounded to a smaller value. For $q>4$, we then use the points
$-1/\omega,1/\omega,-\omega,\omega$ in that order, for $\omega$ powers of 2.
We use powers of 2 because then the resulting coefficients are powers of 2 and
thus values can be scaled by them via a logical bit shift. We prioritize unit
fractions over whole numbers because it helps reduce the size of the linear
combinations: the most-significant chunk (e.g. $x_{k-1}$) will in general be a
few bits longer than all the rest, because $n$ will in general not be evenly
divisble by $k$. Thus we would like to scale it by the smallest coefficient,
which occurs with unit fraction $w_{\ell}$ (see Appendix B). We prioritize
negative coefficients for the same reason as above with $-1$, that it helps
reduce the absolute value of the linear combination which reduces the ancilla
counts.

For pairs of $w_{\ell}$ that are the negative of each other, every negative
term is summed together first, for a value we may denote $x_{1}^{\prime}$.
Then $x_{-w_{\ell}}=x_{w_{\ell}}-2x_{1}^{\prime}$, requiring a single
addition/subtraction. Using this strategy, computing and uncomputing both
$x_{w_{\ell}}$ for a pair $\pm w_{\ell}$ requires a total of $2k-1$ additions,
while computing and uncomputing a single unpaired $w_{\ell}$ requires $2k-2$
additions. This strategy requires a few extra ancilla qubits, to hold the
overflow of both the intermediate sum register and the register holding the
final sum result.

In all cases, the register we choose to overwrite with the linear combinations
is the most significant (e.g. $x_{k-1}$), because it may be a few bits longer
than the others and thus require fewer ancilla qubits to store the value.

Next, we describe certain optimizations that are specific to the individual
applications described in the main text.

###  E.3 2048-bit classical-quantum multiplication

For non-modular arithmetic, the output register for this operation is twice as
long as the input. Thus, in that case we implement the whole PhaseProduct as
two simultaneous PhaseProduct operations, corresponding to the phase rotations
for the low and high halves of the output $z$ register respectively. We
perform both simultaneously so that we only have to compute the linear
combinations on the input $x$ register once rather than twice. This uses a
small amount more ancillas because the same ancillas cannot be reused for each
of the two halves of the output register; we consider this tradeoff worth the
reduction in gate counts.

The base case we use in this construction is very simple: we simply apply
$n_{\mathrm{base}}^{2}$ $CR_{\phi}$ gates, as in Eq. 8. As discussed in Sec. V
of the main text, it is possible to reduce the number of $CR_{\phi}$ gates
further by performing the base case product explicitly into an ancilla
register and applying $R_{\phi}$ gates on the result.

###  E.4 $x^{2}\bmod N$ proof of quantumness

A powerful feature of the proof of quantumness protocol is that “garbage bits”
can be discarded for free via measurement [6]. We use this fact to reduce gate
counts at the expense of a moderate number of ancillas, by storing the
$(A\mathbf{x})_{i}$ etc. in separate ancilla registers and then measuring them
away instead of using adder circuits to uncompute them. This corresponds to
the “fast” version of our circuits, listed in Table 3. The “balanced” version
computes and uncomputes the sums in-place, requiring a few more gates but many
fewer ancillas.

The “narrow” version uses the decomposition
$x^{2}z=z_{0}(x_{0}^{2}+x_{0}x_{1}+x_{1}^{2})+z_{1}(x_{0}^{2}+x_{0}x_{1}+x_{1}^{2})$
to reduce qubit counts quite a bit further, at the expense of some more gates.
In particular, it first computes the low half of the bits of the output
$x^{2}\bmod N$ via the first term above (that proportional to $z_{0}$),
performing half of the final QFT and measuring the result. Then, the same
$n/2$ qubits can be reused to compute the high half of the output bits and the
remaining half of the final QFT. This yields a total qubit count of roughly
only $3n/2$ qubits because all of the bits of the $n$-bit long output never
need to be stored at the same time. For the products themselves, we compute
the additions in-place as in the “balanced” construction to keep the qubit
counts low. (It is possible to take decompositions like this further to reduce
qubit counts even more, but in preliminary exploration we find that the
increase in gate counts is not worth the trade off).

For this construction, we use a more involved base case than in the classical-
quantum case. Notationally, denote the phase rotation we desire to implement
in the base case as $\phi^{\prime}x^{\prime 2}z^{\prime}$. Broadly, as
described in Sec. V of the main text, the idea is to replace the
$n_{\mathrm{base}}^{3}$ $CCR_{\phi}$ gates that would be required to implement
the schoolbook PhaseTripleProduct of Eq. 9 with $2n_{\mathrm{base}}^{2}$
$CR_{\phi}$ gates, by explicitly computing the product $x^{\prime 2}$ and then
performing $CR_{\phi}$ gates between the bits of that product and the bits of
$z^{\prime}$. However, because of the ability to freely uncompute garbage
bits, we can avoid having to allocate $2n_{\mathrm{base}}$ qubits to store
$x^{\prime 2}$ in its entirety. Instead, we decompose $x^{\prime
2}=\sum_{i}2^{i}\sum_{j}x^{\prime}_{j}x^{\prime}_{i-j}$—that is, we group the
bit products by the power of 2 they are scaled by. Following this grouping, we
first compute the sum for $i=0$, apply the $CR_{\phi}$ gates that involve the
least significant bit of the result, and then immediately measure it away to
uncompute it. Next we compute $i=1$, summing into any existing carry bits from
$i=0$, and again perform controlled rotations off of the least significant bit
before measuring it away. We continue this process for all $i$. Throughout
this process, we will need to store at most $\log_{2}n_{\mathrm{base}}$ carry
bits, a substantial improvement over the $2n_{\mathrm{base}}$ ancillas that
would be naively required.

Finally, we note that in the quantum-quantum case, the $k=2$ algorithm scales
as $\mathcal{O}(n^{2})$ just like the base case, but is observed to have worse
constants. Therefore we always move directly to the base case from $k=3$,
never using the $k=2$ construction.

##  Appendix F Parallel sequences

In Tables 4, 5, and 6, we provide examples of sequences of operations that
allow for the parallel computation of branches of the recursive tree, reducing
the overall depth of the circuit.

Operation | Register 0 | Register 1 | Register 2
---|---|---|---
(start) | $\ket{x_{0}}$ | $\ket{x_{1}}$ | $\ket{x_{2}}$
Add reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{1}+x_{2}}$ | $\ket{x_{2}}$
Add reg. 0 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{0}+x_{1}+x_{2}}$ | $\ket{x_{2}}$
Product on all registers | $\mathbf{\ket{x_{0}}}$ | $\mathbf{\ket{x_{0}+x_{1}+x_{2}}}$ | $\mathbf{\ket{x_{2}}}$
Invert sign of reg. 1 | $\ket{x_{0}}$ | $\ket{-x_{0}-x_{1}-x_{2}}$ | $\ket{x_{2}}$
Add reg. 0 to reg. 1 | $\ket{x_{0}}$ | $\ket{-x_{1}-x_{2}}$ | $\ket{x_{2}}$
Add $2\times$ reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{-x_{1}+x_{2}}$ | $\ket{x_{2}}$
Add reg. 1 to reg. 0 | $\ket{x_{0}-x_{1}+x_{2}}$ | $\ket{-x_{1}+x_{2}}$ | $\ket{x_{2}}$
Add reg. 0 to reg. 1 | $\ket{x_{0}-x_{1}+x_{2}}$ | $\ket{x_{0}-2x_{1}+2x_{2}}$ | $\ket{x_{2}}$
Add $2\times$ reg. 2 to reg. 1 | $\ket{x_{0}-x_{1}+x_{2}}$ | $\ket{x_{0}-2x_{1}+4x_{2}}$ | $\ket{x_{2}}$
Product on regs. 1 and 0 | $\mathbf{\ket{x_{0}-x_{1}+x_{2}}}$ | $\mathbf{\ket{x_{0}-2x_{1}+4x_{2}}}$ | $\ket{x_{2}}$
Invert sign of reg. 1 | $\ket{x_{0}-x_{1}+x_{2}}$ | $\ket{-x_{0}+2x_{1}-4x_{2}}$ | $\ket{x_{2}}$
Add $2\times$ reg. 2 to reg. 1 | $\ket{x_{0}-x_{1}+x_{2}}$ | $\ket{-x_{0}+2x_{1}-2x_{2}}$ | $\ket{x_{2}}$
Add reg. 1 to $2\times$ reg. 0 | $\ket{x_{0}}$ | $\ket{-x_{0}+2x_{1}-2x_{2}}$ | $\ket{x_{2}}$
Add reg. 0 to reg. 1 | $\ket{x_{0}}$ | $\ket{2x_{1}-2x_{2}}$ | $\ket{x_{2}}$
Add $2\times$ reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{2x_{1}}$ | $\ket{x_{2}}$
Divide reg. 1 by two | $\ket{x_{0}}$ | $\ket{x_{1}}$ | $\ket{x_{2}}$
Table 4:  $\mathbf{k=3}$ low-depth classical-quantum multiplication algorithm. This table lists the quantum operations performed to implement the unitary $\tilde{\mathcal{U}}_{q\times c}(a)$. The registers are divided into subregisters as $\ket{x}=\ket{x_{0}}\ket{x_{1}}\ket{x_{2}}$ and $\ket{z}=\ket{z_{0}}\ket{z_{1}}\ket{z_{2}}$ (using little-endian notation, so $x_{0}$ is the least-significant subregister). In this table only the state of the $x$ sub-registers are shown; the same operations are applied to the $z$ register. “Product on registers” means to apply a phase corresponding to the product of the respective $x$ and $z$ registers, usually by recursively calling the same algorithm again. Registers containing values upon which the algorithm is applied recursively are highlighted in bold. The linear combinations used here for the products correspond to the evaluation points $w_{\ell}\in\\{0,\infty,\pm 1,-2\\}$.  Operation | Register 0 | Register 1 | Register 2
---|---|---|---
(start) | $\ket{x_{0}}$ | $\ket{x_{1}}$ | $\ket{x_{2}}$
Add reg. 0 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{0}+x_{1}}$ | $\ket{x_{2}}$
Add reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{0}+x_{1}+x_{2}}$ | $\ket{x_{2}}$
Product on all | $\mathbf{\ket{x_{0}}}$ | $\mathbf{\ket{x_{0}+x_{1}+x_{2}}}$ | $\mathbf{\ket{x_{2}}}$
Add $-1\times$ reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{0}+x_{1}}$ | $\ket{x_{2}}$
Add $2\times$ reg. 0 to $-1\times$ reg. 1 | $\ket{x_{0}}$ | $\ket{x_{0}-x_{1}}$ | $\ket{x_{2}}$
Add reg. 1 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}-x_{1}}$ | $\ket{x_{0}-x_{1}+x_{2}}$
Add reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{2x_{0}-2x_{1}+x_{2}}$ | $\ket{x_{0}-x_{1}+x_{2}}$
Add $2\times$ reg. 0 to reg. 1 | $\ket{x_{0}}$ | $\ket{4x_{0}-2x_{1}+x_{2}}$ | $\ket{x_{0}-x_{1}+x_{2}}$
Product on regs. 1 and 2 | $\ket{x_{0}}$ | $\mathbf{\ket{4x_{0}-2x_{1}+x_{2}}}$ | $\mathbf{\ket{x_{0}-x_{1}+x_{2}}}$
Add reg. 0 to $2\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{4x_{0}-2x_{1}+x_{2}}$ | $\ket{3x_{0}-2x_{1}+2x_{2}}$
Add reg. 2 to $-2\times$ reg. 1 | $\ket{x_{0}}$ | $\ket{-5x_{0}+2x_{1}}$ | $\ket{3x_{0}-2x_{1}+2x_{2}}$
Add $2\times$ reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{0}-2x_{1}+4x_{2}}$ | $\ket{3x_{0}-2x_{1}+2x_{2}}$
Add reg. 1 to $-4\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}-2x_{1}+4x_{2}}$ | $\ket{-11x_{0}+6x_{1}-4x_{2}}$
Add $2\times$ reg. 1 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}-2x_{1}+4x_{2}}$ | $\ket{-9x_{0}+2x_{1}+4x_{2}}$
Add $8\times$ reg. 0 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}-2x_{1}+4x_{2}}$ | $\ket{-x_{0}+2x_{1}+4x_{2}}$
Add $2\times$ reg. 0 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}-2x_{1}+4x_{2}}$ | $\ket{x_{0}+2x_{1}+4x_{2}}$
Product on regs. 1 and 2 | $\ket{x_{0}}$ | $\mathbf{\ket{x_{0}-2x_{1}+4x_{2}}}$ | $\mathbf{\ket{x_{0}+2x_{1}+4x_{2}}}$
Add reg. 2 to $-1\times$ reg. 1 | $\ket{x_{0}}$ | $\ket{4x_{1}}$ | $\ket{x_{0}+2x_{1}+4x_{2}}$
Add $-1\times$ reg. 0 to reg. 2 | $\ket{x_{0}}$ | $\ket{4x_{1}}$ | $\ket{2x_{1}+4x_{2}}$
Add $-1\times$ reg. 1 to $2\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{4x_{1}}$ | $\ket{8x_{2}}$
Divide reg. 1 by 4 | $\ket{x_{0}}$ | $\ket{x_{1}}$ | $\ket{8x_{2}}$
Divide reg. 2 by 8 | $\ket{x_{0}}$ | $\ket{x_{1}}$ | $\ket{x_{2}}$
Table 5:  $\mathbf{k=3}$ quantum-quantum multiplication sequence. In this table only the state of the $x$ sub-registers are shown; the same operations are applied to the $z$ register. Registers containing values upon which the algorithm is applied again recursively are highlighted in bold. The linear combinations used here for the products correspond to the evaluation points $w_{\ell}\in\\{0,\infty,\pm 1,\pm 2,-1/2\\}$.  Operation | Register 0 | Register 1 | Register 2 | Register 3
---|---|---|---|---
(start) | $\ket{x_{0}}$ | $\ket{x_{1}}$ | $\ket{x_{2}}$ | $\ket{x_{3}}$
Add reg. 0 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{1}}$ | $\ket{x_{0}+x_{2}}$ | $\ket{x_{3}}$
Add reg. 3 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{1}+x_{3}}$ | $\ket{x_{0}+x_{2}}$ | $\ket{x_{3}}$
Add reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{0}+x_{1}+x_{2}+x_{3}}$ | $\ket{x_{0}+x_{2}}$ | $\ket{x_{3}}$
Add reg. 1 to $-2\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}+x_{1}+x_{2}+x_{3}}$ | $\ket{-x_{0}+x_{1}-x_{2}+x_{3}}$ | $\ket{x_{3}}$
Product on all regs. | $\mathbf{\ket{x_{0}}}$ | $\mathbf{\ket{x_{0}+x_{1}+x_{2}+x_{3}}}$ | $\mathbf{\ket{-x_{0}+x_{1}-x_{2}+x_{3}}}$ | $\mathbf{\ket{x_{3}}}$
Add reg. 1 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}+x_{1}+x_{2}+x_{3}}$ | $\ket{2x_{1}+2x_{3}}$ | $\ket{x_{3}}$
Add $-3\times$ reg. 3 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}+x_{1}+x_{2}+x_{3}}$ | $\ket{2x_{1}-x_{3}}$ | $\ket{x_{3}}$
Add $3\times$ reg. 0 to reg. 1 | $\ket{x_{0}}$ | $\ket{4x_{0}+x_{1}+x_{2}+x_{3}}$ | $\ket{2x_{1}-x_{3}}$ | $\ket{x_{3}}$
Add reg. 2 to $2\times$ reg. 1 | $\ket{x_{0}}$ | $\ket{8x_{0}+4x_{1}+2x_{2}+x_{3}}$ | $\ket{2x_{1}-x_{3}}$ | $\ket{x_{3}}$
Add $3\times$ reg. 3 to $2\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{8x_{0}+4x_{1}+2x_{2}+x_{3}}$ | $\ket{4x_{1}+x_{3}}$ | $\ket{x_{3}}$
Add $-1\times$ reg. 1 to $2\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{8x_{0}+4x_{1}+2x_{2}+x_{3}}$ | $\ket{-8x_{0}+4x_{1}-2x_{2}+x_{3}}$ | $\ket{x_{3}}$
Product on regs. 1 and 2 | $\ket{x_{0}}$ | $\mathbf{\ket{8x_{0}+4x_{1}+2x_{2}+x_{3}}}$ | $\mathbf{\ket{-8x_{0}+4x_{1}-2x_{2}+x_{3}}}$ | $\ket{x_{3}}$
Add reg. 1 to $-1\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{8x_{0}+4x_{1}+2x_{2}+x_{3}}$ | $\ket{16x_{0}+4x_{2}}$ | $\ket{x_{3}}$
Divide reg. 2 by 4 | $\ket{x_{0}}$ | $\ket{8x_{0}+4x_{1}+2x_{2}+x_{3}}$ | $\ket{4x_{0}+x_{2}}$ | $\ket{x_{3}}$
Add $6\times$ reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{32x_{0}+4x_{1}+8x_{2}+x_{3}}$ | $\ket{4x_{0}+x_{2}}$ | $\ket{x_{3}}$
Add $-15\times$ reg. 0 to $4\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{32x_{0}+4x_{1}+8x_{2}+x_{3}}$ | $\ket{x_{0}+4x_{2}}$ | $\ket{x_{3}}$
Add $15\times$ reg. 3 to reg. 1 | $\ket{x_{0}}$ | $\ket{32x_{0}+4x_{1}+8x_{2}+16x_{3}}$ | $\ket{x_{0}+4x_{2}}$ | $\ket{x_{3}}$
Divide reg. 1 by 2 | $\ket{x_{0}}$ | $\ket{16x_{0}+2x_{1}+4x_{2}+8x_{3}}$ | $\ket{x_{0}+4x_{2}}$ | $\ket{x_{3}}$
Add $-15\times$ reg. 0 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{0}+2x_{1}+4x_{2}+8x_{3}}$ | $\ket{x_{0}+4x_{2}}$ | $\ket{x_{3}}$
Add reg. 1 to $-2\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}+2x_{1}+4x_{2}+8x_{3}}$ | $\ket{-x_{0}+2x_{1}-4x_{2}+8x_{3}}$ | $\ket{x_{3}}$
Product on regs. 1 and 2 | $\ket{x_{0}}$ | $\mathbf{\ket{x_{0}+2x_{1}+4x_{2}+8x_{3}}}$ | $\mathbf{\ket{-x_{0}+2x_{1}-4x_{2}+8x_{3}}}$ | $\ket{x_{3}}$
Add $6\times$ reg. 1 to $-2\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{x_{0}+2x_{1}+4x_{2}+8x_{3}}$ | $\ket{8x_{0}+8x_{1}+32x_{2}+32x_{3}}$ | $\ket{x_{3}}$
Divide reg. 2 by 8 | $\ket{x_{0}}$ | $\ket{x_{0}+2x_{1}+4x_{2}+8x_{3}}$ | $\ket{x_{0}+x_{1}+4x_{2}+4x_{3}}$ | $\ket{x_{3}}$
Add $-3\times$ reg. 2 to $2\times$ reg. 1 | $\ket{x_{0}}$ | $\ket{-x_{0}+x_{1}-4x_{2}+4x_{3}}$ | $\ket{x_{0}+x_{1}+4x_{2}+4x_{3}}$ | $\ket{x_{3}}$
Add $3\times$ reg. 0 to $4\times$ reg. 1 | $\ket{x_{0}}$ | $\ket{-x_{0}+4x_{1}-16x_{2}+16x_{3}}$ | $\ket{x_{0}+x_{1}+4x_{2}+4x_{3}}$ | $\ket{x_{3}}$
Add $12\times$ reg. 3 to reg. 2 | $\ket{x_{0}}$ | $\ket{-x_{0}+4x_{1}-16x_{2}+16x_{3}}$ | $\ket{x_{0}+x_{1}+4x_{2}+16x_{3}}$ | $\ket{x_{3}}$
Add $-3\times$ reg. 0 to $4\times$ reg. 2 | $\ket{x_{0}}$ | $\ket{-x_{0}+4x_{1}-16x_{2}+16x_{3}}$ | $\ket{x_{0}+4x_{1}+16x_{2}+64x_{3}}$ | $\ket{x_{3}}$
Add $48\times$ reg. 3 to reg. 1 | $\ket{x_{0}}$ | $\ket{-x_{0}+4x_{1}-16x_{2}+64x_{3}}$ | $\ket{x_{0}+4x_{1}+16x_{2}+64x_{3}}$ | $\ket{x_{3}}$
Product on regs. 1 and 2 | $\ket{x_{0}}$ | $\mathbf{\ket{-x_{0}+4x_{1}-16x_{2}+64x_{3}}}$ | $\mathbf{\ket{x_{0}+4x_{1}+16x_{2}+64x_{3}}}$ | $\ket{x_{3}}$
Add reg. 2 to reg. 1 | $\ket{x_{0}}$ | $\ket{8x_{1}+128x_{3}}$ | $\ket{x_{0}+4x_{1}+16x_{2}+64x_{3}}$ | $\ket{x_{3}}$
Divide reg. 1 by 8 | $\ket{x_{0}}$ | $\ket{x_{1}+16x_{3}}$ | $\ket{x_{0}+4x_{1}+16x_{2}+64x_{3}}$ | $\ket{x_{3}}$
Add $-4\times$ reg. 1 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{1}+16x_{3}}$ | $\ket{x_{0}+16x_{2}}$ | $\ket{x_{3}}$
Add $-1\times$ reg. 0 to reg. 2 | $\ket{x_{0}}$ | $\ket{x_{1}+16x_{3}}$ | $\ket{16x_{2}}$ | $\ket{x_{3}}$
Divide reg. 2 by 16 | $\ket{x_{0}}$ | $\ket{x_{1}+16x_{3}}$ | $\ket{x_{2}}$ | $\ket{x_{3}}$
Add $-16\times$ reg. 3 to reg. 1 | $\ket{x_{0}}$ | $\ket{x_{1}}$ | $\ket{x_{2}}$ | $\ket{x_{3}}$
Table 6:  $\mathbf{k=4}$ quantum-quantum multiplication sequence. In this
table only the state of the $x$ sub-registers are shown; the same operations
are applied to the $y$ and $z$ registers. Registers containing values upon
which the algorithm is applied again recursively are highlighted in bold. The
linear combinations used here for the products correspond to the evaluation
points $w_{\ell}\in\\{0,\infty,\pm 1,\pm 1/2,\pm 2,\pm 4\\}$.