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 (Harvey and Hoeven 2021): 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 (Nussbaumer and Quandalle 1978; Nussbaumer 1980).

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 \[H_0(u,v)=\bigl((u+v)/2,(u-v)/2\bigr)\] 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.

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: \[G(D_vg_1,\ldots,D_vg_t)=D_vG(g_1,\ldots,g_t).\] 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.

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.

Other arithmetic operations

The classical Newton reductions described by Brent (Brent 1976, 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.

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})\) (Karatsuba and Ofman 1962). 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\) (Toom 1963). 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 (Schönhage and Strassen 1971, sec. 1, p. 282).

In the multitape Turing model, Schönhage and Strassen proved the bound \(O(n\log n\log\log n)\) in 1971 (Schönhage and Strassen 1971). Fürer improved it to \(n\log n\,2^{O(\log^* n)}\) in 2007, with a full journal account in 2009 (Fürer 2007, 2009); 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 (De et al. 2013). Later improvements made the exponential factor \(8^{\log^* n}\) and then \(4^{\log^* n}\) (Harvey et al. 2016; Harvey and Hoeven 2019). Harvey and van der Hoeven proved the unconditional \(O(n\log n)\) bound (Harvey and Hoeven 2021).

The multitape bounds just cited concern that fixed model. Schönhage obtained linear-time multiplication on storage-modification machines (Schönhage 1980). 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 (Groff 2019, sec. 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 (Nussbaumer and Quandalle 1978; Nussbaumer 1980; Harvey and Hoeven 2021). The Gaussian-gridding work of Dutt and Rokhlin is an earlier conceptual precursor (Dutt and Rokhlin 1993); Harvey and van der Hoeven discuss the connection in (Harvey and Hoeven 2021, sec. 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 (Harvey and Hoeven 2021, 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 (Bluestein 1970) and the synthetic-root convolution algorithm are retained; the latter follows the polynomial-transform construction in (Harvey and Hoeven 2021, sec. 2.4 and Section 3). Signed coefficient packing is a form of Kronecker substitution, with the fixed-tape version of (Harvey and Hoeven 2021, 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.

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 (Frankl and Wilson 1981) gives low-degree representations of complementary intersection graphs over different fields (Alon 1998, sec. 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 (Bar-Yossef et al. 2006, sec. 3, proof of Theorem 5). Lubetzky and Stav treat decoding over other fields and use related intersection representations to separate their possible ranks (Lubetzky and Stav 2009, 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 (Buhrman et al. 2014, sec. 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 (Grigor’ev 1982, 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.

Lower-bound scope and matrix transposition

Schönhage and Strassen conjectured \(n\log n\) as the optimal order of multiplication cost (Schönhage and Strassen 1971, sec. 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 (Harvey and Hoeven 2025, Theorem 1.2). Combining their quantitative reduction (Harvey and Hoeven 2025, 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 (Harvey and Hoeven 2025, sec. 1.4) for the distinction.

Companion application to exact Fourier transforms

The complex phase network of Proposition 8 is also used in the companion manuscript (OpenAI 2026, sec. 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.

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 \[ b_v=\lfloor\log_2v\rfloor+1,\qquad \Gamma(v)=0^{b_v-1}\operatorname{bin}_{b_v}(v),\qquad |\Gamma(v)|=2b_v-1\quad(v\ge1). \tag{1}\] 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.

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 \[I=3v^2=78440670000\] 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 \(j\)th 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 (Buhrman et al. 2014, sec. 3). The forward schedule is \[ \begin{array}{c|l|l} \text{row}&\text{operation}&\text{gate grouping}\\ \hline 0&\text{subtract side injection}&\text{one gate per target}\\ 1&\text{subtract scatter}&\text{one gate on targets and center}\\ 2&\text{copy}&\text{one gate per source}\\ 3&\text{gather}&\text{one gate on sources and center}\\ 4&\text{scatter}&\text{one gate on targets and center}\\ 5&\text{side injection}&\text{one gate per target}\\ 6&\text{undo gather}&\text{one gate on sources and center}\\ 7&\text{undo copy}&\text{one gate per source}. \end{array} \tag{2}\]

A fixed triple has \[z_{\rm b}=3\binom{97}{2}=13968,\qquad z_{\rm c}=\binom{97}{3}+3\cdot97=147731\] 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 \[\begin{align*} W_{\rm b}&=2N+I(vz_{\rm b}+c_{\rm b}) =177176569091445000000,\tag{3}\\ W_{\rm c}&=2N+I(vz_{\rm c}+c_{\rm c}) =1873807244643542670000. \tag{4}\end{align*}\] All these numbers are fixed, independently of any eventual input length.

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 \[D_{\rm in}=\bigoplus_w D_{\mathrm{source}(w)},\qquad D_{\rm out}=\bigoplus_w D_{\mathrm{sink}(w)}.\] 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 \[ D_{\rm out}\,\mathsf S_\Omega\,D_{\rm in}^{-1}. \tag{5}\] 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.

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 \[\langle x,y\rangle=x^{\mathsf T}(I-J/9)y,\] 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 \[\langle t_S,t_T\rangle=|S\cap T|-1, \qquad \langle t_T,t_T\rangle=2.\] 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 (Alon 1998, sec. 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 \[u_a=t_{a_1}\otimes t_{a_2}\otimes t_{a_3},\qquad U_a=\langle u_a\rangle.\] The prescribed terminal labels are \[\begin{array}{c|cc} \text{role}&\text{source label}&\text{sink label}\\ \hline X_a&U_a&\mathcal F\\ Y_a&0&U_a^\perp\\ \text{each scratch role}&0&\mathcal F. \end{array}\] 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 \[A=F^{\otimes(j-1)},\qquad P=\langle t_{a_1}\otimes\cdots\otimes t_{a_{j-1}}\rangle, \qquad B=P^\perp\subset A.\] 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 \(j\)th triple: \[\begin{array}{c|c|c} &\text{incoming label}&\text{outgoing label}\\ \hline X&A\otimes\langle t\rangle\otimes Q&(A\otimes F)\otimes Q\\ Y&B\otimes\langle t\rangle\otimes Q& \bigl((B\otimes F)\mathbin{\perp}(P\otimes t^\perp)\bigr)\otimes Q. \end{array}\] 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 \[A'=A\otimes F,\qquad P'=P\otimes\langle t\rangle,\qquad B'=(P')^\perp_{A'} =(B\otimes F)\mathbin{\perp}(P\otimes t^\perp).\] 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 \[A\otimes F=(B\otimes F)\mathbin{\perp}(P\otimes F),\] 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 \[P\otimes\langle t_X\rangle \ \subset\ P\otimes t_Y^\perp \ \subset\ P\otimes F.\] 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 \[Y\text{-side},\quad Y\text{-center},\quad X\text{-side},\quad X\text{-center},\quad Y\text{-center},\quad Y\text{-side},\quad X\text{-center},\quad X\text{-side}.\] 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.

\[ \begin{array}{c|l|l} \text{time}&\text{physical wires}&\text{common gate label}\\ \hline 0&Y,\ \text{its side wires}&B\otimes\langle t\rangle\\ 1&Y,\ \text{all center wires}&B\otimes F\\ 2&X,\ \text{its side wires}&(B\otimes F)\mathbin{\perp}(P\otimes\langle t\rangle)\\ 3&X,\ \text{all center wires}&A\otimes F\\ 4&Y,\ \text{all center wires}&B\otimes F\\ 5&Y,\ \text{its side wires}&(B\otimes F)\mathbin{\perp}(P\otimes t^\perp)\\ 6&X,\ \text{all center wires}&A\otimes F\\ 7&X,\ \text{its side wires}&A\otimes F. \end{array} \tag{6}\] 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 \[\frac{L_{\rm b}}{N}=\frac{100}{539},\qquad \frac{L_{\rm c}}{N}=\frac{101}{539}.\] In particular, both are smaller than \(1/2\). We now turn the dimension count into the two interfaces needed later.

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.

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 \[q_U(x)=\operatorname{wt}(P_Ux)\pmod4.\] The projection here is over \(\mathbb F_2\). Let \[H=\frac1{\sqrt2}\begin{pmatrix}1&1\\1&-1\end{pmatrix},\quad C=H\begin{pmatrix}1&0\\0&i\end{pmatrix}H =aI+bX,\qquad a=\frac{1+i}{2},\quad b=\frac{1-i}{2},\] where \(X\) interchanges the two coordinates. On functions of \(x\in\mathbb F_2^m\) define \[\mathcal C_U=H^{\otimes m}\operatorname{diag}_x(i^{q_U(x)})H^{\otimes m}.\] 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.

Explicit rational bounds for the two recurrences

The relative rank deficits are \[\eta_{\rm b}=\frac{W_{\rm b}m-s_{\rm b}}{W_{\rm b}m} =\frac{339}{22587335000000},\qquad \eta_{\rm c}=\frac{W_{\rm c}m-s_{\rm c}}{W_{\rm c}m} =\frac{73}{19906842167500}.\] For both recurrences we may use the fixed rational exponent \[ \tau=\sigma=1-2^{-50}. \tag{10}\] 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 \[m^{1-2^{-50}} =m\exp(-2^{-50}\log m) >m(1-20\cdot2^{-50})>m(1-\eta)\] for either \(\eta=\eta_{\rm b}\) or \(\eta=\eta_{\rm c}\). Consequently \[\frac{s_{\rm b}}{W_{\rm b}}<m^\tau, \qquad \frac{s_{\rm c}}{W_{\rm c}}<m^\sigma, \qquad 0<\tau=\sigma<1.\] Here \(\sigma\) is a rational upper bound for the logarithmic exponent of the complex recurrence.

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 (Grigor’ev 1982, Appendix A.2, Proposition 14(c)). We give its elementary elimination proof in this order.

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 \[A_e=M_{\mathrm{head}(e)}-M_{\mathrm{tail}(e)}, \qquad s=\sum_e\operatorname{rank}_{\mathbb Q}A_e.\] The required conditions are \[ M_{\mathrm{out}(\rho(w))}-M_{\mathrm{in}(w)}=I_m \quad(1\le w\le W),\qquad s<Wm. \tag{11}\] 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.

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 \[(h_0,H,g,D)\longmapsto(h_0,D,g,H) \longmapsto(D,g,h_0,H).\] 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 \[(H,g,d_0,D)\longmapsto(d_0,H,g,D) \longmapsto(d_0,D,g,H).\] 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.

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 \[H_0=\frac12\begin{pmatrix}1&1\\1&-1\end{pmatrix},\qquad H=\sqrt2 H_0,\qquad \mathsf S=\operatorname{diag}(1,i),\qquad Z=\operatorname{diag}(1,-1).\] For the bit-flip matrix \(X=\bigl(\begin{smallmatrix}0&1\\1&0\end{smallmatrix}\bigr)\), define \[ C=H\mathsf SH=aI+bX,\qquad a=\frac{1+i}{2},\quad b=\frac{1-i}{2}. \tag{13}\] Direct multiplication gives \[ H_0=b\mathsf SC\mathsf S,\qquad C^2=X,\qquad C^4=I,\qquad C^{-1}=bI+aX=-iZCZ. \tag{14}\] 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 \[\mathcal C_{U_+}\mathcal C_{U_-}^{-1} =\prod_{j=1}^{\rho_{\rm edge}}(aI+bX_{v_j})^{\epsilon_j}, \qquad \epsilon_j\in\{1,-1\}.\] 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.

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.

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 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 \[E=64(W+m+1)^3,\qquad B=s+E, \qquad \nu=\min\{j\ge1:m^j\ge B\},\] and set \[C_1=\max\{20,\nu+3\},\qquad C_0=128mB^2, \qquad \Delta=\lceil C_0d^{C_1}\rceil.\] For the stated \(h=100\) construction, the complex network has \[\begin{split} m&=1{,}000{,}000,\\ W&=1{,}873{,}807{,}244{,}643{,}542{,}670{,}000,\\ s&=1{,}873{,}807{,}244{,}636{,}671{,}267{,}308{,}000{,}000. \end{split}\] 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 \[A(e)\le sA(e/m)+E \quad\hbox{at an internal node},\qquad A(e)\le8e \quad\hbox{at a leaf}.\] 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 \[B^{k+1}\le B^{\log_m d+2} =B^2d^{\log_m B}\le B^2d^\nu.\] Thus \[A(e)\le 8dB^{k+1}\le8B^2d^{\nu+1}.\] 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: \[\begin{aligned} 8mB^2d^{\nu+2}+8d+t+9 &\le 8mB^2d^{\nu+2}+18d\\ &\le 26mB^2d^{\nu+2}\\ &<128mB^2d^{\nu+3}\le C_0d^{C_1}. \end{aligned}\] 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.

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 \[(F_t^-z)_j=\frac1t\sum_{k=0}^{t-1}\omega_t^{-jk}z_k, \qquad (F_t^+z)_j=\frac1t\sum_{k=0}^{t-1}\omega_t^{jk}z_k.\] 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 \[H_0=\frac12\begin{pmatrix}1&1\\1&-1\end{pmatrix}.\] 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\}\), \[ (F_t^-z)_{2h+b} =\frac1{t/2}\sum_{k=0}^{t/2-1} (\omega_t^2)^{-hk}\omega_t^{-bk} \frac{z_k+(-1)^b z_{k+t/2}}2. \tag{27}\] 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 \[k=\sum_{v=0}^{h-1}2^v x_v,\qquad 0\leq k<u/2.\] 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 \[ B_tF_t^-,\qquad F_t^+B_t^{-1}, \tag{28}\] respectively.

For \(D\) axes of lengths \(t_1,\ldots,t_D\), put \[M=\prod_{i=1}^D t_i,\qquad F_{\boldsymbol t}^\pm=\bigotimes_{i=1}^D F_{t_i}^\pm, \qquad B=\bigotimes_{i=1}^D B_{t_i}.\] 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 \[ y^{E},\qquad E=-\sum_{i\ \mathrm{participating}} \frac{2r}{u_i}\,b_i k_i\pmod{2r}. \tag{29}\] 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.

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.

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 \[\mathcal M_r(f,g)=\frac{f*g}{rM}=\frac{f*g}{T}.\] 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: \[\|\mathcal M_r(f,g)\|_\infty \leq\|f\|_\infty\|g\|_\infty.\] If the other operand is a disk array, replacing one operand changes the output by at most the norm distance between its two values.

The one-dimensional interface

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.

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.

A fixed rational choice

The two motif exponents were fixed in (10). Take the following fixed rational numbers: \[ \begin{gathered} \tau=\sigma=1-2^{-50},\qquad \beta=\frac12,\qquad \delta=\frac1{16},\\ \lambda=1-2^{-52},\qquad c=2^{-56},\qquad \lambda'=1-2^{-54},\\ C_1=20,\qquad \epsilon=2^{-75},\qquad \kappa=2^{-182}. \end{gathered} \tag{38}\] These constants are deliberately conservative; no optimization is claimed.

The layer conditions hold with room to spare: \[\max\{\tau,\sigma\}<\lambda<1,\qquad \tau(1+c/\beta)<1-31\cdot2^{-55}<\lambda,\] and \[\sigma+\beta(1-\sigma)=1-2^{-51}<\lambda<\lambda'<1.\] The strict final inequality absorbs the logarithmic number of groups in a base-\(m\) partition. The dimension exponent also satisfies \[\begin{gathered} \epsilon<\frac1{12},\quad \epsilon C_1<1,\quad \epsilon(2-\tau)<1-\tau,\quad \frac34+\delta+\frac32\epsilon<1,\\ \epsilon(1+c)<1,\qquad \epsilon+\delta<1. \end{gathered}\]

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.

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 (Agarwal and Cooley 1977). 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 \[P_i=\prod_{j<i}s_j,\qquad \mu_i=P_i^{-1}\pmod{s_i}.\] Here \(P_i\) denotes an integer prefix product, not a Fourier-coordinate permutation. The assignment \[\Phi:\mathbb C[x]/(x^S-1)\longrightarrow \mathbb C[x_1,\ldots,x_d]/(x_1^{s_1}-1,\ldots,x_d^{s_d}-1), \qquad x\longmapsto\prod_i x_i^{\mu_i}\] 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 \[k=\sum_{j=1}^d a_jP_j,\qquad 0\le a_j<s_j.\] The coordinate of \(\Phi(x^k)\) in axis \(i\) is \[ b_i=a_i+\mu_i\sum_{j<i}a_jP_j\pmod{s_i}. \tag{42}\] 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 \[\mathcal P=\prod_i[0,t_i),\qquad \mathcal V=\prod_i[0,s_i)\] 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 \[O\bigl(Tp\,d^2(1+\ell^\tau)\bigr).\] 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\).

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 \[\Theta:\mathbb C[x_d]/(x_d^r-1)\longrightarrow\mathcal R, \qquad x_d\longmapsto\zeta y\] 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 \[\zeta^{k+l}(-1)^h=\zeta^j(\zeta^r)^h(-1)^h=\zeta^j.\] 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 \[(\iota f)_{\boldsymbol j} =\sum_{k=0}^{r-1}f_{\boldsymbol j,k}y^k, \qquad \boldsymbol j\in\prod_{i<d}\mathbb Z/t_i\mathbb Z.\] 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 \[\mathcal M_{\mathbb C}(f,g) :=D_-\iota^{-1} \mathcal M_r(\iota D_+f,\iota D_+g) =\frac{f*g}{T}.\] 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.

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 \[-\sum_i s_i j_i^2(r/t_i)\pmod{2r}.\] 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 \[\begin{aligned} x_0&=Q_p(\overline{\widetilde a}\mathbin{\cdot}x),\\ y_0&=\widetilde{\mathcal M}_{\mathbb C}(\widetilde a,x_0),\\ \widetilde{\mathscr C}(x) &=Q_p(\overline{\widetilde a}\mathbin{\cdot}y_0). \end{aligned}\] 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.

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. 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. Compute \(\widetilde{\mathscr F}(u)\) and \(\widetilde{\mathscr F}(v)\).

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. 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. 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. 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.

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\).

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: \[ \begin{split} g_1&=1-\epsilon(1+c),\qquad g_2=\epsilon c(1-\tau),\qquad g_3=\epsilon(1-\lambda'),\\ g_4&=1-\tau-\epsilon(2-\tau),\qquad g_5=1/4-\delta-3\epsilon/2,\\ g_6&=1-\delta-\epsilon,\qquad g_7=\epsilon. \end{split} \tag{45}\] Direct substitution gives \[g_1>\frac12,\quad g_2=2^{-181},\quad g_3=2^{-129},\quad g_4>2^{-51},\quad g_5>\frac18,\quad g_6>\frac12,\quad g_7=2^{-75}.\] 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 \[O(p/K)\,K^\tau=O(pK^{\tau-1}).\] 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 \[\sum_{i=1}^d\frac{T}{t_i} O(t_i p^{3/2+\delta}\alpha) =O(dTp^{3/2+\delta}\alpha).\] 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.