Integer multiplication below n log n

1 Introduction

Schönhage and Strassen proved in 1971 that two n-bit integers can be multiplied in O(n\log n\log\log n) time on a multitape Turing machine, and proposed n\log n as the optimal order of growth [25, Section 1]. After a sequence of improvements, Harvey and van der Hoeven reached the unconditional O(n\log n) upper bound [15]. We obtain a strict power saving in the logarithmic factor, in the same fixed finite-tape model. The algorithm uses linear combinations of intermediate data to reduce the cost of both address rearrangement and Fourier-transform layers.

For a binary string x of length n, let \operatorname{val}(x) be its nonnegative integer value, with the leftmost bit most significant. For 0\le z<2^k, let \operatorname{bin}_k(z) denote its binary representation padded on the left to length k. We consider exact multiplication: on input x\#y, with x,y\in\{0,1\}^n, the output must be \operatorname{bin}_{2n}(\operatorname{val}(x)\operatorname{val}(y)). Write \lg n=\max\{\lceil\log_2 n\rceil,1\}.

The bound is asymptotic through all input lengths. The constants and thresholds in the construction are extremely large; the algorithm is intended to establish a bit-complexity bound. The theorem refutes an eventual \Omega(n\log n) lower bound for every fixed multiplication machine in this model.

1.1 The two tape operations

Splitting the inputs into radix digits reduces multiplication to convolution of short integer coefficients. Our starting point is the Gaussian-resampling reduction of Harvey and van der Hoeven [15]: it replaces the Fourier transforms needed for this convolution by multidimensional transforms whose axis lengths are powers of two. After a further conversion to convolution, one coordinate is stored as the coefficient list of a polynomial modulo y^r+1, with r a power of two. Powers of y then serve as roots of unity in the remaining axes, and multiplication by such a root is a signed shift of the coefficient list. These are the synthetic polynomial transforms of Nussbaumer and Quandalle and Nussbaumer [21, 20].

In this representation, the transform calculation consists of address rearrangements, signed shifts, and two-point butterfly operations. Scanning an array of \Theta(n) stored bits at each of \Theta(\log n) transform levels already costs \Theta(n\log n). We therefore need savings in both movement and arithmetic. Two tape procedures supply them; in both bounds, V denotes the stored bit volume.

The first procedure interchanges two disjoint contiguous address fields of u bits in an array whose addresses form a complete Cartesian product. Its cost is O(Vu^\tau) for a fixed \tau<1. Although its final effect is an exact permutation, its intermediate steps form pointwise XOR combinations of data. The final permutation is exact for arbitrary initial auxiliary values.

The second procedure applies

simultaneously on selected coordinate bits of the numerical arrays. Each record is a polynomial with complex fixed-point coefficients, and H_0 acts coefficientwise. On the layouts used in the multiplication algorithm, one selects the same bit position in each of at most d equal-width address chunks, and the cost is O(Vd^{\lambda'}) for a fixed \lambda'<1. Intermediate values are exact Gaussian dyadics, whose real and imaginary parts have power-of-two denominators. A single truncation follows the completed contractive layer. Section 4 proves the address bound, and Section 5 gives the numerical layout, precision hypotheses, and layer bound.

1.2 Linear networks and the recursive saving

Each procedure arises from its own fixed linear network that exchanges two banks of values, with prescribed signs, and restores every auxiliary value. A wire’s fixed place in the network is called its role. To use the network on an array, regard a row as the complete ordered suffix, including all active address fields, at a fixed prefix and row index. Distribute rows cyclically among the W roles, padding the row count for each prefix to a multiple of W. Let V now denote the node volume after padding. All role streams have the same local address set and each contains V/W bits. A scalar gate is then applied pointwise to matching addresses on the roles that it touches.

The representation of a role’s array may change along the network. At a source, gate, or sink v, choose an invertible linear operator D_v on the space of arrays indexed by the common address set. We call D_v its frame: the stored array is D_vg when the logical array is g. All incidences of one gate have the same frame. Along an edge from u to v, the operator D_vD_u^{-1} changes the stored array from D_ug to D_vg. At a pointwise linear gate G, the common frame commutes with the gate:

The supplied input needs no preliminary change of frame: its logical interpretation is defined by applying the inverse input frame. Thus the intermediate frames cancel, and a route from input role w to output role \rho(w) applies D_{\mathrm{sink}(\rho(w))}D_{\mathrm{source}(w)}^{-1}, together with its prescribed scalar sign. Undoing the role exchange and those signs leaves the endpoint operator on every role array, including the roles that carried auxiliary scalar values.

For the bit procedure, that endpoint operator is an address shear: two fields, written as vectors H,D over a common residue ring, are updated by (H,D)\mapsto(H+D,D). Two further linear-cost address updates complete their interchange. For the numerical procedure, the endpoint operator becomes the simultaneous H_0 layer after diagonal phase operations and a dyadic scale. The intermediate frames are chosen so that every edge can be implemented by smaller calls to the same procedure.

The finite construction in Section 3 attaches subspaces of a fixed m-dimensional space to its gates and terminals. The subspaces at adjacent vertices are nested. Their dimension changes count the smaller operations needed on an edge, with an additional source contribution in the bit construction. The construction uses three-element subsets: gather and scatter create coefficients depending on intersection size, and auxiliary wires cancel the unwanted coefficients. Orthogonality of the corresponding indicator vectors permits the subspace assignments. The decisive property is that the total number s of smaller calls satisfies s<Wm.

A recursive call at parameter mf therefore uses s calls at parameter f, each on a stream of volume V/W. The parameter counts field digits for address interchange and selected axes for a numerical layer. After dividing cost by volume, multiplying the parameter by m increases the recursive contribution by s/W<m. This strict inequality is the source of the power saving. The later tape and arithmetic estimates preserve it, including the work of changing numerical address coordinates. Those changes are implemented by whole-chunk moves and controlled shifts, followed by a deterministic correction of the explicitly identified exceptional addresses.

1.3 From transforms to the exact product

The remaining proof places these procedures inside the multiplication reduction and accounts for every change of representation. First, the radix-digit coefficient list is placed on a product of distinct prime cyclic axes by a Chinese-remainder map that preserves convolution. Section 9 implements this map through triangular controlled shifts and charged axis moves. Gaussian resampling then expresses each source Fourier transform through a power-of-two tensor transform. Section 8 retains the source and target coordinate permutations in the transform identity, so they need not be executed as separate full sorts.

Next, Bluestein’s chirp identity converts the middle transform, including its target permutation, into normalized cyclic convolution. A phase twist changes the last cyclic coordinate into a polynomial modulo y^r+1; its negative wrap gives precisely the original cyclic wrap after the twist is removed. Sections 6 and 7 compute convolution on the remaining axes over this polynomial ring. The two forward synthetic transforms retain the same frequency order through their pointwise product, and the opposite transform consumes that order. Their polynomial products use signed packing and the established O(n\log n) multiplier as an ordinary subroutine.

Finally, the source permutation retained during resampling cancels across the two forward source transforms, their pointwise product, and the opposite-sign transform. Section 9 removes the known normalizations and digit scales, proves that the recovered coefficient error is strictly below 1/2, and propagates carries to obtain the exact binary product. A compatible choice of rational parameters makes all terms in its cost table fit Theorem 1. Figure 1 summarizes these reductions.

1.4 Other arithmetic operations

The classical Newton reductions described by Brent [5, Lemmas 2.1–2.4] compute reciprocal and square-root approximations at geometrically increasing precisions. Their multiplication costs are bounded by a geometric sum controlled by the last precision, and exact remainder and square tests recover the required integer answers. Section 10 gives the complete fixed-tape implementation.

1.5 Historical context and inherited methods

Karatsuba’s recursive construction, published jointly with Ofman in 1962, gave a binary multiplication circuit of size O(n^{\log_2 3}) [18]. Toom’s subsequent use of polynomial evaluation and interpolation gave circuit size n\exp(O(\sqrt{\log n})), and hence O(n^{1+\eta}) for every fixed \eta>0 [26]. These early results were formulated in terms of automata and logical networks. Cook adapted Toom’s method to the Turing-machine setting; the contemporary account of Schönhage and Strassen credits this adaptation explicitly [25, Section 1, p. 282].

In the multitape Turing model, Schönhage and Strassen proved the bound O(n\log n\log\log n) in 1971 [25]. Fürer improved it to n\log n\,2^{O(\log^* n)} in 2007, with a full journal account in 2009 [10, 11]; here \log^* n counts iterated logarithms until a constant threshold is reached. De, Kurur, Saha and Saptharishi developed a modular-arithmetic route to the same scale [7]. Later improvements made the exponential factor 8^{\log^* n} and then 4^{\log^* n} [17, 14]. Harvey and van der Hoeven proved the unconditional O(n\log n) bound [15].

The multitape bounds just cited concern that fixed model. Schönhage obtained linear-time multiplication on storage-modification machines [24]. Groff later gave a smaller-than-\log n factor in a unit-cost RAM model with constant-time arbitrary memory access, after polynomial preprocessing whose cost is excluded from that bound [13, Sections 2 and 4]. Neither result establishes the bound in the fixed finite-tape model of Theorem 1.

The transform methods behind these bounds also underlie the present construction. Schönhage–Strassen multiplication uses roots of unity in Fermat-type residue rings, where certain multiplications become shifts. Fürer combines cheap polynomial roots of unity with more general roots, reducing the frequency of expensive multiplications. Harvey and van der Hoeven use Gaussian resampling to pass to multidimensional transforms with power-of-two dimensions, then evaluate them using polynomial transforms developed by Nussbaumer and Quandalle and subsequently by Nussbaumer [21, 20, 15]. The Gaussian-gridding work of Dutt and Rokhlin is an earlier conceptual precursor [8]; Harvey and van der Hoeven discuss the connection in [15, Section 4.4.3]. The quantitative resampling and bit-cost estimates used here are those of the latter paper.

More specifically, the proof adapts Harvey and van der Hoeven’s Gaussian-resampling identity, contraction bounds, and short-convolution procedures [15, Theorem 4.2, Proposition 4.7 and Lemmas 4.8–4.12]. Their framework already permits arbitrary dimension. Here the source and target permutations remain in the transform identity, while the faster tape procedures perform the required axis movements. Bluestein’s conversion to convolution [4] and the synthetic-root convolution algorithm are retained; the latter follows the polynomial-transform construction in [15, Section 2.4 and Section 3]. Signed coefficient packing is a form of Kronecker substitution, with the fixed-tape version of [15, Lemma 2.5] as a direct antecedent. Calls to the established O(n\log n) multiplier are ordinary subroutine calls; they do not invoke the improved theorem recursively.

1.6 The finite-network ingredients

The subset-intersection construction in Section 3 has an earlier source outside multiplication. Alon’s presentation of the Frankl–Wilson construction [9] gives low-degree representations of complementary intersection graphs over different fields [2, Section 3, remark after the proof of Theorem 1.1]. At the prime p=2, its subsets have size three and the two polynomial evaluations are the intersection count modulo two and the intersection count minus one. These are the two pairing formulas used here. We change the ground-set size and express the second pairing through a rational bilinear form.

The gather and scatter operations, together with the side wires, have the form of linear index-code decoding: a low-rank linear summary gives the wanted coordinate plus contributions from known neighboring coordinates, and those contributions are subtracted. Bar-Yossef, Birk, Jayram and Kol establish the fitting-matrix formulation and this decoder [3, Section 3, proof of Theorem 5]. Lubetzky and Stav treat decoding over other fields and use related intersection representations to separate their possible ranks [19, Proposition 2.1 and Section 2.3]. These comparisons explain the algebraic role of the two fields here; they do not supply the subspace labels on the subsequent network.

Restoring scratch registers with arbitrary initial contents also has a precise precedent. The transparent-computation construction of Buhrman, Cleve, Koucký, Loff and Speelman subtracts an initial register value, performs a computation, adds the resulting value, and reverses the computation [6, Section 3]. Our scratch schedule is the linear-readout version of that construction. The rank reduction used for address shears is a singular-matrix Bruhat factorization: after reversing the coordinate order, Grigor’ev’s upper-triangular factors become the two lower-triangular factors needed here [12, Appendix A.2, Proposition 14(c)].

The finite networks and their labels are proved directly in Section 3. The additional requirements are the three-stage bank exchange, nested subspaces with a strict total rank or dimension deficit, and endpoint frames acting on every role. The tape implementations then turn that deficit into the volume-normalized recurrences of Sections 4 and 5, including the costs of addressing, padding, coefficient growth, and finite precision. The cited decoding, restoration, and factorization results provide ingredients; these further properties are established in the present proof.

1.7 Lower-bound scope and matrix transposition

Schönhage and Strassen conjectured n\log n as the optimal order of multiplication cost [25, Section 1]. In the present fixed-tape model, an eventual lower bound for every multiplication machine means that, for each fixed machine A, there are constants c_A>0 and N_A such that T_A(n)\ge c_A n\lg n whenever n\ge N_A. Its negation requires one machine satisfying \liminf_{n\to\infty}T_A(n)/(n\lg n)=0. The power saving proved here gives a limit of zero through all input lengths, which is a stronger conclusion.

Harvey and van der Hoeven consider a binary transposition machine supplied an integer m\ge1 and a k\times k binary matrix, where k=\lfloor\sqrt m\rfloor. They proved that if every fixed such machine has worst-case cost \Omega(m\lg m) through all sufficiently large supplied values of m, then every fixed multiplication machine has the corresponding eventual lower bound [16, Theorem 1.2]. Combining their quantitative reduction [16, Corollary 6.2] with Theorem 1 gives a deterministic fixed-tape algorithm that transposes a row-major k\times k binary matrix in O(k^2(\lg k)^{1-\kappa}) time, with \kappa as in Theorem 1. Section 11 proves the more general rectangular-matrix bound.

This comparison concerns unrestricted computation on a fixed number of one-dimensional tapes. Our address procedure forms XOR combinations of intermediate data, even though its final effect is a permutation. It therefore lies outside models restricted to moving data without such computations; see [16, Section 1.4] for the distinction.

1.8 Companion application to exact Fourier transforms

The complex phase network of Proposition 8 is also used in the companion manuscript [22, Section 1.1, Theorem 1.1 and Corollary 1.2]. That manuscript gives a single algorithm for the exact discrete Fourier transform at every positive transform length n. As n\to\infty, its cost is O(n(\log n)^{1-10^{-13}})=o(n\log n).

The cost is measured in an exact complex-register model that assigns unit cost to exact complex field operations. The algorithm accepts every x\in\mathbb C^n, with no bound on its entries or on intermediate complex magnitudes. Its computation on the input uses only addition, subtraction, and multiplication by prepared input-independent scalars. The magnitudes of those scalar coefficients are unrestricted as well.

Integer operations and random access on O(\log(n+2))-bit words have unit cost, with a fixed number of words for each address or integer. The model supplies one specified Fourier root \zeta_{D_*} of explicitly computable order D_*<1024n^3. Scalar preparation, schedule construction, and index work are charged. This is a separate application of the finite network, with its own uniform array implementation. It is neither a bit-complexity consequence nor a stable floating-point FFT consequence of the fixed-tape multiplication theorem, Theorem 1.

2 The machine model and elementary streams

A machine has one fixed finite alphabet and one fixed finite number of one-dimensional tapes. A step reads and writes at most one cell on each tape and moves each head by at most one cell. Fixed tracks and finitely many work tapes change time by a fixed factor. We distinguish moving one digit in a fixed radix, which needs a bounded number of scans, from interchanging fields whose numbers of digits grow. The latter permutation will be realized in Section 4 by computing XOR combinations of intermediate data streams. The present section supplies the elementary scans it uses. Every scan includes its counters and the return of its tape heads.

2.1 Array order and elementary scans

For a positive integer a write [a]=\{0,\ldots,a-1\}. An array with address set [a_1]\times\cdots\times[a_t] is stored in lexicographic order, with the first coordinate most significant. Its entries are bits unless a record width is specified. A trailing record of B bits can equivalently be regarded as a final address field [B]. An address is implicit in an entry’s position in this ordered stream; it is not stored beside every data bit. The logical volume V is the number of stored data bits. Work tapes and duplicate streams are charged in time but do not increase V in recurrences. All the fields discussed below are nonempty. A call with an empty array is handled directly.

An address permutation f moves the entry at x to f(x). Thus its action on an array A is (fA)(f(x))=A(x). Adjacent fields which are irrelevant to an operation may be replaced by a single field whose length is their product. In particular, the elementary procedures below have a fixed number of formal address fields even when a surrounding algorithm has a growing list of coordinates. Lengths and positions are supplied in canonical binary on descriptor tapes; no additional tape head is attached to an address coordinate. Each elementary call receives a fixed number of collapsed canonical lengths. Constructing them from a longer surrounding list is work of the caller; recursive calls below construct their own collapsed descriptors within their stated bounds.

When one address field is designated as a row field, collapse the fields before it to a prefix [P] and those after it to a suffix [B]. If the row field has range [R], a row at prefix p and row index r is the complete consecutive block \{p\}\times\{r\}\times[B]. Every row therefore contains B bits in the same suffix order.

For headers that concatenate a list of positive integer parameters, we use the self-delimiting code

The initial zeroes specify how many further bits follow the first one, so concatenated codes can be read in order. We use this code when an interface specifies a header; elementary calls still receive the fixed number of canonical binary lengths described above.

The next operations are linear because their target field follows their control fields. This order will be preserved whenever they are invoked.

When coefficients have a denominator d prime to q, their reductions modulo q^b can be computed by the extended Euclidean algorithm using ordinary tape arithmetic. The dimension and all original coefficients are fixed, so these preparations have polynomial cost in b. Nothing in Lemma 4 assumes unit-cost arithmetic.

For later use, fixed-width signed fixed-point binary addition, subtraction, sign change, and multiplication or division by a fixed power of two are streaming operations of linear bit cost, provided enough integer and fractional guard bits have been allocated for exact results. If an input record is most-significant-bit first, reverse it onto the arithmetic work tape. With least significant bits first on that tape, propagate carries using finite control. Reverse the completed record if its required external order is opposite. This convention separates the order of bits inside a record from the lexicographic order of records. Each required reversal is a linear pass and is charged.

3 Two finite networks with a rank saving

We construct two fixed linear networks that exchange two banks of values while restoring every auxiliary value. When the values are replaced by arrays, a common representation on all wires at a gate lets that gate act pointwise without disturbing the representation. We explain this invariant before constructing the subspace labels that make changes of representation cheap. The first network uses bit values and rational subspaces; the second uses Gaussian dyadic values and binary subspaces. The scalar domain specifies the values combined at a gate, while the subspace labels are used to construct operators on each wire’s address-indexed array.

Put

Let \mathcal T be the set of three-element subsets of [h]. There are two data banks X_a,Y_a, indexed by a=(a_1,a_2,a_3)\in\mathcal T^3. A wire carries one scalar. A gate is a fixed linear map on a specified set of wires; untouched wires retain their values. The vertices on a wire are its source, the gates that touch it in time order, and its sink; its edges join consecutive vertices. We call the name of a wire its role, so that we can distinguish a fixed place in the network from the value currently carried there.

3.1 Scalar operations and restored auxiliary values

There are three stages. At stage j, fix the other two coordinates and operate on the v values in each bank obtained by varying coordinate j. Thus there are v^2 separate invocations at each stage. Each such invocation has its own auxiliary wires, distinct from all other invocations, including those in other stages. There are

invocations. In a forward invocation, write x_T for the data used as controls and y_S for the data to be updated; S,T\in\mathcal T are their varying jth coordinates. These are the logical source and target banks of the invocation. For each ordered pair (S,T) that is neighboring under the rule below, there is an auxiliary wire A_{ST}. There are also central wires C_i, i\in[h], and in the second construction a wire C_\ast. We call all these auxiliary wires scratch wires; their initial values may be arbitrary.

The following maps specify the construction.

• •

  In the bit construction, scalars lie in \mathbb F_2, and S,T are neighbors when |S\cap T|=1. The copy map adds x_T to A_{ST} for each neighbor pair. The side injection adds \sum_{T:\,|S\cap T|=1}A_{ST} to y_S. The gather adds \sum_{T\ni i}x_T to C_i, and the scatter adds \sum_{i\in S}C_i to y_S.

• •

  In the complex construction, scalars lie in \mathbb Z[i,1/2], and neighbors are distinct triples with even intersection. Copy is unchanged. Side injection adds

  -\sum_{T\text{ neighboring }S}\frac{|S\cap T|-1}{2}A_{ST}

  to y_S. Gather adds \sum_{T\ni i}x_T to C_i and \sum_Tx_T to C_\ast. Scatter adds (\sum_{i\in S}C_i-C_\ast)/2 to y_S.

Each update leaves its control wires unchanged. The first two rows of the following schedule remove the effect of the old scratch values; the next four introduce the source values and add their effect to the target; the last two restore the scratch values. This cancellation of arbitrary initial scratch values is a linear instance of the restoration construction for transparent computation in Buhrman, Cleve, Koucký, Loff and Speelman [6, Section 3]. The forward schedule is

A fixed triple has

neighbors in the two constructions, respectively. The second count separates intersection zero and intersection two. Ordered pairs must be counted: there are vz side wires in each invocation, not vz/2. With c_{\rm b}=100 and c_{\rm c}=101 central wires per invocation, the total wire counts are

All these numbers are fixed, independently of any eventual input length.

3.2 From scalar values to arrays

We record the representation principle before choosing the labels. Let R be the scalar ring of either network and let \Omega be a finite address set common to all wires. Replace each scalar by an array in R^\Omega, and apply each scalar gate pointwise at every address. Write \mathsf S_\Omega for the resulting network on the vector of wire arrays.

At each source, gate and sink v, assign an invertible R-linear operator D_v on R^\Omega, called its frame. A frame acts on the values indexed by addresses within one wire’s array and may mix them; a gate instead combines wire values at each fixed address. All incoming and outgoing incidences of one gate use the same D_v. On an edge from u to v, apply D_vD_u^{-1} to that wire’s array. Define the operators on all input and output roles by

For a supplied physical input f, define its logical interpretation to be g=D_{\rm in}^{-1}f; this definition requires no physical preprocessing. At a vertex v, the invariant is that the stored array equals D_v applied to the logical array there. It holds at the sources by definition. On an edge whose logical array is \varphi, the edge operator changes D_u\varphi to D_v\varphi. At a gate with logical input arrays \varphi_j, R-linearity gives \sum_j c_jD_v\varphi_j=D_v(\sum_j c_j\varphi_j) for each output combination, with c_j\in R. Thus the common frame commutes with the pointwise gate, and the invariant continues to the sinks. The implemented operator is therefore

In particular, if the scalar network is a signed role permutation that sends input role w to output role \rho(w) with scalar sign \varepsilon_w, that route applies \varepsilon_wD_{\mathrm{sink}(\rho(w))} D_{\mathrm{source}(w)}^{-1} to its input array. All roles, including the auxiliary roles, enter this identity. Once the scalar routing is fixed, the endpoint frames determine the array operator on each route; the intermediate frames can then be chosen to make each edge change inexpensive.

Figure 2 summarizes this invariant.

3.3 Subspaces attached to the gates

The scalar network already has the required signed exchange. We now choose its intermediate labels to make the changes of frame inexpensive. A label is a nondegenerate subspace of a fixed bilinear space. On every edge one label will contain the other. For such an edge, the orthogonal residual is the orthogonal complement of the smaller label inside the larger one; its dimension is the absolute change in label dimension. The saving will come from keeping the total dimension lost on decreasing edges small relative to N.

For the bit network, let F=\mathbb Q^h with bilinear form

where J is the all-one matrix. This form has eigenvalues 1 and 1-h/9=-91/9, so it is nondegenerate. For a triple T let t_T be its indicator vector. Then

For the complex network, take instead F=\mathbb F_2^h with the ordinary dot product; here t_T\cdot t_T=1. Thus every triple line is nondegenerate, and neighboring triple lines are orthogonal, in the label field belonging to either construction. The scalar domain need not equal the label field. These degree-one intersection formulas specialize the construction described by Alon [2, Section 3, remark after Theorem 1.1] to the prime 2, with the ground-set size changed to h=100.

The ambient label space is \mathcal F=F^{\otimes3} with its tensor bilinear form. It has dimension m and is nondegenerate. Write

The prescribed terminal labels are

Write W for the wire count in the construction under consideration. The source dimensions sum to N and the sink dimensions to Wm-N. Thus the signed dimension changes will sum to Wm-2N, because the contributions at gates cancel. If decreasing edges lose a total dimension L, the sum of absolute changes is Wm-2N+2L. We will obtain L<N/2. The rational interface below will require an additional N source ranks, so this inequality will leave a strict saving in both constructions.

At stage j, fix an invocation. Put

An empty tensor product is the one-dimensional ground field; thus B=0 at stage 1. Let Q=\langle t_{a_{j+1}}\otimes\cdots\otimes t_{a_3}\rangle \subset F^{\otimes(3-j)}, again using the empty tensor product at stage 3. Since P is nondegenerate, A=B\mathbin{\perp}P. We prescribe the following transition for each data wire in the stage, where t is the indicator of that wire’s jth triple:

The outgoing Y label can also be written (P\otimes\langle t\rangle)^\perp_{A\otimes F}\otimes Q. These transitions expose one tensor factor at each stage. To see that they fit together, for j<3 write Q=\langle t_{a_{j+1}}\rangle\otimes Q'. The next stage has

Thus the prescribed outgoing labels are A'\otimes\langle t_{a_{j+1}}\rangle\otimes Q' on X and B'\otimes\langle t_{a_{j+1}}\rangle\otimes Q' on Y, exactly the next incoming labels. At stage 1 they start at U_a on X_a and zero on Y_a; at stage 3 they end at \mathcal F on X_a and U_a^\perp on Y_a, as prescribed at the terminals.

The gate labels are organized by the decomposition

with the common summand B\otimes F and the remaining summand P\otimes F of dimension h. Suppress the one-dimensional factor Q for the moment. On a side wire joining neighboring physical triples with indicators t_X,t_Y, the common component grows to B\otimes F, while the remaining component passes through

The first inclusion is possible because t_X\perp t_Y. For the central gates we assign A\otimes F on the physical X bank and B\otimes F on the physical Y bank. A central wire returns once from the former label to the latter, losing exactly P\otimes F. The gate order below confines all decreases to this return, while the other edges realize the prescribed stage transitions.

Number the successive times in an invocation by r=0,\ldots,7. At stages 1 and 3, time r uses row r of the scalar schedule (2). At stage 2 it uses the inverse of row 7-r, with the logical banks exchanged. In either case the gates touch the physical wires in the order

A side wire still joins neighboring physical triples in stage 2 because the neighbor relation is symmetric. The following table assigns the common gate labels. It suppresses the tensor factor Q; for a gate indexed by a physical bank triple, t denotes that triple’s indicator.

All labels are tensored by Q and viewed as subspaces of \mathcal F. Every displayed summand is nondegenerate.

All scratch sources have label zero and all scratch sinks have label \mathcal F, as prescribed above. The next lemma verifies the stage boundaries and all intervening edge changes, including these scratch edges.

The loss ratios have the exact values

In particular, both are smaller than 1/2. We now turn the dimension count into the two interfaces needed later.

3.4 The rational matrix interface

For a nondegenerate rational subspace U\subset\mathcal F, write P_U for the projection onto U along U^\perp. Assign P_U to each gate or sink with label U. Assign the source matrix -P_{U_a} to X_a and zero to every other source. These are fixed rational m\times m matrices; write M_v for the matrix at vertex v. Define \rho(X_a)=Y_a, \rho(Y_a)=X_a, and let \rho fix every auxiliary wire.

In the tape application, two address chunks are written as vectors H,D of m fields. For a matrix M, the frame \Phi_M moves the array entry at (H,D) to (H+MD,D), leaving the other address fields unchanged. The fields use a radix coprime to the fixed denominators, and addition is modulo their common range. Hence \Phi_M is an invertible \mathbb F_2-linear permutation of array entries, and \Phi_{M'}\Phi_M=\Phi_{M'+M}. Using D_v=\Phi_{M_v} in (5), the endpoint matrix difference on each routed role is exactly the shear applied to that role. The following identity makes it the same shear H\gets H+D on every role; Section 4 turns that shear into a field interchange and implements an edge change with one smaller interchange per unit of its rational rank.

3.5 The binary phase interface

Use the complex scalar network and binary labels. For a binary vector x, let \operatorname{wt}(x) denote its integer Hamming weight, and write [b]\in\{0,1\} for the integer representative of b\in\mathbb F_2. For a label U, let

The projection here is over \mathbb F_2. Let

where X interchanges the two coordinates. On functions of x\in\mathbb F_2^m define

The use of H in this definition is an identity of matrices. Each entry of \mathcal C_U and its inverse is a Gaussian integer divided by 2^m, so both preserve arrays over \mathbb Z[i,1/2]. The kernels implemented below have Gaussian dyadic coefficients.

To apply the scalar network to these functions, replace the scalar on each wire w by an array f_w:\mathbb F_2^m\to\mathbb Z[i,1/2]. At each address, a gate applies its scalar linear map to the values on the wires it touches. The frame in (5) is D_v=\mathcal C_U at a vertex with label U. On an edge from label U to label V, apply \mathcal C_V\mathcal C_U^{-1} to that wire’s array. All frames are diagonalized by the same matrix H^{\otimes m}. For an integer k\ge1, an address for k columns is (x^{(1)},\ldots,x^{(k)})\in(\mathbb F_2^m)^k; every frame and edge operator is then the tensor product of its one-column version over these k coordinates. Each scalar gate is still applied once at each joint address, to the vector of values on the roles it touches.

3.6 Explicit rational bounds for the two recurrences

The relative rank deficits are

For both recurrences we may use the fixed rational exponent

Both exact deficits displayed above are greater than 20\cdot2^{-50}. Since m<2^{20} and \log 2<1, we have \log m<20, whence

for either \eta=\eta_{\rm b} or \eta=\eta_{\rm c}. Consequently

Here \sigma is a rational upper bound for the logarithmic exponent of the complex recurrence.

4 Faster interchange of address chunks

We use the bit network to interchange two equal address chunks on tape. An edge matrix of rank a will require a smaller interchanges. Because the sum of these ranks is less than the number of wires times m, the resulting recurrence saves a power of the chunk width.

4.1 A matrix decomposition adapted to tape order

Lemma 4 updates a later field using earlier fields in linear time. Lower triangular changes of coordinates can be performed entirely in this direction. The required factorization is the Bruhat factorization for arbitrary matrices, with coordinate order reversed; see Grigor’ev [12, Appendix A.2, Proposition 14(c)]. We give its elementary elimination proof in this order.

4.2 Transferring a finite circuit to address permutations

Fix integers m,W\ge2. The circuit has W wires, including its scratch wires, with separate input and output terminals. Its fixed, pointwise bit gates induce a permutation \rho of the W scalar inputs at completion: the input value at role w appears at output role \rho(w). Assign a rational m\times m matrix to each terminal and to each gate, the latter assignment being common to every wire at that gate. A wire segment between consecutive assigned vertices is called an edge. For an edge e put

The required conditions are

All these data are finite and fixed. Proposition 7 supplies these conditions with W=W_{\rm b} and s=s_{\rm b}. Its negative source projection is essential to the endpoint identity and has already been included in the rank sum.

We use the rational value \tau=1-2^{-50} established in (10). Its verification uses fixed integer comparisons and the displayed exponential inequality; it needs no real-arithmetic oracle.

4.3 Removing width and row restrictions

We have proved the fast operation for power widths with enough rows. Rows for a general array can be supplied by a few of its own high digits. Their movement will cost less than the final bound.

The lemma permits arbitrary prefix, gap, and suffix lengths. Thus every later interchange of equal bit chunks is an instance of its stated layout, after multiplying the irrelevant contiguous lengths. For widths differing by one bit, choose the leading bit of the longer chunk as the extra bit. Suppress the outer prefix and suffix, and write g for the intervening address field. If the first chunk is longer, write it as (h_0,H), so H and the second chunk D have equal width. The physical orders are

The first step interchanges H,D by the lemma; the second moves the single bit h_0 immediately before H. If the second chunk is longer, write it as (d_0,D), with H the first chunk. Use

Here the first step moves d_0 immediately before H, and the second interchanges H,D. In each case g remains between the two full chunks in the final order. If H,D are empty, the indicated equal-width interchange is the identity. Each case uses one single-bit move, costing O(V) by Lemma 3, in addition to the equal-width interchange. No unrestricted array transpose is being used in these operations.

5 Fast simultaneous synthetic butterfly layers

We now apply the complex network to a common butterfly layer on several axes. The array has complete address set

with P,S\ge1, in lexicographic order. Its D consecutive K-bit chunks each select the bit at one common offset 0\le\rho<K, counted from the least significant end. Each address carries a polynomial record with r complex coefficients, initially in the unit disk on the grid 2^{-p}\mathbb Z[i], with \Theta(p) bits per real or imaginary component. The intended operation is the tensor product of H_0(u,v)=((u+v)/2,(u-v)/2) on the D selected bits, acting coefficientwise and preserving all other address fields.

The parameters will satisfy

The exponents and comparison constants are fixed, and we impose their further restrictions in Proposition 14. Two consequences explain the layout. The polynomial length r grows faster than every fixed power of p, so a scan of a record pays for polynomial work on its address. Also K/\log p\to\infty, which will make the exceptional set in a packed change of address bits small enough to repair by sorting. We will obtain cost O(Vd^{\lambda'}), with \lambda'<1 and V the logical bit volume, whereas separate butterfly factors would require D whole-array passes.

We first construct tensor products of the auxiliary two-point kernel C from Proposition 8. Its phase frames express a C tensor product using smaller groups of selected bits, separated by binary changes of basis. We implement those changes by moving whole address chunks. Splitting complete rows among the network’s wire roles then reduces the volume of each recursive call. Finally, we bound the intermediate dyadic values and convert the completed C layer to the normalized H_0 layer. The entire layer is evaluated exactly before one final truncation.

5.1 Phase frames for a simultaneous butterfly layer

Hold the unselected address bits fixed and regard the stored complex values on one wire role as a function of the selected bits. The phase identities below express an edge’s residual operator as a product of two-point kernels on selected directions.

Put

For the bit-flip matrix X=\bigl(\begin{smallmatrix}0&1\\1&0\end{smallmatrix}\bigr), define

Direct multiplication gives

In particular, C=iZC^{-1}Z. These identities use only Gaussian dyadic coefficients. The Hadamard matrices in the definition of \mathcal C_U are used to prove identities; they need not be executed by the machine.

Consider an edge from label U_- to label U_+, and choose an orthonormal basis v_1,\ldots,v_{\rho_{\rm edge}} of its residual. Proposition 8 gives the frame change

Here (X_vf)(x)=f(x+v), and we write C_v=aI+bX_v. The number of directional kernels is precisely \rho_{\rm edge}, including on decreasing edges. The next step puts those directions into coordinate positions without performing a general transpose.

Implementing an arbitrary residual basis. For M\in\operatorname{GL}(m,\mathbb F_2), let \mathcal P_M be the address permutation that moves the value at x to Mx. Thus (\mathcal P_Mf)(x)=f(M^{-1}x), and

Extend v_1,\ldots,v_{\rho_{\rm edge}} to a basis of \mathbb F_2^m and take these vectors as the first \rho_{\rm edge} columns of M. The displayed edge identity then becomes

where e_j is the jth standard basis vector. The matrix M need not preserve the dot product: equation (15) only uses its invertibility.

For the array layout used below, take a group of e=mf consecutive axis chunks and divide it into m consecutive slots of f chunks, numbered in physical order. Number the chunks within each slot by j=0,\ldots,f-1, starting at its least significant end. Write x_{h,j} for the selected bit in chunk j of slot h, where 1\le h\le m. Column j is the vector (x_{1,j},\ldots,x_{m,j}). A frame on this group is the tensor product of the same phase frame \mathcal C_U of Proposition 8 on these f columns.

To implement (16) on the group, first apply \mathcal P_M^{-1} separately to each column. For each 1\le h\le\rho_{\rm edge}, apply (C^{\epsilon_h})^{\otimes f} to the f bits x_{h,0},\ldots,x_{h,f-1}. This is one child operation on the f axes in slot h. Finally apply \mathcal P_M to every column. A negative directional kernel still requires only one forward child, because

The rightmost phase is applied first. All inverse wrappers are address signs and constant fourth-root phases, so they add no recursive calls.

Gaussian elimination expresses any M\in\operatorname{GL}(m,\mathbb F_2) as a product of at most m(m-1)+3m elementary row additions: one clears at most m-1 other entries per pivot, and each row swap is three row additions. Reversing the elimination gives the required product since each row addition is an involution. The number is a constant because m is fixed. Adding row h' to row h in every column applies x_{h,j}\leftarrow x_{h,j}\mathbin\oplus x_{h',j} for every 0\le j<f. Thus one row addition changes corresponding selected positions in two whole slots. The packed address operation proved below performs all of these XORs together.

The exact recursive contract. Apply the finite network with these column frames to a group of e=mf selected axes. Summing its residual dimensions gives exactly s_{\rm c} child calls, each on the f=e/m axes of one slot. The source and sink Z_{a,f} signs, together with the scalar i^{27f} at each affected sink, turn every endpoint kernel into a forward one. Negating the prescribed physical bank and undoing its exchange returns each result to its original role. Proposition 8 therefore gives C^{\otimes e} on every role, including roles whose initial scratch values were arbitrary. Formula (17) supplies each inverse child from one forward call, and none of the corrections adds a child. The remaining operations are a fixed number of binary row additions and coefficientwise scalar operations.

Returning to normalized butterflies. For an address x in a full layer, write x_1,\ldots,x_D for its selected bits. The diagonal operator \mathcal S_D multiplies the value stored at an address by the phase determined by those bits:

Here the unselected address bits are unchanged, and the same phase is applied to every coefficient stored at that address. Applying (14) on all selected coordinates gives

The scalar is

so it is a binary shift and fourth-root phase, with at most one further factor b. All phases, inverse wrappers, and scalar corrections above are exact Gaussian dyadic operations. Thus, once the physical binary basis changes and recursive stream layout are supplied, the complete parallel butterfly layer is exact up to its prescribed final truncation.

5.2 Packed changes of selected address bits

We implement a binary row addition on corresponding selected bits in two address chunks. A third chunk supplies temporary address bits; its arbitrary initial value must be restored. The construction first uses a constant number of cyclic shifts of whole chunk ranges, then corrects an explicitly bounded set of addresses. The correction depends only on addresses, so the time bound is worst-case over all payloads.

We use three previously established streaming operations. Lemma 12 exchanges two disjoint address chunks of equal width L in O(VL^\tau) time, where V is the current stream volume and 0<\tau<1 is fixed. The construction in Lemma 4 cyclically shifts the complete range of a target chunk by an offset determined by preceding fields in O(V) time apart from calculating the offsets. Finally, Lemma 3 applies any fixed permutation to the two-bit value at two specified address positions in O(V) time. Each operation uses a fixed number of one-dimensional tapes and restores the order of all spectator fields.

The asymptotic simplification in the following lemma describes a fixed family of parameters as p\to\infty, after the exponents and comparison constants have been fixed. The statement that R exceeds every fixed polynomial in p means that R/p^C\to\infty for each fixed C. Implied constants and eventual cutoffs may depend on the fixed choices and on when the family’s asserted bounds begin to hold. In particular, each explicit bound is uniform over all of its declared inputs; its asymptotic simplification is uniform over the allowed row counts and payloads within a fixed family.

Application to a fixed binary change of basis. An invertible matrix over \mathbb F_2 of fixed order m\ge3 is a product of a fixed number of elementary row additions; a row exchange is three such additions. In one invocation of the finite network on e=mf axes, slot h consists of f consecutive K-bit axis chunks. It therefore has address-bit width L=fK, and its selected bits x_{h,j} occupy positions \rho+jK for 0\le j<f. To add row h' to row h, use slot h' as x, slot h as y, and any third, distinct slot as z in Lemma 13. The lemma performs x_{h,j}\leftarrow x_{h,j}\mathbin\oplus x_{h',j} for every j in one operation and restores the third slot, whatever its initial address value.

Under the asymptotic hypotheses in the final assertion of that lemma, the complete change of basis and its inverse therefore cost O(V((eK)^\tau+1)) on a stream of volume V. Here e counts axis chunks, while eK counts address bits; the fixed factor m between eK and the slot width does not affect the bound. The same conclusion holds after the row padding and splitting below, because every role stream retains a complete suffix with the same fixed field lengths in every row.

5.3 Batching a common layer on a fixed number of tapes

We now turn the finite complex network into an algorithm for a whole layer. Throughout this subsection its constants are denoted by m,W,s: the network has W wires, its binary label space has dimension m, and the sum of the dimensions of its edge residuals is s<Wm. These are the constants of the complex network; they need not equal the constants used to prove the chunk-swap lemma.

Use the complete layout and parameter family specified at the start of the section, with fixed rational exponents c,\epsilon to satisfy the recurrence and guard conditions below. Store each polynomial in coefficient order with a common width w=\Theta(p) per real or imaginary component, and assume the complete address and its binary shape descriptor have length O(p). Thus an input record has 2rw bits. Write V for the current logical bit volume, excluding temporary work streams. We first count the work on C^{\otimes D} assuming one common exact coefficient format of width O(p); the guard calculation below then supplies that format.

An individual selected-bit kernel. Write each input coefficient once in this common exact format and retain it until the final truncation. The input width is \Theta(p), so the initial copy costs O(V) and the guarded stream has volume O(V).

To apply one C kernel, use the split from Lemma 3 at the selected address bit. Call the two streams 0 and 1 according to that bit. Because that bit has its complete two-element range, the two streams have the same remaining address set in the same order. Their aligned record addresses therefore differ only in the selected bit, and their coefficients are still ordered by the same polynomial index. Scan each pair of records in that order. For corresponding complex coefficients u from stream 0 and v from stream 1, compute

Addition, subtraction and halving use the fixed-point conventions of Section 2; multiplication by i swaps the real and imaginary components and changes one sign. The common guarded width makes these operations exact. The coefficient scan costs O(rp) per pair of records, linear in their guarded bit length and including any reversal of bits inside a coefficient. Write the first result to stream 0 and the second to stream 1, both in coefficient order, then merge according to the selected bit. The merge restores the original lexicographic address order. The elementary split and merge include their counters and cleanup and cost O(V); the paired arithmetic scan also costs O(V). Thus e individual kernels cost O(Ve) on any stream with these complete ranges and record format.

First consider an internal node containing e=m^k consecutive active chunks, with k\ge1. Split it into m consecutive slots, each containing e/m chunks. The phase construction supplies one recursive call on a slot for each vector in an edge residual basis, hence exactly s calls in one invocation of the network. If the logical volume entering this node is V, binary changes of basis, their inverses, scalar gates, endpoint corrections and reassembly of the physical banks cost

at this node. Each slot has (e/m)K address bits, so Lemma 13 gives the stated bound for a row addition. The network and all of its m-dimensional residual bases are fixed, so the total number of row additions is fixed as well. At the same guarded width, the scalar gates cost a constant number of scans of the current volume. Address and descriptor work bounded by any fixed polynomial in p per record is also absorbed: its ratio to the record’s \Theta(rp) bits tends to zero.

Where the rows come from. Put

If D\le q_0, use individual C kernels, at cost O(V\log d). Otherwise apply C individually on the selected bits of the first q_0 chunks, in O(V\log d) time. These chunks remain a complete consecutive address block before all the remaining chunks; regard their concatenated bits as one row index. Its range has length R_{\rm row}=2^{q_0K}\ge W^{k_0}, since K\ge1 for sufficiently large p and q_0K\ge k_0\log_2W. Within each preceding prefix, add the same number of zero rows to reach the next multiple of W^{k_0}; the volume grows by a factor at most two. At any recursive depth j\le k_0, the row count is divisible by W^{k_0-j}.

At a node, write a row index as u=Wg+w, with 0\le w<W, and send that entire row to physical wire role w, indexed there by g. Every role stream has the same remaining address shape. Thus a pointwise gate reads the records with the same g and the same remaining address in its touched roles. The scalar network allows these role values to be arbitrary, so the split supplies its scratch roles from existing rows. Each recursive edge call acts on one role stream with volume exactly 1/W of its parent’s logical volume.

Splitting the row index removes no address within a row. Every row in every role stream, including a padded row, retains the complete suffix [2^K]^{D-q_0}\times[S] in the same field order and with the same lengths; inactive chunks remain complete spectator fields. The row counts also remain equal across preceding prefixes, decreasing by the same factor W at each split. The network is invoked only on these equal rows: their common suffix gives the aligned role addresses and exact V/W child volumes above. For each packed call, collapse the preceding prefix and row index into one field, and collapse the spectator intervals around its three slots. This gives the rectangle of Lemma 13, whose guard count and deterministic repair apply to every child.

The W role streams use a fixed number of tapes. Before entering a child, push the other role streams onto a LIFO work tape, copy the active stream to a reserved child input area, and clear or reuse the role areas. On return, copy the result back and pop the parked streams, reversing each popped stream when necessary to restore its order. The stack head is at its top at every call boundary. Pushing, popping, copying, resetting and clearing at this node cost its own volume times a fixed constant; no step traverses an ancestor’s parked data. The call descriptor and program counter use a separate stack. The number of child calls is the fixed integer s, so these costs are included in the node overhead above. The fixed-tape chunk-swap subroutine has its own fixed work areas, reused on every call.

For each packed call, padding has increased the row count by at most a factor two and row splitting can only decrease it. Its active address still has O(p) bits. The four collapsed spectator lengths therefore have O(p)-bit descriptions. Its record-width parameter describes a polynomial of \Theta(rp) bits and has O(\log(rp))=O(p) bits as well. These lengths, together with K,f,\rho and the finite chunk-name permutation, are precisely the local descriptor required by Lemma 13. Both the record count and record width have O(p)-bit lengths, so its bound uses A_{\rm rect}=\lceil\log_2(2V)\rceil=O(p). A selected-position mask has at most DK=O(p) bits.

The binary descriptions of all growing shape parameters have polynomial length in p. Each nonempty child stream still contains a complete length-r polynomial record. Its volume therefore absorbs polynomial descriptor setup at that node, as well as polynomial address work per record. Even if every node is counted separately, there are only d^{O(1)} nodes: the branching factor s is fixed and the depth is at most \lceil\log_m d\rceil. Any additional fixed power of p from this count is still dominated by r.

Every completed network invocation applies C^{\otimes e} separately to each input row. Its selected bits lie in the remaining chunks, disjoint from the bits used by the earlier preprocessing. On the genuine rows, those two completed operations therefore commute. Each padded zero row is again zero when an invocation finishes, so the added rows can be deleted after the invocations on the remaining chunks have completed. Intermediate scratch values may be arbitrary.

Unrolling the recurrence. Fix \beta=1/2, and choose real exponents satisfying

The exponents can all be chosen rational, with a strict upper bound in the condition on \sigma. Stop recursion when e<d^\beta and apply the e C kernels individually. If F(e) denotes time divided by the current logical stream volume, then

Indeed, the s child calls are executed on streams of volume V/W, so their total cost is s(V/W)F(e/m). Dividing by the current volume V gives the factor s/W. At an internal node e\ge d^\beta, hence K\le d^c\le e^{c/\beta}. The overhead is therefore O(e^{\tau(1+c/\beta)}+1)=O(e^\lambda). Put a=s/W. Since a\le m^\sigma<m^\lambda, the internal costs sum to at most

If the leaves after j levels contain u=e/m^j<d^\beta active chunks each, their combined normalized contribution is at most

The first inequality uses a^j\le m^{j\sigma}=(e/u)^\sigma. It also covers a<1 by taking \sigma\ge0.

Write D-q_0=\sum_j a_jm^j in base m, with 0\le a_j<m. Partition the remaining chunks into a_j consecutive pieces containing m^j chunks for each j. There are at most (m-1)(1+\lfloor\log_m d\rfloor) pieces. Adding their costs and the row preprocessing gives

for any fixed

The case D\le q_0, in which preprocessing would consume all active chunks, was handled by individual kernels. Thus all constants and all tape counts are independent of p,d,r and the input data.

5.4 An explicit guard-width bound

The accelerated network is evaluated exactly until an entire normalized butterfly layer is complete. Its intermediate operators need not be contractions. We now give explicit constants that bound denominator and magnitude growth within one layer.

For the complex network, let

and set

For the stated h=100 construction, the complex network has

These values give m^{10}<B<m^{11}, hence \nu=11 and C_1=20. In fact the bound below needs only exponent fourteen; the choice twenty leaves an explicit additional margin.

We measure arithmetic dependency depth using complex additions, subtractions, halvings, sign changes and multiplication by i. For a group of e axes, let A(e) bound the number of these operations on any chain from an input coefficient to an intermediate coefficient in its recursive C-layer. A chain follows an operand into a later arithmetic result. Work on other coefficients or records that do not feed this value lies on other chains, even though the fixed-tape machine executes that work sequentially. To obtain a conservative bound, we concatenate all scalar gates and all recursive edge calls within one network invocation.

The scalar schedule (2) touches a data wire at most four times in each of the three stages. A scratch wire belongs to one invocation and is touched at most four times there. Thus every physical wire occurs in at most twelve scalar gates, and there are at most 12W gates. A gate has at most W inputs and outputs and coefficients in \{0,\pm1,\pm1/2\}. One output is a sum of at most W terms, each requiring at most one halving and one addition or subtraction. Evaluating all outputs from saved inputs therefore takes at most 2W^2 elementary operations.

For an inverse child, each Z^{\otimes f} is one sign chosen by the parity of the selected address bits, and the scalar fourth-root phase takes at most one sign change and one multiplication by i. This gives at most 4s operations for the child corrections. The source, sink and bank corrections take at most 4W+4 more. Calculating the address parities and phase exponents is address or descriptor work, not arithmetic on coefficients. Since s<Wm, the total 24W^3+4s+4W+4 is less than E=64(W+m+1)^3. Binary address permutations add no coefficient arithmetic at their interfaces.

The paired evaluation in (25) uses at most eight elementary operations: form the sum and difference, multiply the difference by i, then form the two half-sums. We may therefore use

The depth estimate concatenates the s calls, so it uses s. The factor 1/W in the earlier time recurrence accounted for their smaller stream volumes.

There are at most k\le\lceil\log_m d\rceil internal levels. For any bound x\ge1, sx+E\le(s+E)x=Bx. Starting with 8e\le8d\le8dB at a leaf and applying this inequality through the levels gives the first bound below. For the second, use B\le m^\nu and

Thus

For the sufficiently large sizes under consideration, the earlier piece count is at most md. Concatenating their depths contributes at most 8mB^2d^{\nu+2}. In either preprocessing branch, at most d individual kernels are used, contributing at most 8d. Splitting rows, padding, copying and removing padding add no coefficient arithmetic.

Each of the two outer diagonals applies a fourth-root phase and costs at most two elementary phase operations. Put t=\lfloor D/2\rfloor. Since b^2=-i/2, the scalar b^D is (-i)^t2^{-t}b^{D\bmod2}. Its phase costs at most two operations, and the possible factor b is computed as (z-iz)/2 in three operations. The shift by t consumes t fractional bits, so we charge it t units of denominator growth even though one tape scan performs the shift. Charging these same units to magnitude growth gives a common upper bound for both:

Here d\ge1, mB^2\ge1, and C_1\ge\nu+3.

To see how this count controls numerical size, start with disk-valued coefficients in 2^{-p}\mathbb Z[i]. A value at depth \ell lies in 2^{-(p+\ell)}\mathbb Z[i] and has modulus at most 2^\ell. This follows by induction: addition or subtraction can at most double the maximum modulus of its operands; halving adds at most one denominator bit; and signs and multiplication by i preserve both bounds. The charge of t for the final shift is the same as t successive halvings. Thus reserving p+\Delta fractional bits and \Delta+3 integer and sign bits per real component represents every intermediate coefficient exactly and prevents overflow. Shifts consume the reserved fractional zeros and discard no bits.

An address-swap subroutine may temporarily hold XOR combinations of the fixed-width encoding bits. These strings are not used as coefficients during the swap. It returns an exact permutation of the encodings before coefficient arithmetic resumes, so these internal bit operations cause no numerical rounding or denominator growth.

Choose \epsilon C_1<1. Since d=\Theta(p^\epsilon), \Delta=o(p), and the full coefficient width is O(p), as required in the time analysis. Addition, signed shifts, copying and phase changes consequently take O(p) bit operations per coefficient.

6 Synthetic transforms and their tape layout

We now implement the normalized Fourier transforms needed for convolution. Their entries will be short polynomials in a quotient where the required roots of unity act by exact signed permutations of coefficient lists. We also specify where each frequency entry is stored, because the opposite transform consumes that same order after a pointwise product.

The polynomial roots and radix-two butterflies belong to the fast polynomial-transform method of Nussbaumer and Quandalle [21, 20]; we use the form developed in [15, Section 2.4 and Lemma 3.2]. The work here is to apply the simultaneous layers on a suitable tape layout and retain the resulting frequency order through convolution.

For an integer r\geq 2 that is a power of two, put

We always use representatives of degree less than r. Multiplication by y^a, for any integer a, permutes the coefficients and changes some signs; consequently it is an isometry for this norm. For arrays over \mathcal R_r, the same notation denotes the maximum over all coefficients in all entries. A disk array has norm at most one. A p-bit grid array has real and imaginary coefficient parts in 2^{-p}\mathbb Z. For any array, let Q_p truncate each real and imaginary coefficient part toward zero to this grid. Then

6.1 Roots, butterflies, and stored frequency order

For a power of two t dividing 2r, define \omega_t=y^{2r/t}. For z\in\mathcal R_r^t and 0\leq j<t, set

Every power of \omega_t is an isometry, and each transform entry is an average of t such images. Thus both transforms are contractions. The ring has zero divisors, so we verify character cancellation directly instead of dividing by \omega_t^j-1.

Write

The length-one procedure is the identity. For t\geq2, the decimation-in-frequency step pairs z_k with z_{k+t/2} for 0\leq k<t/2. The two outputs of H_0 occupy these same slots, whose high address bit is b=0 or b=1, respectively. Multiply the b=1 output by \omega_t^{-k}, and then apply the length-t/2 procedure separately within the two branches. This computes F_t^- because, for b\in\{0,1\},

To follow the recursion, name the original address-bit slots of an axis of length t=2^a by a-1,\ldots,0, from most to least significant. For 0\leq s\leq a, after s levels the values of the first s named slots are branch bits b_0,\ldots,b_{s-1}. Equation (27) identifies these as frequency bits 0,\ldots,s-1, numbered from the least significant bit. The remaining slots encode an index in [0,u), where u=t/2^s. When s<a, the next level uses slot h=a-1-s to distinguish the paired inputs and then to hold b_s. Writing x_v for the bit in a still unprocessed slot v, the index used by its twiddle is

Previously written branch bits do not enter k. The same description holds in each shorter branch by (27), which proves the assertion inductively.

After all a rounds, the named slots from high to low therefore contain the frequency bits from low to high. If \operatorname{rev}_a(j) reverses the a binary digits of j, the entry (F_t^-z)_j occupies slot \operatorname{rev}_a(j). Let B_t denote this record permutation: it moves the entry at j to \operatorname{rev}_a(j). The recursion writes its outputs directly in this order; it does not move whole frequency indices.

The procedure consuming this stored order traverses the recursion in reverse: first process the two smaller branches, then apply the inverse monomial to the second branch, and finally apply H_0. Each forward level is invertible, since both its butterfly and its monomial are invertible. As H_0^2=I/2, the reversed level is half the inverse of the corresponding forward level. There are a=\log_2t levels, so the reverse procedure composed with B_tF_t^- is 2^{-a}I=I/t. Lemma 15 gives the same composite for F_t^+B_t^{-1} with B_tF_t^-. The forward procedure is invertible, so these two reverse operators agree. Thus the two exact procedures implement

respectively.

For D axes of lengths t_1,\ldots,t_D, put

A round applies one current H_0 on each participating axis. The twiddle for axis i depends only on that axis’s address bits. The two entries paired by another axis’s butterfly have the same axis-i address, so that twiddle commutes with the other butterfly. Consequently the individual axis steps may be grouped as a parallel butterfly \mathcal H, followed by one recordwise monomial operation J_E. At each address, J_E multiplies the record by

Here u_i is that axis’s current recursive length, b_i its newly written branch bit, and 0\leq k_i<u_i/2 is encoded by its still unprocessed named slots. The exponent depends on the address, not on coefficient values. In the opposite direction the grouped round is \mathcal HJ_E^{-1}: apply y^{-E} first and the parallel butterfly second. Each \mathcal H is a contraction and each J_E is an isometry, so every completed round is a contraction.

6.2 A layout on which the parallel layers apply

To apply Proposition 14 in a transform round, we must place the selected bits of the participating axes at a common offset in consecutive K-bit chunks. We now construct that layout, keeping each polynomial’s coefficients in a contiguous suffix. The layout uses the single-bit moves of Lemma 3 and the chunk interchanges of Lemma 12; each movement is included in the transform cost. The layer then supplies one truncated contraction per round, and its temporary rows and work streams are removed before the next round begins.

Pointwise multiplication commutes with applying the same record permutation LB to both operands. We may therefore multiply the forward outputs in their retained order and pass that order directly to the opposite procedure. The following polynomial product implements this operation with an established integer multiplier.

7 Exact normalized polynomial products

To multiply the polynomial records, we use the fixed-tape integer multiplier of Harvey and van der Hoeven [15, Theorem 1.1], with cost O(N\log(2N)) on N-bit inputs. This is an independent subroutine. Its inputs below have O(rp) bits, and we include the conversions between coefficient streams and signed integers in the cost. The reduction is a signed form of Kronecker substitution; compare [15, Lemma 2.5]. We specify the packing and recovery of coefficients to include negative values and the required tape format.

7.1 The normalized convolution interface

For arrays on \prod_i\mathbb Z/t_i\mathbb Z with values in \mathcal R_r, let * denote cyclic convolution in the displayed axes and define

For a fixed output coefficient, each of the M first cyclic addresses determines one second address, and each of the r first polynomial indices determines one second index modulo r. The negacyclic wrap supplies the sign. Thus that coefficient averages rM signed scalar products, and \mathcal M_r is a bilinear contraction:

If the other operand is a disk array, replacing one operand changes the output by at most the norm distance between its two values.

8 Resampling with the permutations left in the transform

The transform lengths needed for the integer product are pairwise coprime, whereas the synthetic transforms use powers of two. Gaussian resampling connects these two shapes. Its usual implementation sorts twice on each coordinate. We retain both permutations in the transform formula and compute that formula directly. This section proves the resulting interface, including its tape cost when the number of coordinates grows.

For an integer q\geq 1, write

All subscripts on a vector of length q are taken modulo q. The vector norm is the maximum modulus of its coordinates. All operator norms are induced by this norm: \|L\|=\sup_{\|u\|\leq1}\|Lu\|. Each normalized Fourier transform above is a contraction, since each output averages inputs multiplied by complex numbers of modulus one. Coordinate permutations preserve the norm. For an integer precision p\geq1, put

For storage between completed scalar routines, put w_{\rm c}=p+2. Represent 2^{-p}(a+ib)\in\mathbb D_p by a signed w_{\rm c}-bit two’s-complement word for a, followed by one for b, with each word written most significant bit first. A coordinate list is the concatenation of these fixed-length records in its stated order. Since |a|,|b|\leq2^p<2^{w_{\rm c}-1}, this format includes both signs of every permitted numerator, including -2^p and 2^p. An approximation \widetilde L:\mathbb D_p^a\to\mathbb D_p^b to a linear contraction L:\mathbb C^a\to\mathbb C^b has error E if 2^p\|\widetilde L u-Lu\|<E for every input in \mathbb D_p^a. The approximation itself need not be linear. Composing such algorithms adds their errors: at each step, insert the exact map at the computed input and use its contraction bound on the accumulated input error.

We use the established O(N\log N) integer multiplier of Harvey and van der Hoeven [15, Theorem 1.1] for the shorter integer records that occur below. We also use their Lemma 2.12: for integer precision p>100 and integers k>1, 0\leq j<k, one can compute \widetilde\zeta\in\mathbb D_p with 2^p|\widetilde\zeta-e^{2\pi i j/k}|<2 in time O(\max(p,\log k)^{1+\delta}), for every fixed 0<\delta<1/8. Throughout this section, numbered citations to that paper refer to its author-hosted 45-page manuscript.

8.1 The one-dimensional interface

8.2 Growing dimension and the padded array

For the tensor product, apply the one-dimensional maps to one coordinate line at a time. The chunk-swap algorithm supplies the movement of the physical axis fields, including its cost when the number of axes grows.

8.3 The middle transform and the final permutation

The factor QF_{\boldsymbol t} in (35) can be computed from a cyclic convolution and two pointwise phase multiplications. The chirp identity below performs the frequency change in Q without applying Q as a separate permutation of the array. The second identity shows how the retained source permutation R cancels when computing a convolution. Here * denotes cyclic convolution on the product of the coordinate groups and a dot denotes coordinatewise multiplication.

9 Exact integer multiplication and its cost

We now use the preceding algorithms to multiply two n-bit integers. An equal-width address swap on logical bit volume V costs O(Vu^\tau) for width u. A completed parallel normalized butterfly layer costs O(Vd^{\lambda'}) and has coefficient error less than \sqrt2\,2^{-p}. The ordered Gaussian resampling maps are contractions; their disk-grid approximations have error less than p^2 2^{-p} and cost O(tp^{3/2+\delta}\alpha) on a line of length t. Address permutations and intermediate dyadic arithmetic are exact. We choose compatible parameters for these procedures, recover the integer coefficients from their approximations, and account for the remaining tape operations.

9.1 A fixed rational choice

The two motif exponents were fixed in (10). Take the following fixed rational numbers:

These constants are deliberately conservative; no optimization is claimed.

The layer conditions hold with room to spare:

and

The strict final inequality absorbs the logarithmic number of groups in a base-m partition. The dimension exponent also satisfies

9.2 Input and transform sizes

The input is the promised string x\#y, with both operands written most significant bit first. In one left-to-right pass, copy x and y to two fixed work tapes and count the bits of x with a binary counter. The promise gives both operands this counted length n. The counter bound in Section 2 and the two sequential copies give cost O(n), including all leading zeroes. The work-copy heads finish just after the least significant bits, ready to move leftward when the digits are formed.

9.2.1 The padded box and working precision

After obtaining n, compute the size parameters below with the separate setup costs stated in this subsection. For sufficiently large n, put

Take the power of two T in [4n/b,8n/b), and let r=2^{\lceil(\log_2T)/d\rceil} and \ell=\log_2r. Choose d\ell-\log_2T of the d axis lengths to be r/2 and the others to be r. This gives

Indeed d\ell-\log_2T is an integer in [0,d), so at least one axis has length r; choose it as the last axis. These definitions give

The comparison constants can be fixed for every input length at once. For sufficiently large b, the bounds on n,T give b/2\leq\log_2T\leq b. Since p=6b, d=\lfloor b^\epsilon\rfloor and d\leq b, they imply

Thus we fix a_d=a_r=1/12, b_d=1, b_r=1/3, and C_\ell=1 in the layer and transform interfaces. Their pointwise bounds cover every computed pair (d,r), including all choices arising at the same p. These bounds also give that r exceeds every fixed power of p, K/\log p\to\infty, and \ell/K\to\infty. There are eventually enough axes for the O(\log d) reserved row axes and enough bits for all work blocks.

Set

The ceiling does not affect the needed bounds: for A\ge1, \lceil A\rceil\le2A, and hence

It follows that 2\le\alpha<\sqrt p eventually and \gamma=o(b). For example b\ge2^{24} implies \gamma\le b/4.

9.2.2 Prime source lengths

The distinct prime lengths s_i must be close to the powers of two t_i: their product must occupy a fixed fraction of the padded box, while each ratio t_i/s_i must leave enough separation for resampling. We verify both requirements for the parameters just chosen.

The prime-interval result [15, Lemma 5.1], based on explicit estimates of Rosser and Schoenfeld [23], gives the following bound: for 0<\eta<1/4 and x>e^{2/\eta}, the interval

contains at least \eta x/(2\log x) primes. Here \log is the natural logarithm. The size bounds in (39) give d=\Theta(p^\epsilon) and \log r=\Theta(p/d). For large n we have d\geq2, so 0<\eta<1/4, and

Thus both x=r/2 and x=r satisfy \log x>8d=2/\eta. For either choice, the stated lower bound for the number of primes is x/(8d\log x) and exceeds d: indeed, \log x=\Theta(p^{1-\epsilon}), whereas \log(8d^2\log x)=O(\log p)=o(\log x), so x>8d^2\log x eventually. The two intervals are disjoint because (1-\eta)r/2<(1-2\eta)r for \eta<1/4, and both intervals lie above 2. We may therefore choose distinct odd primes satisfying

Every \gcd(s_i,t_i) equals one because t_i is a power of two. Put S=\prod_i s_i. Multiplying the lower bounds on s_i and applying strict Bernoulli inequality gives

Together with s_i<t_i, this gives T/2<S<T. Writing \theta_i=t_i/s_i-1, the upper bound on s_i gives \theta_i\geq\eta/(1-\eta)>\eta. Therefore

The choice of \alpha gives \alpha=\Theta(p^{1/4+\epsilon/2}), and (40) already ensures 2\leq\alpha<\sqrt p for large n. We also have t_i\leq r\leq T<8n/b<2^{6b}=2^p. Thus every hypothesis of Lemma 19 holds for all sufficiently large input lengths. The superpolynomial growth of r established above also supplies the setup condition in Lemma 20.

The machine finds these primes by scanning the two intervals deterministically and testing each candidate by trial division through its square root. Testing every candidate would visit O(r) candidates and at most O(\sqrt r) trial divisors per candidate. Elementary division on O(\log r)-bit integers and management of a list of at most d selected primes give the conservative bound O(dr^{3/2}\operatorname{poly}(p)). Its logarithm is

Since p=6\lceil\log_2 n\rceil=\Theta(\log n), the search takes n^{o(1)}=o(n) time. All remaining scalar setup, including comparisons for fixed rational powers, modular inverses, prime products and shape descriptors, has cost polynomial in p. For example, d=\lfloor b^{2^{-75}}\rfloor is found by binary search using exact comparisons d^{2^{75}}\le b; the exponent is fixed. The same procedure computes K and all other fixed rational powers. All primes and length parameters are therefore generated by the machine; they are not advice.

9.2.3 The radix-digit polynomials

Write the two inputs as polynomials in 2^b, each with q=\lceil n/b\rceil nonnegative digits. Their product has degree at most

Thus convolution modulo x^S-1 has no nonzero wraparound. Every product coefficient obeys

since q\le n\le2^b.

9.3 An explicit cyclic-convolution layout

To use the source tensor, we first place the ordinary coefficient stream in its tensor order. Agarwal and Cooley used the Chinese remainder theorem to convert one cyclic convolution into multidimensional cyclic convolution [1]. The following map uses that algebraic split, with coordinate powers that make the change of order a sequence of controlled cyclic shifts.

Let s_1,\ldots,s_d be the pairwise coprime source lengths, let S=\prod_i s_i, and put

Here P_i denotes an integer prefix product, not a Fourier-coordinate permutation. The assignment

respects x^S=1, so it defines an algebra map. The target monomial group has size S, and the image of x has order S, because each \mu_i is a unit and the s_i are pairwise coprime. Its powers therefore run through all target basis monomials exactly once. Hence \Phi is an isomorphism and preserves cyclic convolution.

Every ordinary coefficient index has a unique mixed-radix expression

The coordinate of \Phi(x^k) in axis i is

Indeed terms with j>i are divisible by s_i, and \mu_iP_i=1 modulo s_i. Numerical order of k is lexicographic order of (a_d,\ldots,a_1). Insert the padding in one nested scan of the full box \prod_i[0,t_i) in this lexicographic order, with a_d the outermost counter and a_1 the innermost. At an address satisfying a_i<s_i for every i, copy the next original coefficient record; at every other address, write a zero record. The valid addresses occur in the original coefficient order, so each record is copied to its intended mixed-radix address. To delete the padding after the inverse axis reorder, use the same enumeration, copying the record at each valid address and skipping the record at each invalid address.

Store the nested binary counters consecutively on a fixed work tape. Their total width is \sum_i\log_2t_i=\log_2T=O(p), and the stored bounds s_i also have total width O(p). At each of the T addresses, updating the counters, testing validity, and handling the O(p)-bit record therefore cost O(p) steps, including counter-head returns. Each insertion or deletion scan costs O(Tp).

Now reorder the complete padded axes to (a_1,\ldots,a_d). In this physical order, the coordinates with indices below i precede coordinate i. Update the targets in the order i=d,d-1,\ldots,1 using (42). When target i is updated, none of its lower-index controls has yet changed, so the offset is computed from the original digits a_j required by the formula.

To specify the action on padded addresses, let

be sets of integer addresses. On a fixed control prefix, rotate the valid target interval [0,s_i) by the offset in (42) and fix [s_i,t_i). If any lower-index control is invalid, use the identity. This defines a bijection of \mathcal P: for fixed controls and spectators it is a cyclic permutation of the valid target interval and the identity elsewhere. Each step preserves \mathcal V and its complement. Consequently zero padding remains zero after every step. An inverse is explicit: process i=1,\ldots,d, subtracting the same expression formed from the already recovered lower coordinates; then undo the axis reorder and perform the deletion scan just described.

Each triangular update is a prefix-controlled cyclic shift. The controls precede its target in the current physical order. For each prefix, prepare the offset, split the valid target interval into its two consecutive pieces, and merge those pieces in the opposite order. Every target position carries its full spectator suffix as a block; the invalid interval is copied unchanged. There are at most T/t_i control prefixes for target i. For any fixed setup degree C, all offsets therefore cost O(dTp^C/r)=o(Tp) to prepare, since t_i\ge r/2 and r exceeds every fixed power of p.

The axis widths are \ell or \ell-1. To exchange two adjacent axes of unequal width, move the longer axis’s leading bit outside their equal-width pair, swap that pair by the proved chunk-swap procedure, and put the leading bit at its new position. This adds a fixed number of linear passes to the equal-width swap cost. Reversing all axes uses O(d^2) such moves. Therefore the map \Phi, its inverse, and their padding operations cost

This construction describes payload movement from address a to its new address. It is distinct from the operator convention (R_cz)_j=z_{cj} used for Fourier dilations, which moves a payload from address j to c^{-1}j.

9.4 Complex convolution by a coefficient twist

The power-of-two complex convolution used in the chirp identity is obtained from the synthetic convolution by a coefficient twist. Put \mathcal R=\mathcal R_r=\mathbb C[y]/(y^r+1) and \zeta=e^{\pi i/r}. The map

is well defined, since (\zeta y)^r=(-1)(-1)=1. On the bases (x_d^k)_{0\le k<r} and (y^k)_{0\le k<r} it is diagonal, multiplying coefficient k by the nonzero number \zeta^k. It is therefore an isomorphism, with inverse multiplying that coefficient by \zeta^{-k}. Both exact maps preserve the coefficient sup norm. For an explicit check on the wrap, write k+l=j+hr, where h\in\{0,1\}. A product term after twisting and negacyclic reduction has factor

Untwisting therefore recovers the cyclic coefficient with the correct sign. We next specify how this linewise map meets the ring-convolution interface. Put M=\prod_{i<d}t_i=T/r. For a complex array f on \prod_{i=1}^d\mathbb Z/t_i\mathbb Z, with t_d=r, define

This identifies the scalar tensor with an array of M records in \mathcal R. As maps of values, \iota and \iota^{-1} preserve norms and grids. Write (D_\pm f)_{\boldsymbol j,k}=\zeta^{\pm k}f_{\boldsymbol j,k}. The linewise extension of \Theta is \iota D_+, with inverse D_-\iota^{-1}. Consequently, with \mathcal M_r from Section 7, normalized complex cyclic convolution is

Here \mathcal M_r convolves the first d-1 axes over \mathcal R, whereas * on the right convolves all d scalar axes. The division by T is already part of \mathcal M_r. The maps D_\pm and \iota^{\pm1} are isometries, so \mathcal M_{\mathbb C} is also a bilinear contraction.

9.4.1 Computed complex convolution

For the tape representation, use the scalar format of Section 8, with w_{\rm c}=p+2 bits per component. For p\geq2, this satisfies p+2\leq w_{\rm c}\leq2p, so fix C_w=2 and instantiate the synthetic component width as w=w_{\rm c}. Both representations then concatenate exactly the same real and imaginary words in the same coefficient order: the last coordinate is already the contiguous polynomial suffix. Hence \iota and \iota^{-1} change only the shape descriptor and require no payload movement when these component formats agree.

For the synthetic transform calls, the machine writes the header

For each normalized polynomial product it supplies \Gamma(p)\Gamma(r)\Gamma(w_{\rm c}). These are precisely the headers of Lemmas 16 and 17, formed from the computed parameters. Their lengths are O(p); forming the transform header costs \operatorname{poly}(p)=o(Tp), and copying a product header for each of the M=T/r polynomial records costs O(Mp)=O(Tp).

For 0\le k<r, obtain disk-grid approximations \widetilde\theta_k^+ and \widetilde\theta_k^- to \zeta^k and \zeta^{-k}, each with error less than 2\,2^{-p}. The exponential routine uses denominator 2r and costs O(p^{1+\delta}) per factor. Write \widetilde{\mathcal M}_r for the procedure in Proposition 18. For disk-grid inputs f,g, one of norm at most .7, the computed complex convolution is the following sequence, with the twist arrays repeated on every last-coordinate line:

Each phase product is computed exactly as a dyadic product before Q_p truncates its two components. The factors and disk operands have numerators of magnitude at most 2^p; four signed integer products and their sums therefore use O(p) bits. Truncation shifts the absolute component numerators right by p places and restores their signs, then writes the resulting disk-grid value in the w_{\rm c}-bit component format. It costs O(p^{1+\delta}) per coefficient. Generating every factor anew therefore still gives total twist cost O(Tp^{1+\delta}). The precision calculation below proves that these calls return disk-grid arrays.

At the h' step, Proposition 18 returns its computed grid numerators after the exact scale M, in the known width w_M=w_{\rm c}+\log_2M. The disk margin proved below gives \|h'\|<.701 in every use here, so each component numerator has magnitude at most 2^p. Before applying \iota^{-1}, scan each w_M-bit word and remove only leading sign-extension bits, writing the same integer in w_{\rm c} bits. The retained sign bit agrees with every removed bit by that magnitude bound. This exact normalization costs O(Tw_M)=O(Tp) and supplies the matching scalar format for the output untwist; it introduces no numerical truncation.

9.5 Chirp and source-transform procedures

We next use the complex convolution to evaluate the middle factor QF_{\boldsymbol t} in (35), and then compose it with the Gaussian maps. For the chirp array a in (36), the phase numerator at tensor address j is computed modulo 2r as

All integers have O(p) bits and r/t_i\in\{1,2\}. Scan the stored axis descriptors and compute the d squares and products with the previously established integer multiplier. Additions, reductions modulo a power of two, and restarting local counters cost O(p) per term; multiplication costs O(p\log p)=O(p^{1+\delta}). Exponential evaluation then has the same O(p^{1+\delta}) bound. Thus chirp construction costs O(Td p^{1+\delta}) in total. In particular, this cost includes the scan of the axis descriptors at every coefficient.

Let \widetilde a be the resulting disk-grid array, so \|\widetilde a-a\|_\infty<2\,2^{-p}. For a disk-grid target tensor x of norm at most .6, the computed chirp transform forms the arrays

Thus \widetilde a is the convolution kernel and x_0 is the data operand after the first endpoint phase. Replacing the approximations by exact phases and exact convolution gives \overline a\mathbin{\cdot}\mathcal M_{\mathbb C} (a,\overline a\mathbin{\cdot}x)=QF_{\boldsymbol t}x by (36). Below we show that \|x\|_\infty\le.6 implies \|x_0\|_\infty\le.6, which supplies the smaller operand required by the computed convolution.

9.5.1 The source transform

For a source disk-grid tensor of norm at most .5, write \widetilde{\mathscr F} for the following procedure: apply the linewise approximation to \mathcal A, apply the computed chirp transform \widetilde{\mathscr C}, apply the linewise approximation to \mathcal B_0, and multiply the grid numerators exactly by 2^\gamma. Its exact counterpart is

Here an exact counterpart replaces every numerical subroutine by the exact map it approximates and omits truncations. The source arrays are stored in the padded power-of-two box, with zeros outside the valid source coordinates. The compressed scalar output uses w_{\rm c}-bit components. Appending \gamma zero bits to each word gives the exact scaled numerator in the known signed width w_{\rm c}+\gamma. For the stated input norm, the disk estimate below proves that every scaled coordinate still belongs to \mathbb D_p. Its component numerators therefore have magnitude at most 2^p, so a scan removes only leading sign-extension bits and returns the same integers in w_{\rm c} bits. This scan costs O(T(w_{\rm c}+\gamma))=O(Tp) and leaves the padded zero records zero. Define the opposite-sign procedure by

its exact counterpart is RF_{\boldsymbol s}^{+} by Lemma 21.

9.6 The multiplication algorithm

The source-transform procedures now supply the three transforms needed for the radix-digit convolution. Read each stored operand backward from its least significant end in groups of b bits. This produces the radix-2^b digits a_0,\ldots,a_{q-1} in increasing coefficient order. Buffer each group on a fixed work tape and read it back most significant bit first, reversing its least-significant-first input order; pad the final short group with leading zeroes. For digit a_j, write the exact grid numerator a_j2^{p-b-2} in a signed w_{\rm c}-bit real word, followed by a zero w_{\rm c}-bit imaginary word. This exponent is nonnegative because p=6b, and the represented value is exactly a_j2^{-b-2}. Append S-q zero coefficient records. Buffering and ordering the groups costs O(n), including the one partial group’s padding, and writing the S records costs O(Sp). Together with the initial count and copies, this input conversion costs O(n+Sp)=O(Tp); all n input bits, including leading zeroes, are retained before the final group is padded.

With these streams prepared, perform the following operations.

1. 1.

   Apply the coefficient-address map \Phi, with its padded layout, to obtain source arrays u,v. Their coefficients are exact p-bit grid values of norm less than 1/4.

2. 2.

   Compute \widetilde{\mathscr F}(u) and \widetilde{\mathscr F}(v).

3. 3.

   Multiply corresponding complex coefficients exactly and truncate each real and imaginary part toward zero to p fractional bits, writing the result in the common w_{\rm c}-bit component format. Denote the resulting source array by z.

4. 4.

   Compute \widetilde{\mathscr F}^{+}(z) and multiply its grid numerators exactly by S, retaining the denominator 2^p and writing each component in signed width w_{\rm c}+b. Denote the represented dyadic array by \widetilde w.

5. 5.

   Multiply the grid numerators of \widetilde w exactly by 2^{2b+4}S, again retaining the denominator 2^p, and write each component in signed width w_{\rm c}+4b+4. Round the represented real values to the nearest integers. At each valid address, write the rounded coefficient in one signed (3b+1)-bit word, most significant bit first; at every padded address, write a zero word of that width. Supply this local payload width to \Phi^{-1} and apply it, including removal of the address padding, to restore the ordinary coefficient order.

6. 6.

   Propagate radix-2^b carries in increasing coefficient order. Write the resulting bits from least to most significant on a work tape, then reverse that stream and output exactly 2n bits in most-significant-first order.

The two occurrences of S in this algorithm undo different normalizations. By (37), the exact counterpart after the three transforms is (u*v)/S^2, so after the first scale it is w=(u*v)/S. The second scale gives 2^{2b+4}(u*v), whose inverse image under \Phi is the original integer coefficient convolution. The next subsection verifies the disk bounds needed by every call and proves that the final error is less than 1/2; it also justifies the rounding, carries, and output length.

9.7 Precision and exact recovery

Proposition 14 uses C_1=20, so \epsilon C_1<1 makes its coefficient width O(p). Lemma 16 therefore applies to every synthetic transform here, with error less than \sqrt2\ell\,2^{-p}<L2^{-p}, where L=\log_2T, for all sufficiently large n.

9.7.1 Scales and disk margins

Put h_p=2^{-p}, and use the coefficient sup norm throughout this subsection. An error of a will mean absolute error less than ah_p. A linear contraction propagates an input error without increasing it. For a bilinear contraction \mathcal M, replacing its operands gives

When all four operands are in the disk, their two errors therefore add.

We also need the norm as well as the error of a computed phase product. Let e be a unit phase, let \widetilde e be a disk-grid approximation with |\widetilde e-e|<2h_p, and let x be a disk operand. The truncation bound in (26) gives

The second inequality follows from |\widetilde e|\le1 and the fact that Q_p does not increase the norm. It preserves every smaller operand norm through a twist or an endpoint phase multiplication. The same bounds hold coordinatewise for arrays of phases.

For disk-grid ring arrays f,g, put M=T/r. Two normalized synthetic forwards, pointwise polynomial multiplication divided by r, and the normalized opposite transform give (f*g)/(rM^2). The exact multiplication by M in Proposition 18 then gives \mathcal M_r(f,g)=(f*g)/T. Before that multiplication, the three transform errors and the polynomial-product truncation total less than 3L+2. Thus the error of the returned normalized convolution is bounded by

If one operand has norm at most .7, bilinear contraction gives \|\mathcal M_r(f,g)\|\le.7. Its computed grid output consequently has norm less than .7+5TLh_p<.701 for large n. This proves the disk margin after the scale M, without clipping the output.

Consider now \widetilde{\mathcal M}_{\mathbb C}(f,g) for disk grid operands, one of norm at most .7. Its two computed input twists f',g' are disk grid arrays, and the smaller one still has norm at most .7 by the phase inequality. The identification \iota preserves these bounds. The preceding ring margin therefore gives \|h'\|<.701. Neither \iota^{-1} nor the computed output untwist increases this norm. In particular, \widetilde{\mathcal M}_{\mathbb C}(f,g) is a disk grid array.

To compare it with \mathcal M_{\mathbb C}(f,g), insert the exact normalized \mathcal M_r at \iota f' and \iota g'. Its bilinear contraction bound propagates the two input twist errors, each less than 4. Its approximation contributes E_{\rm ring}, and the output untwist contributes less than 4. All these comparisons are at the normalized \mathcal M_r interface, after its factor M. They give

For a disk-grid chirp input x with \|x\|\le.6, the previously defined data operand x_0=Q_p(\overline{\widetilde a}\cdot x) has norm at most .6. The kernel \widetilde a is a disk grid array. Thus this pair satisfies the smaller-operand condition just proved: \|y_0\|<.701, and the last endpoint multiplication gives \|\widetilde{\mathscr C}(x)\|<.701. In particular, the chirp output is a legal disk-grid input to the compression maps.

For its error, the exact operands a and \overline a\cdot x are in the disk because every entry of a has modulus one. The kernel approximation contributes less than 2 and the first endpoint multiplication contributes less than 4 when inserted into the exact bilinear contraction \mathcal M_{\mathbb C}. The computed convolution contributes less than 6TL, and the last endpoint multiplication less than 4. Equation (36) identifies the exact result, so

Finally let the input to a source transform have norm at most .5. The computed expansion has error less than dp^2 and norm less than .5+dp^2h_p<.51. It is therefore a permitted chirp input. The chirp output is in the disk by the preceding argument, so each compression call is defined and returns another disk-grid array. Insert the exact contractions \mathcal A, QF_{\boldsymbol t}, and \mathcal B_0 at the computed inputs. The error before the source scale is less than 2dp^2+8TL. Multiplication by 2^\gamma therefore gives error less than 2^\gamma(2dp^2+8TL)h_p, which is bounded by E_sh_p with

because 2dp^2<TL eventually. The exact source transform RF_{\boldsymbol s} is a contraction, so the scaled approximation has norm at most .5+E_sh_p<.501. Its exact numerator scaling therefore also leaves a disk-grid output without clipping. Conjugation preserves the grid and norm, so the opposite procedure has the same error and disk guarantee. All these margins may be enforced by one fixed cutoff. For example, b\ge2^{24}, together with the other eventual construction conditions, makes dp^2h_p, 8TLh_p, E_sh_p and 28\,2^\gamma T^2Lh_p smaller than 10^{-3}.

9.7.2 The final integer coefficients

The scaled radix digits in u,v are exact p-bit grid values of norm strictly less than 1/4. Put

The exact arrays have norm less than 1/4, and the computed arrays have norm less than 1/4+E_sh_p<.26. Hence z=Q_p(\widetilde U\cdot\widetilde V_1) has norm at most .26^2<.5. To estimate its error, write

The factors \widetilde V_1 and U are in the disk. The two terms have total norm less than 2E_sh_p, and truncation adds less than 2h_p. Thus \|z-U\cdot V_1\|<(2E_s+2)h_p.

Let G:=RF_{\boldsymbol s}^+ be the exact opposite contraction. Lemma 21, applied to U=RF_{\boldsymbol s}u and V_1=RF_{\boldsymbol s}v, gives

The computed opposite transform may be applied to z because \|z\|<.5. Insert G at that computed input:

The first exact numerator multiplication by S gives \widetilde w, approximating w=(u*v)/S with error less than E_wh_p, where

Each coefficient of w is an average of S products of coefficients of u and v, so \|w\|<1/16. The earlier cutoff therefore gives

This proves the disk bound for the computed first scaled result.

In tensor order, the unscaled integer coefficient array is 2^{2b+4}S w. This second occurrence of S reverses the normalization still present in w, and 2^{2b+4} reverses the two input scales. Both multiplications act exactly on grid numerators. Since S\le T<2^b and L<b eventually, the final absolute error is less than

The true coefficients are real integers. If H_j is the real numerator after the second scaling, its represented value is H_j/2^p. The strict error below 1/2 makes the nearest integer to that value the unique correct coefficient.

These exact operations use known fixed record widths. If N is a signed w_0-bit integer, then -2^{w_0-1}\leq N<2^{w_0-1}. Since 0<S<2^b, its product satisfies -2^{w_0+b-1}<SN<2^{w_0+b-1} and therefore fits in signed width w_0+b. The first product by S starts at width w_{\rm c} and is written in width w_{\rm c}+b. The next product by S is written in width w_{\rm c}+2b; appending 2b+4 least significant zero bits then gives the exact second scale in width w_{\rm c}+4b+4. Each complex record still places its most-significant-bit-first real word before its imaginary word. The known widths determine every boundary in the next scan and are O(p) because p=6b.

The established integer multiplier performs the two products by S in O(p\log p) time per entry. To serialize each exact signed product at its stated width, zero-pad its magnitude and complement and increment for a negative sign. The proved bounds ensure that no significant bit is discarded. These conversions and the power-of-two shift use linear scans. To round H_j/2^p, recover the real magnitude, inspect its low p bits, and apply the rounding carry with the sign restored. The strict error bound makes the result the correct nonnegative coefficient, and (41) puts it in [0,2^{3b}). It therefore fits in one signed (3b+1)-bit two’s-complement word, written most significant bit first. Scan the padded box in the physical axis order recorded by its shape descriptor, using the same O(p)-bit nested counters and validity tests as above. Consume the imaginary word and write this real coefficient word at each valid address, with an all-zero word at each padded address. This serialization costs O(Tp) in total. It keeps the tensor address order and supplies payload width 3b+1 to the local descriptors for \Phi^{-1}, so the address routine receives a stream of known fixed-length records.

Lemma 17 proves exact signed packing, centered extraction and negacyclic reduction. Its only numerical error is the stated truncation after division by r, already included above. Undo the coefficient-address isomorphism \Phi before propagating radix carries.

To propagate carries, let Q_0=2^b. The recovered ordinary coefficients are the c_j in (41); they are zero for 2q-2<j<S. With q=\lceil n/b\rceil as above, use the recurrence for 0\le j<S:

The carries satisfy 0\le k_j\le q(Q_0-1)<Q_0^2. If this holds at j, then

which proves the next carry bound. The defining relation c_j+k_j=d_j+Q_0k_{j+1} also gives the exact invariant

It holds for j=0, and the defining relation proves its next instance. The unprocessed terms are multiples of Q_0^j, so the identity fixes the first j radix digits of the full product. The value k_j is the carry from the processed coefficients into the remaining positions.

After \Phi^{-1} the coefficients occur in increasing ordinary index, so this carry scan proceeds from least to most significant. Carry storage and work per coefficient are O(b) bits and steps. Write each digit’s bits from low to high on a work tape, and drain the final O(b)-bit carry in the same order. Copy the completed stream backward once to obtain most-significant-first order. The stream has length O(Sb)=O(Tp), including chunk padding and the final carry, so this reversal is a charged linear tape pass. The exact product is less than 2^{2n}; all bits above position 2n-1 are therefore zero. Remove only those zero padding bits and retain leading zeros as needed to produce exactly 2n bits. Zero inputs obey the same invariant and output rule.

9.8 The complete time bound

For the cost calculation put V=Tp=\Theta(n). The fixed-width coefficient formats use \Theta(p) bits per coefficient, so a full padded array has \Theta(V) stored bits. Table 1 lists costs after division by V; a fixed number of repeated transforms and convolutions is included in its constants.

The seven dominant powers in Table 1 have the following margins below one:

Direct substitution gives

Thus \min_i g_i=2^{-181}=2\kappa. This explicit slack will cover any remaining fixed power of \log p.

Here d^2 is absorbed by d^2\ell^\tau, and \log(rp)=O(p^{1-\epsilon}+\log p)=O(p^{1-\epsilon}). The first three rows are the terms in (30). The prefix-slot moves and individual butterfly rounds in Lemma 16 use O(dK) linear scans. Its positional layout uses O(p/K) exchanges of width-K chunks, because d\ell=O(p). Each exchange costs O(VK^\tau) by Lemma 12, so the normalized total is

For the at most \ell simultaneous rounds, Proposition 14 supplies O(Vd^{\lambda'}) per round. Its bound already includes the binary basis changes, deterministic repair, row padding and removal, all base-m groups, and fixed-tape cleanup.

Axis exposure and restoration require O(d^2) adjacent axis moves. Each costs O(V(1+\ell^\tau)): peel a single excess bit when widths differ by one, exchange their equal-width parts, and restore the bit. The d CRT rotations copy valid and invalid intervals in O(dV) time. Their offset setup was bounded above by O(dTp^C/r)=o(Tp), using at most T/t_i control prefixes for target i.

For Gaussian resampling, axis i has at most T/t_i processed lines. The line-machine costs therefore sum to

Dividing by V=Tp gives the fifth row. The factor \alpha=\Theta(p^{1/4+\epsilon/2}) explains its displayed power of p.

Chirp numerators have O(p) bits and can be computed modulo 2r as sums of d terms s_i j_i^2(r/t_i). Computing their squares and products with the previously established O(N\log(2N)) multiplier on N-bit inputs, then generating the exponentials, costs O(dp^{1+\delta}) per entry. The smaller twist and scalar pointwise-product costs fit the same row: they use a fixed number of exact products on O(p)-bit component numerators, followed by linear truncation scans. The two final products by S also cost O(p\log p) per entry, so their normalized cost is O(\log p)=O(p^\delta) as listed separately. There are T/r polynomial records in each pointwise ring product. Four integer products on O(rp) bits per record therefore cost O(Tp\log(rp)) in total. Signed packing, centered extraction, negacyclic wrapping, and copying are linear stream passes. These calls use the previously established multiplier on records of length O(rp)=o(n); they do not call the improved algorithm being constructed. The exact source scaling by 2^\gamma is a shift of O(p)-bit numerators by \gamma=O(p) places, so it is a linear scan per entry and fits the final-scaling row. The scalar-format adapters and the exact removal of sign-extension bits after the ring and source scales also use linear scans of their O(p)-bit components, within the listed scan costs.

Setup and exceptional streams across all calls. We next separate setup at recursion nodes from work repeated for every polynomial record or every Gaussian line. Here a visit means a node of the layer or address-layout recursions, or one of their round and group invocations. The internal recursion of the established integer multiplier and the iterations of the quoted Gaussian algorithms are already included in their separately charged running times. The layer and address-layout recursions have constant branching and depth O(\log p), and the number of their rounds and base-m groups is polynomial in p. Choose a fixed exponent A so that their total number of visits is at most p^A, and a fixed C bounding descriptor setup performed once per visit by p^C steps. This setup then costs at most p^{A+C}=o(Tp).

At a layer or synthetic-transform visit that processes polynomial records, there are at most O(T/r) such records, including the bounded row padding. Every relevant child keeps the whole length-r polynomial suffix. Consequently O(p^C) descriptor or address work per record, over these visits, costs at most O(Tp^{A+C}/r)=o(Tp). Setup before each Gaussian line call has a different count: there are at most 2dT/r lines per tensor map and only a fixed number of such maps, so O(p^C) setup per line costs O(dTp^C/r)=o(Tp). These estimates all use that r exceeds every fixed power of p.

The general address-swap subroutine also covers one-bit payloads, where no polynomial suffix is available. For root volume V_0 and current child volume V_j, the descriptor estimate in the proof of Proposition 11 gives (\log(2V_0))^{O(1)}=O(V_j). It uses the invariant that both original chunk ranges remain present in every child, and therefore includes polynomial descriptor work even for a one-bit payload. Elementary counter updates use local length copies, so they do not rescan a growing descriptor for each data bit. The chirp computation does scan the axis descriptors at every coefficient; its O(Tdp^{1+\delta}) cost was included in the chirp row.

For one packed change of basis, let V_{\rm call} be its current logical volume. The bad-address fraction is O(d2^{-K}), independently of the payload. Each polynomial record has \Theta(rp) payload bits, and an exceptional record adds an O(p)-bit destination key. Thus the compact exceptional stream has length

Stable bit sorting makes O(p) scans of this stream, for one-call cost O(V_{\rm call}dp2^{-K}). Extraction and reinsertion use a constant number of scans of the full current volume; those scans are already in the node overhead. Since there are at most p^A visits, each of volume O(V), the compact sorting costs sum to at most

for another fixed exponent A'. Here d=O(p^\epsilon) and K=\Theta(p^{\epsilon c}) with \epsilon c>0. Predicate and key arithmetic is the polynomial work per record counted above. This is a deterministic worst-case bound.

Finally, all role, arithmetic, and stack tapes form fixed finite families. A child parks and restores only its parent’s streams, with the stack head at the active top; copying and cleanup are charged to that logical volume. No operation repeatedly scans parked ancestors. Growing axis lists and call states are stored data, not extra tapes or alphabet symbols. Global setup and prime search cost o(V), as proved earlier.

Each of the seven displayed powers is at most 1-2\kappa by (45). Any fixed power of \log p is at most p^\kappa eventually. The complete large-input bound is consequently

Fix one cutoff beyond all the size, precision, margin and layout conditions, and use schoolbook multiplication below it. The finitely many smaller lengths are included by enlarging the implicit constant. With \lg n=\max\{\lceil\log_2n\rceil,1\}, the resulting single fixed-alphabet, fixed-tape deterministic machine has worst-case time O(n(\lg n)^{1-\kappa}) for every input length and produces precisely the required 2n-bit product.

10 Exact division and integer square root

Let n be the supplied input width, let a be the nonnegative dividend or radicand, and let b>0 be the divisor for the division task, with the values and output conventions of Corollary 2. Here p\ge1 denotes an integer target precision for a Newton routine.

11 A consequence for matrix transposition

For positive integers n_1,n_2,b, a row-major n_1\times n_2 array of b-bit strings stores entry (i,j), where 0\le i<n_1 and 0\le j<n_2, at position n_2i+j; transposition places this string at position n_1j+i. Harvey and van der Hoeven prove that any multiplication machine with worst-case cost M(m) yields a fixed multitape transposition machine with cost [16, Corollary 6.2]

The cited transposition interface receives m as well as n_1,n_2,b and the matrix stream [16, Section 1.1]. For the matrix problem defined here, compute m=n_1n_2b from the three binary descriptors and supply all four parameters. The four descriptors have total length O(\log(2m)), and the two schoolbook products that compute m cost O(\log^2(2m))=O(m) steps.

The cited multiplication interface likewise receives m explicitly and two nonnegative operands smaller than 2^m; it uses standard binary representations of the operands and product. On each such call, read the supplied m and form two words of exactly m bits, most significant bit first, by adding leading zeroes. Feed their x\#y concatenation to Theorem 1. Its final output pass writes the 2m-bit product consecutively. Scan this block to remove leading zeroes, emitting one 0 bit if every bit is zero. This fixes a canonical binary representation of the product. If a tape orientation requires reversal, read the stored word backward while writing it forward on a spare tape. Padding, normalization, and copying therefore cost O(m) steps on a fixed number of tapes.

Suppose the multiplication bound is O(m(\lg m)^{1-\kappa}) for a fixed 0<\kappa<1. Positivity of the dimensions and entry width gives \max(n_1,n_2)\le n_1n_2b\le m. Hence

The final relation uses 1-\kappa>0; the bounded cases are absorbed using the stated convention for \lg. Thus the cited transposition cost is O(m(\lg m)^{1-\kappa}). Taking m=n_1n_2b expresses this bound in matrix bits. In particular, an n\times n binary matrix can be transposed in O(n^2(\lg n)^{1-\kappa}) time. Dimensions equal to one cause no exception, since our convention gives \lg1=1.

A Packed additions on grouped rectangles

The simultaneous-layer construction uses the rectangular layout of Lemma 13. We record here a more general version in which a supplied descriptor lists groups with different spectator lengths. Each group is a complete rectangle, so the full rectangular algorithm, including its deterministic repair, can be applied to one group at a time. The only additional work is reading the group descriptors and passing consecutive groups through the same local tapes.