大卫·哈维评测近期关于 o(n log n) 整数乘法算法的论文

On the recently proposed o(n log n)integer multiplication algorithm

David Harvey, UNSW Sydney 20th October 2026

Abstract

I offer some thoughts about the recently proposed integer multiplication algorithm whose complexity is claimed to be o(n log n) bit operations.

1 Executive summary

I have not checked most of the details, but I think that the claimed o(n log n) bound is plausible.

2 Declaration

This note was written without AI assistance, except for (a) typesetting (i.e., TeX wizardry), (b) proofreading, and most importantly (c) extensive discussions to help me understand the new algo-rithms. The main system used throughout was Claude Fable 5.1.

3 Introduction and background

A certain American technology company (henceforth “the Company”) recently dumped a large pile of manuscripts on GitHub. This “math(s) dump” contains many remarkable claims and results. I have struggled to find the words to convey the stupendousness of these discoveries to my non-mathematical associates. One paper in particular caught my attention [TC26a]. It proposes a new algorithm for multiplying

n-bit integers, claiming that it runs in o(n log n) bit operations, in either the multitape Turing model or the Boolean circuit model. The small “ o” is not a typo. More precisely, the paper claims a complexity bound of the form O(n(log n)1−ε) for some small positive constant ε.This bound is of course utterly ridiculous. Previously the best bound was O(n log n) (big “ O”), and this was widely believed to be optimal. The dump contains many other ludicrous results, but most of them were at least supposed to be true . This one was not. Nevertheless, this is a serious Company with a track record of doing serious things, so we should make an effort to take the claim seriously. 1I would like to congratulate the Company, or more precisely its employees, for creating a Model with these astonishing capabilities. My understanding is that none of these people actually know how the Model works. Nor do they, for the most part, understand the content of the papers that they just publicly released. I would also like to congratulate the unnamed Model itself. Both for its brilliant ideas, and also for its perseverance against unfathomable odds. I assume that by now the Company’s employees are aware of the pain and profound sense of loss that this event has caused in my professional community. I believe that they are decent and hard-working people, whose efforts are sincerely aimed at making the world a better place. I hope and pray that they succeed. Some of my colleagues believe that this is a nefarious Company, whose only goal is profit, if necessary at our community’s expense. Perhaps there is some truth to this, but I think it is missing the larger point. The fact is that there now exists an automated system that can perform mathematical research — at least the part that involves “solving hard problems” — far better than any human. Now that such a system has been demonstrated, it is clear that results like these were bound to start appearing soon anyway. I do not blame the Company for getting there first. Who got there first is an irrelevant detail, a historical footnote. I am confident that this is not a one-off. At some point soon, the Model, or one of equivalent capability, possibly developed by a different Company, will be released publicly. At first it will be expensive, but then it will become cheap, and perhaps even free. Its capabilities will continue to improve. More mathematical results will appear, with increasing frequency, produced not only by the Company and its peers, but by ordinary people with access to these Models. Everyone in my community must grapple with the consequences. For instance, I doubt that it is possible for humans to be taken seriously as authors of mathematics research papers for much longer, for the same reason that we no longer accept written assignments from our undergraduates. It is possible that we have already reached this point; we will soon find out. This is not the place to dig into my views on the wider impacts of AI, but I can recommend the excellent blog by Zvi Mowshowitz ( https://thezvi.substack.com/ ), whose perspective on AI is broadly similar to my own. I learn a huge amount from his voluminous and highly addictive writing. Zvi, if you are reading this, many thanks for your efforts, they are greatly appreciated. Returning now to the o(n log n) integer multiplication result. The new work builds directly on top of the O(n log n) algorithm described by myself and Joris van der Hoeven [HvdH21] (henceforth “the 2021 algorithm”). Much commentary on the interweb has focused on how insanely small the (log n)ε improvement is over the previous O(n log n) record, harping especially on the tiny numerical value of ε. In the initial writeup ε is 2 −182 , meaning that the exponent of log n is roughly 0.99999999999999999999999999999999999999999999999999999984 .

But from where I stand, this is actually a gobsmackingly enormous improvement. To put the gap in perspective, recall that when F¨ urer established his O(n log n · Klog ∗ n) bound in 2007, the improvement that he had achieved over the previous state of the art, namely the epochal Sch¨ onhage– Strassen O(n log n log log n) bound from 1971, was smaller than log log n. And the subsequent 2improvement made by my work with van der Hoeven (and several other papers along the way) was much, much, much smaller than this. The claimed (log n)ε improvement in the new work is

exponentially larger than any of these. The actual value of ε is irrelevant. This is a huge gap. English has many synonyms for “big” but not enough for what I am trying to say. I suspect that very few people (i.e., humans) worldwide understand all the details of the 2021 algorithm, but I happen to be one of them. This puts me in the rare position of having the background to be able to comment immediately on the new work.

4 Executive summary (longer version)

The bottom line is that I think the claimed result is plausible. I have not actually read the paper. In fact, there are three closely related papers [TC26a], [TC26b], [TC26c], and I have not read any of them. I tried, but soon gave up. The math(s) dump does not appear to be strong on exposition. Instead, I asked Fable to explain [TC26a] to me. Fable seems to understand the argument and thinks that it is correct. I expect that thousands of frontier LLM instances around the world have come to the same conclusion over the last two weeks. If not, surely we would have heard by now. At first I did not understand Fable’s explanations. But I asked questions, and then asked more questions, and drilled down towards the core of the problem. Gradually I started to understand parts of the argument. I have seen enough now to arrive at the tentative belief that the main result is possibly correct (notice that I choose my words very carefully), and that the bound O(n(log n)α)has the right shape. I stress that I have not checked the proof in any serious level of detail, and there are many parts that I do not yet understand properly. In the sections below I will try to explain what I have learned so far, and why I lean towards thinking it might be correct. Incidentally, when I say that the math(s) dump “does not appear to be strong on exposition”, all I mean is that I, as an human expert in this field, have difficulty understanding it (specifically [TC26a]). I have heard similar comments from other human experts regarding other papers in the dump. But as far as Fable seems to be concerned, the exposition in [TC26a] is perfectly clear. So perhaps the problem is not with the exposition, but with us humans. Before the advent of AI, it was already the case that different styles or levels of exposition were appropriate for different audiences. Perhaps we need to ask the LLMs to “explain it to me like I’m human” (ELIH). There is one other aspect of the situation that I would like to emphasise. After two weeks of very hard work, I have managed to understand only a small fraction of [TC26a]. When I ask Fable about the paper, it seems to read and understand the whole thing in just a few minutes, even starting from a fresh context. I know that some readers may be uncomfortable with my anthropomorphic terminology, but I believe that “understand” is accurate, in the sense that the LLM seems to be able to correctly and lucidly answer all kinds of questions about the paper, including about some very fiddly technical details. That is exactly how I would tell that a human colleague has “understood” something. This chasm between human and LLM “speed of thought” is unsettling to say the least. It would be good to have a formalised proof of the o(n log n) result, for example in Lean. As far as I can tell, one does not yet exist. Such a formalisation seems to be achievable by current technology, but it is still not easy. If the result is correct, I expect a formalisation to be available fairly soon. My confidence in the result would be increased greatly if such a formalisation existed, assuming 3that the initial definitions and statement of the main theorem had been checked by a (human?) Lean expert. I apologise that in this note I have not maintained my usual level of care and scholarship. Everything is happening so fast. In particular I apologise to those people whose work is used in this note without adequate attribution, due to my ignorance and/or negligence. Maybe one day I will make a more careful attempt, though I suspect that an LLM will beat me to it. Before getting to the explanations, I want to mention some consequences of the main result, as-suming that it is correct:

ˆ There is a now a quantity that I would suggest calling the “exponent of integer multiplication”. It is the infimum of real numbers α ⩾ 0 such that there exists an integer multiplication running in time O(n(log n)α). We previously believed that α = 1, but the Model’s claim is that α < 1. Already some people have found improvements in the value of α. At the time of writing the best value appears to be α ≈ 0.9983, but this number will surely be out of date by the time you are reading this. In fact, there might be two different versions of α: a “Boolean circuit” version that only takes into account the cost of arithmetic, and a “Turing machine” version that also includes the cost of data rearrangement. It is not at all clear to me that these exponents should take the same value. Some people have asked me, how low could α go? Of course I have no idea, but I want to men-tion one obstruction. In [HvdH21], the Gaussian resampling step contributes altogether time

O(n(log n)3/4). As sketched in [HvdH21, §4.4.2], this step can be improved to O(n(log n)1/2). So without further improvements here, or a completely different approach that does not rely on Gaussian resampling, it does not seem possible to do better than α = 1 /2. One method that does not depend on Gaussian resampling is the variant in [HvdH22], but this unfortunately requires some unproved number-theoretic hypotheses. There has been so much progress in analytic number theory in recent weeks that it might be worth taking another look. Either way, we are still very far away from the 1 /2 limit, so I would advise against thinking about it too hard yet.

ˆ Together with the Boolean circuit lower bound result of [AFKL19], the new result disproves the Li–Li network coding conjecture [LL04].

ˆ The new result does not directly imply new complexity bounds for polynomial multiplication over a fixed finite field, say multiplication in F2[x]. The state of the art for this problem is an algorithm that multiplies polynomials of degree n in time O(n log n 4log ∗ n) [HvdH19]. I do not immediately see how to adapt the new method to that algorithm. However, the method of [HvdH22] also works for polynomials over finite fields, and achieves the better bound O(n log n) under certain number-theoretic hypotheses. The structure of that algorithm suggests that the new techniques are applicable, and that they would lower its complexity to

O(n(log n)α) (still conditional however).

ˆ My recent paper on fast enumeration of primes [Har26] relies heavily on fast multiplication in F2[x]. Maybe the results of that paper can be improved by plugging in the (conditional) multiplication algorithm just mentioned. I had a very quick look, and it seems that the 4bottleneck becomes the computation of the “restricted products”. So I am not sure, but it is certainly tempting to investigate. Some care is needed, because the current framework in the paper more or less assumes that the cost of multiplication is Ω( n log n); one would have to audit the entire complexity analysis from scratch without this assumption.

ˆ When the O(n log n) integer multiplication result was first announced in March 2019, Fredrik Johansson wrote on Reddit:

The result is of extreme practical importance. Not for actually multiplying integers (as usual with these algorithms it’s probably not faster than existing algorithms for integers that can be stored in the observable universe), but for writing papers. It has always been a hassle to write down the complexity of integer multiplication or algorithms based on integer multiplication by introducing soft-O notation, little-o exponents, epsilons greater than 0, or iterated logarithms. From now on I can just write O(n log n) in my papers and be done with it!

Alas, after seven years’ sojourn in this typographical Utopia, we are now condemned to wander forever in a wasteland of superscripts and curly braces. Or, if you prefer, you can continue to write O(n log n). No harm will come to you.

5 The 2021 algorithm

To understand the new integer multiplication algorithm, one must first understand essentially everything in the 2021 algorithm. It is impossible for me to explain all of this here. I will give instead a brief summary of the points that are relevant for understanding how the new algorithm changes the analysis. Readers interested in the details should consult [HvdH21], which by all accounts is a superbly written paper. We start with two n-bit integers. After a series of reductions that I will not attempt to explain (Agarwal–Cooley CRT rearrangement, Gaussian resampling, etc), we reduce to the problem of computing a Nussbaumer-style DFT (discrete Fourier transform) of a certain ( d − 1)-dimensional array of size r × · · · × r, whose entries lie in the ring R = C[y]/(yr + 1). Here d is a fixed dimension parameter (independent of n) and r is a power of two, of size roughly n1/d . The coefficient ring R

contains a principal 2 r-th root of unity, namely y, and this is the “synthetic” root of unity used to define the DFTs. Thus, the DFT may be computed by an FFT (fast Fourier transform) algorithm, which amounts to performing a series of “butterflies”, i.e., maps ( a, b ) 7 → (a + b, a − b) operating on pairs of complex numbers, together with multiplications by powers of y, which are just shifts and negations. The complex numbers themselves are represented as fixed point approximations with

p = O(log n) bits of precision. The number of butterflies (over C) is easily calculated to be 12 rd log 2(rd−1). One can already see the O(n log n) shape of the final complexity bound in this butterfly count. These reductions lead to a recursive bound for M(n) (the cost of multiplying n-bit integers) of the form

M(n) < Kn n′ M(n′) + O(n log n), n′ = n 1

d+o(1)

.

5The key point is that the “expansion factor” K is an absolute constant, independent of d. In the paper we show that one may take K = 1728; this arises from analysing the zero-padding and encoding steps that occur during the various reductions. (In §5.4 we sketch various techniques that would reduce it to K = 8.) The first ( Kn/n ′) M(n′) term arises from the cost of the recursive multiplications in the Nussbaumer convolution. Provided we choose d > K , for example d = 1729, then unrolling the recursion leads to a geometric progression and keeps the cost of the subsequent recursion levels under control. For simplicity, in [HvdH21] we keep d constant, since all we need is d > K . In the new algorithm one actually needs to take d increasing with n. When d grows with n, one must be more careful with certain estimates, but it introduces no essential difficulties. (In fact, in the initial stages of writing [HvdH21], we had d increasing with n, and only later in the project did we decide to fix d.) The crucial question now is where the O(n log n) term comes from. In the Boolean circuit model, it arises entirely from the cost of the butterfly arithmetic; everything else is negligible. In the Turing model, there are additional data rearrangement costs, such as transposing arrays, which also contribute O(n log n). If you had asked me two weeks ago for my subjective estimate of the probability that either of these O(n log n) contributions could be improved, I think I would have guessed roughly 0 .01% for the data rearrangement steps, and precisely zero for the arithmetic steps. In the past, I have tried, perhaps half-heartedly, to invent faster matrix transposition algorithms. But it has never, ever, even remotely occurred to me to launch an attack on the butterfly arithmetic. Since I first learned about FFTs several decades ago, this barrier has always seemed to be woven into the fabric of the universe, a hard physical limit like the speed of light. But the Model now asserts that we can in fact go faster than light. In the following sections I will try to explain its reasoning. Once you believe that the butterfly arithmetic can be improved, then the rearrangement side becomes less surprising. We already know that there is a close relationship between these two problems [HvdH25]. Therefore, in my study of the new algorithm so far, I have completely ignored the data rearrangement part of the problem, and focused entirely on the arithmetic side, which I think is more fundamental. If you like, I am only thinking about the Boolean circuit model, and ignoring the Turing situation for the moment. This will be the context that I work in for the rest of this note.

6 Linear algebra and amplification

I am not really an algebra person. Tensor products are hard. Deep breath and repeat after me: “growth mindset”. Let A be a 2 × 2 complex matrix. It acts in a natural way on C2. For example, A could be the butterfly map ( 1 11 −1 ). Now let m be a positive integer and consider ( C2)⊗m. This is the m-th tensor power of C2. You can think of an element of ( C2)⊗m ∼= C2m

as an m-dimensional array of complex numbers of size 62 × · · · × 2 ( m factors). If A1, . . . , A m are 2 × 2 matrices, then A1 ⊗ · · · ⊗ Am acts on ( C2)⊗m in a natural way:

A1 ⊗ · · · ⊗ Am = ( A1 ⊗ I ⊗ · · · ⊗ I)( I ⊗ A2 ⊗ · · · ⊗ I) · · · (I ⊗ I ⊗ · · · ⊗ Am).

This decomposition is just algebra, but it also expresses an algorithm for applying A1 ⊗ · · · ⊗ Am:first act on slices in the m-th direction by Am, then act on slices in the ( m − 1)-th direction by

Am−1, and so on. And because the slice operators commute, you can do these m steps in any order. Each Ai is applied 2 m−1 times, so the total number of 2 × 2 matrix applications is m · 2m−1.Let’s now focus on the case A1 = · · · = Am = A, so we are studying the tensor power A⊗m of a single A. Applying A⊗m to a vector in ( C2)⊗m using the above algorithm requires m · 2m−1

applications of A. For example, if A is the butterfly matrix, then this is just the standard number of butterflies needed to compute the DFT of a 2 × · · · × 2 array. Now it is conceivable that for certain A, one might be able to evaluate A⊗m using fewer than

m · 2m−1 applications of A. This is trivially true for example if A = 0 or A = I! Suppose however that A is not so trivial, and that someone magically hands us an algorithm for evaluating A⊗m for some m, that uses only t applications of A, where t is strictly smaller than m · 2m−1. How can we use this algorithm to speed up evaluations of A⊗n, where n is much larger than m?Here is a fairly standard method. Suppose for simplicity that n = mℓ for some integer ℓ ⩾ 1. Then

A⊗n = A⊗mℓ ∼= A⊗m ⊗ · · · ⊗ A⊗m (ℓ factors). Essentially we have split up the original n axes of the array into ℓ groups of m axes. So what we could do is apply our fast algorithm for A⊗m

along the first group of axes, then along the second group, and so on. This reduces the number of applications of A from the original n · 2n−1 to ( t/ (m · 2m−1)) · n · 2n−1. This saving is not very interesting, because it is by a factor of only t/ (m · 2m−1), regardless of the size of n.What we really want is a saving that grows with n. This can be achieved by the following method, which Fable refers to as “amplification”. (The paper [TC26a] does not use this term. It is used in the companion papers. Fable points to prior art in Alman (2021); Alman, Guan and Padaki (2023); Alman and Li (2025). I have not followed up these references. The amplification trick is certainly new to me. See also comments about Strassen’s matrix multiplication method below.) For amplification to work, we have to assume that our fast algorithm for A⊗m has some additional structure. Let us denote by ( F2)m the space of “addresses” for ( C2)⊗m ∼= C2m

, i.e., the vector

u = ( u1, . . . , u m) ∈ Fm

2

refers to the ( u1, . . . , u m)-th coordinate of ( C2)⊗m. We consider only algorithms given as sequences of the following two types of “instructions”:

ˆ Slice call . For each i = 1 , . . . , m , the operator Ai := I ⊗ · · · ⊗ A ⊗ · · · ⊗ I, i.e., Ai acts on (C2)⊗m by applying A along the i-th direction. This amounts to 2 m−1 invocations of A.

ˆ Relabelling. For each matrix M ∈ GL m(F2), the operator PM that acts on ( C2)⊗m by permuting the entries according to x 7 → M x (for x ∈ Fm

2

). For example, if M is a permutation matrix that swaps two coordinates of Fm

2

, then PM just switches two axes of the array. But M could also involve more complicated transformations 7that “mix up” the axes (i.e., bits of the addresses) in interesting (linear) ways. (Of course, this setup is not general enough to allow arbitrary permutations of the 2 m coordinates of (C2)⊗m.) Note that in the Boolean circuit model, these sorts of relabellings are entirely free. In the Turing model we would have to account for their cost, but I am ignoring this here. An example of such a “program” is the obvious one A⊗m = A1A2 · · · Am that we already discussed earlier. This performs exactly m slice calls, with no relabelling steps. Suppose now that we are magically given a program for computing A⊗m that performs only g slice calls, where g is strictly less than m. Without loss of generality it has the form PM0 Ai1 PM1 Ai2 · · · Aig PMg .(Of course we could get away with just A1 if we modified the M ’s appropriately. But it is convenient to allow ourselves the extra flexibility of the Ai’s.) Let us try to use this program to speed up the computation of A⊗n, where as before n = mℓ for some ℓ ⩾ 1. First observe that

A⊗n = ( A⊗m)⊗ℓ = A⊗m ⊗ · · · ⊗ A⊗m.

In terms of arrays, the formula is splitting up the n axes into ℓ groups of m axes. Now we substitute in the program:

A⊗n = ( PM0 Ai1 PM1 Ai2 · · · Aig PMg )⊗ℓ.

Keep in mind that each of the ℓ copies of (say) Ai1 is actually operating on a different group of m

axes. We can bring the tensor product inside the composition:

A⊗n = ( PM0 )⊗ℓ(Ai1 )⊗ℓ · · · (Aig )⊗ℓ(PMg )⊗ℓ.

I don’t know about you, but at this point my head is spinning. Let’s try to decode what is happening. Consider first one of the factors ( PMj )⊗ℓ. This is applying the same relabelling map PMj (which only acts on m axes) to each of the ℓ groups of m axes. You could do them one group at a time, or all at once; in the Boolean circuit model it makes no difference, it is free. Now consider one of the factors ( Aij )⊗ℓ. This is really ( I ⊗ · · · ⊗ A ⊗ · · · ⊗ I)⊗ℓ, where we are applying A in the ij -th component. In other words, this map is applying A⊗ℓ to a subset of ℓ of the original n axes; it chooses one axis (the ij -th one) out of each group of m axes, and only operates on those. Let’s compare the total work to what would have happened if we had instead used the standard program A1 · · · Am for A⊗m. In effect, we have reduced the number of calls to A⊗ℓ by a factor of

g/m . So if we don’t use any special techniques to evaluate A⊗ℓ, then we save a factor of g/m in the total number of applications of A.But of course, if ℓ is itself a multiple of m, we can (and should!) apply this strategy to A⊗ℓ itself to save another factor of g/m . And then of course we should do it recursively: if we take n = mk

at the beginning, then overall we save a factor of ( g/m )k. This is the amplification: the additional 8structure lets us save further factors of g/m as n grows. The total number of applications of A

becomes (g/m )k · n · 2n−1 = exp

 log n

log m log( g/m )



· n · 2n−1

= nα · 2n−1

where α = log g/ log m. You can see here the N (log N )α shape of the result (where N = 2 n). Notice that if g < m , then α is strictly less than 1. (The amplification trick is reminiscent of Strassen’s matrix multiplication algorithm. But notice that in Strassen’s method, k recursion levels gets you from dimension 2 to dimension 2 k, whereas here k recursion levels gets you from dimension 2 m to dimension 2 mk

, i.e., doubly instead of singly exponential in k. Fable explained the difference to me as follows: in Strassen’s method, you are “tensoring a fixed gadget with the problem”, whereas here we are “substituting the problem into the gadget as a tensor power”. So I think they are algebraically not the same thing.) Applying the preceding construction to the butterfly map A = ( 1 11 −1 ), we see that if we can find a program for A⊗m with g < m for even one single value of m, then we can compute DFTs on (C2)⊗n for large n in essentially O(N (log N )α) amount of work ( N = 2 n), with α < 1. The picture, as I have sketched it so far, is missing two main ingredients. First, as mentioned earlier, in the 2021 algorithm we don’t want to compute DFTs of size 2 × · · · × 2 over C, but rather DFTs of size r × · · · × r over R = C[y]/(yr + 1). I am pretty sure that the method will transfer to this situation. The latter problem can certainly be decomposed into butterflies over C of the right shape. I’m not sure exactly what happens to the twiddle factor multiplications (i.e., synthetic roots of unity), but these should be free in the Boolean circuit model (except possibly for chasing some signs around), and I doubt that they interfere in an essential way with the key amplification idea, which is really about how to deal with high dimensional tensor products. You can also see why we need the dimension parameter d to increase with n; we need room to allow n = mk with growing

k. All of this sounds plausible to me but I haven’t worked through it carefully yet. The second thing we are missing is that the Model has not in fact found a way of obtaining a program with g < m . (Moreover, Fable tells me this is impossible, but I do not follow the argument yet.) Instead, what it has found, as I will discuss in the next section, is a way of doing this on average when evaluating A⊗m on sufficiently many independent input arrays. More precisely, taking A to be the butterfly map, and for an integer W ⩾ 1, consider the problem of evaluating I ⊗ A⊗m on CW ⊗ (C2)⊗m. An element of CW ⊗ (C2)⊗m is just a list of W arrays of size 2 m, and evaluating I ⊗ A⊗m means applying A⊗m to each of the W arrays. The standard number of slice calls needed to evaluate I ⊗ A⊗m is W m , i.e., just transform each of the arrays separately. In the next section I will discuss how to construct a gadget that gets away with slightly less (for suitably chosen m and W ). The gadget also includes some other linear operations that “mix up” the entries of the W arrays. Note that achieving even one fewer slice call, for a fixed

m and W , is enough to make the whole argument go through. (Well, this is not the whole story, because we also have to track the number of “mix up” operations, which I have not thought about carefully yet. The reason I am not too worried is because of the analogy with Strassen’s matrix 9multiplication; the key issue is to reduce the number of recursive calls, and the other arithmetic (like the extra additions in Strassen’s algorithm) should hopefully just sort itself out.) Once we have such a gadget available, it is not hard to modify the calling framework to take advantage of it; after all we do have to perform lots of independent butterflies, and W is fixed. (Fable tells me that the choice of how you allocate these W coefficients can affect the exponent, but I haven’t looked into this yet.) The reader might think that it sounds crazy that one could speed up W independent evaluations of A⊗m without being able to speed up a single evaluation on its own. This was certainly my first reaction. But in fact, there is a compelling precedent. Consider the problem of multiplying an

n × n matrix by a single vector v ∈ Cn. Assuming the matrix has no special structure, this requires Ω( n2) operations. But if you want to multiply the matrix by n independent vectors, this is actually a matrix-matrix multiplication problem, and we learned from Strassen that we can do this faster than n3, i.e., o(n2) operations on average per vector. So I do not see any obvious obstruction here.

7 The gadget

Let B be the butterfly map ( 1 11 −1 ) acting on C2. The goal in this section is to explain the con-struction of positive integers m and W and a “program” that computes the action of I ⊗ B⊗m on

CW ⊗(C2)⊗m in fewer than W m slice calls. As mentioned above, there are also extra additions, but I will ignore these, I will focus only on the slice calls, which is the parameter driving the recursive gain.

7.1 Normalisations

First some normalisations. Let

H := 1

√2 B = 1

√2

1 11 −1



.

So H (for Hadamard) is just a scaling of B to be unitary. Next let X := ( 0 1 1 0 ) and

C := aI + bX =

a bb a



, a := 1 + i

2 , b := 1 − i

2 .

Then |a| = |b| = 1 /√2, and one checks that 1 + i

√2 · H = SCS where S :=

1 00 i



.

So we can think of C as just a version of H where we have changed basis by S on the domain and S−1 on the codomain (or maybe I have these backwards), and then thrown in a scaling by the norm-1 factor (1 + i)/√2. If we can compute C⊗m then we can easily recover the effect of B⊗m,we just have to pre- and post-multiply every coordinate by the correct powers of i, 1 + i and √2which is easy. So for the rest of the note I will forget about B and work with C.Observe that C is unitary. One checks that it is diagonalised by H, and has eigenvalues 1 and i:

C = HSH −1.

10 (Of course H−1 = H but I prefer to write H−1 to make it clear that we are using the same change of basis on both sides.) (I will just pause here to mention something I find very mysterious. The rest of the development depends crucially on the fact that the eigenvalues of C are 1 and i. But there is nothing “natural” about this; we introduced it artificially when we decided how to change bases to get from B to C.The only justification I see at the moment is that it makes everything work later on. But I have no idea why. (Update: Fable tells me this is forced by the requirements that C commutes with X,which is what makes all the Cv simultaneously diagonal in the H⊗m basis (see below), and that C

is the same as the butterfly up to free rescalings. I have not checked this.))

7.2 The Cv operators

As I said above, C is basically our butterfly expressed in terms of different bases on the domain and codomain. Now we want to graduate to C⊗m, and understand what applying C across a “slice” of (C2)⊗m looks like. The simplest case is where we just apply C on the i-th coordinate. But we will actually need something more general that this. As a notational matter, we will sometimes write elements of ( C2)⊗m as functions Fm

2

→ C, where

Fm

2

are the addresses as before. Let v ∈ Fm

2

, v̸ = 0. Think of v as a “direction”. For example if v = ei then this is just the i-th coordinate direction. Let Xv denote translation by v, i.e., it sends f ∈ (C2)⊗m to Xvf ∈ (C2)⊗m given by ( Xvf )( x) =

f (x + v). Since v + v = 0, you can also think of this operation as swapping x and x + v, for every

x ∈ Fm

2

. Then define

Cv := aI + bX v.

Thus Cv acts on ( C2)⊗m by applying C to the 2 m−1 pairs ( x, x + v). In other words it applies C

in “direction v”. For example, Cei means to apply C on all pairs in the i-th coordinate direction. In the previous section, a “slice call” was only allowed to be in a coordinate direction, i.e., of the form Cei . But we can easily reduce a general Cv to this case by inserting appropriate relabellings before and afterwards. So from now on, “slice call” is allowed to be the more general Cv. Also, one easily checks that ( Cv)−1 = CvXv, and since Xv is just a free rearrangement, we can also treat (Cv)−1 as a slice call. Although the actual algorithm works with the standard basis for ( C2)⊗m, there is a bunch of bookkeeping we need to do in a different basis, namely the basis determined by the m-dimensional Hadamard map H⊗m. Our next job is to work out what Xv and Cv look like in this basis. For Xv we have (H⊗m)−1 · Xv · H⊗m = diag(( −1) v·x),

where diag( f (x)) means the diagonal matrix with x-th entry f (x), for each x ∈ Fm

2

. Here v · x is the standard bilinear form on Fm

2

, i.e., v1x1 + · · · + vmxm ∈ F2. Sketch proof: first check it on Xei

by direct calculation, and then observe that translation by v amounts to translating by each ei for 11 which vi = 1. (The result also falls directly out of standard results about Fourier transforms on the finite group ( Z/2Z)m.) For Cv we just use the formula Cv = aI + bX v. This gives (H⊗m)−1 · Cv · H⊗m = diag( a + b(−1) v·x).

If v · x = 0 then the RHS is a + b = 1; if v · x = 1 then we get a − b = i. Therefore (H⊗m)−1 · Cv · H⊗m = diag( i[v·x]), (7.1) where for z ∈ F2 we define [ z] ∈ Z/4Z by the naive lift, i.e., 0 mod 2 gets mapped to 0 mod 4 and 1 mod 2 gets mapped to 1 mod 4. Quick numerical example as a sanity check. Take m = 3. Then Fm

2

= {000 , 001 , 010 , 011 , 100 , 101 , 110 , 111 }

and ( C2)⊗m ∼= C8. Take v = 011. I won’t bother to write out the 8 ×8 matrix for Cv in the standard basis, but to work it out in the H⊗m basis we first compute (v · 000 , v · 001 , . . . ) = (0 , 1, 1, 0, 0, 1, 1, 0) .

So the matrix for Cv is just diag(1 , i, i, 1, 1, i, i, 1). A really nice feature of this setup is that all of the Cv’s are diagonal in this basis. So all the Cv

commute with each other, and when we multiply them together, the diagonal entries just multiply. For example, since C⊗m = Ce1 · · · Cem , we get (H⊗m)−1 · C⊗m · H⊗m = diag( i[e1·x]) · · · diag( i[em·x]) = diag( iwt( x)) (7.2) where wt( x) denotes the weight of x, i.e., the number of entries equal to 1. In our m = 3 example, the matrix for C⊗m would be diag(1 , i, i, −1, i, −1, −1, −i). Having introduced the above machinery, we can now see roughly the shape of the algorithm we are aiming for. Suppose you give me an element of CW ⊗ (C2)⊗m, i.e., a list of W arrays. I want to compute C⊗m of each one. The standard algorithm is to perform Ce1 , . . . , C em on each array; this amounts to W m slice calls. But I have other options. I could apply various Cv’s to various arrays in the list, and I know the effect in the H⊗m basis via (7.1); I just need to keep track of how the exponents of i add up modulo 4. To mix things up a bit, I can also do things like add one array to another, which means that the various Cv are in effect going to be acting on several input arrays “in parallel”. The aim is to somehow get the right pattern of powers of i, i.e., to match those of

C⊗m, in each coordinate, at the very end, and with all the “mixups” having been undone.

7.3 Building up nondegenerate subspaces of Fm

2

To keep track of the effect of applying various Cv’s to ( C2)⊗m, we introduce the following framework. Let U ⊆ Fm

2

be a subspace. We write U ⊥ ⊆ Fm

2

for the complement with respect to the standard dot product, i.e., U ⊥ := {v ∈ Fm

2

: u · v = 0 for all u ∈ U }. The rank-nullity theorem shows that always dim( U ⊥) = m − dim U .12 Recall that U is called nondegenerate with respect to the dot product if U ∩ U ⊥ = 0. This condition is equivalent to saying that U + U ⊥ = Fm

2

. Note that U is nondegenerate if and only if

U ⊥ is nondegenerate. We will work with nondegenerate subspaces throughout. For example, when m = 3, the span of 100 is nondegenerate (its complement is the span of 010 and 001). But the span of 110 is degenerate (because 110 · 110 = 0, so 110 lies in both U and U ⊥). More generally, for any m the span of {e1, . . . , e j } is nondegenerate for any j.When U is nondegenerate, it makes sense to define a projection operator PU : Fm

2

→ U . We just write a given x ∈ Fm

2

uniquely as y + z where y ∈ U and z ∈ U ⊥, and then put PU (x) := y.Now for any nondegenerate subspace U , define qU : Fm

2

→ Z/4Z by

qU (x) := wt( PU x) mod 4 .

For example, q0 is just the zero function, and qFm

2

is the function that returns wt( x) (mod 4), i.e., exactly the exponents appearing in (7.2). Then define ΦU := H⊗m · diag( iqU (x)) · (H⊗m)−1,

so Φ 0 = I, and Φ Fm

2

= C⊗m by (7.2). If U and V are nondegenerate and mutually orthogonal, then it is easy to see that U + V is also nondegenerate, and that PU +V = PU + PV .We now claim that weight modulo 4 is additive on orthogonal vectors in Fm

2

. You can prove this by observing that wt( x + y) = wt( x) + wt( y) − 2| supp x ∩ supp y| (the correction term comes from the coordinates where both xi = 1 and yi = 1), and if x · y = 0 then | supp x ∩ supp y| is even. It follows that that qU +V = qU + qV , since

qU +V (x) = wt( PU x + PV x) mod 4 = qU (x) + qV (x),

using the fact that PU x ⊥ PV x. Finally this means that Φ U +V = Φ U ΦV .So now we can start to see how the game works: we begin with U = 0, i.e., Φ 0 = I, and we try to add a sequence of orthogonal subspaces to U , until we finally reach Φ Fm

2

= C⊗m. As we proceed, we can measure our “progress” towards evaluating C⊗m by the dimension of the intermediate spaces

U that we visit along the way. Let us now investigate what a “single step” looks like. Let U be nondegenerate, and consider a nonzero vector orthogonal to U , say v ∈ U ⊥. We also need to assume that ⟨v⟩ is itself nondegenerate, i.e., we need wt( v) to be odd . Under these conditions, U ′ := U + ⟨v⟩ is also nondegenerate, and ΦU ′ = Φ U Φ⟨v⟩.All that remains is to calculate Φ ⟨v⟩. Observe that (since v · v = 1)

q⟨v⟩(x) = wt(( v · x)v) mod 4 = [ v · x] wt( v) mod 4 .

13 Since we assumed wt( v) is odd, there are two cases. If wt( v) ≡ 1 (mod 4), then we get q⟨v⟩(x) = [v · x] mod 4, and comparing with (7.1) we find that Φ ⟨v⟩ = Cv. If wt( v) ≡ 3 (mod 4), then

q⟨v⟩(x) = −[v · x] mod 4 and so Φ ⟨v⟩ = ( Cv)−1. In summary, Φ⟨v⟩ =

(

Cv, wt( v) ≡ 1 (mod 4) ,

(Cv)−1, wt( v) ≡ 3 (mod 4) .

So now the picture is: if we start with some nondegenerate U with corresponding matrix Φ U ,and then we add in a sequence of vectors v1, . . . , v k of odd weight, each one orthogonal to U and all the previous vectors, then U + ⟨v1, . . . , v k⟩ is also nondegenerate, and Φ U +⟨v1,...,v k ⟩ is simply the product Φ U C±1

v1

· · · C±1

vk

with appropriately chosen signs (depending on the weights of the vi

modulo 4). Each of the C±1

vi

corresponds to a single slice call. Fable points out that my earlier m = 3 example cannot illustrate a nontrivial vector sequence, because the only orthonormal basis for F32 is the standard one. It suggests instead an m = 4 example with orthonormal basis 1110 , 1101 , 1011 , 0111. Then one has C⊗4 = C−11110 C−11101 C−11011 C−10111 .

7.4 Programs, registers and subspaces

In the following sections I will discuss some “programs”, which are a sequence of instructions acting on an array of W “registers”, say ( r1, . . . , r W ). Each register holds an element of ( C2)⊗m. Here are the allowed operations:

ˆ We can apply Cv or C−1

v

to a register, replacing the existing value. Each of these operations costs one slice call (plus data rearrangements, which I remind the reader we are completely ignoring).

ˆ We can add a scalar multiple of one register to another register. For example I could replace

r2 by r2 + 12 r1. As (vaguely) explained above, I am ignoring the cost of these operations. Each register has an associated (nondegenerate) subspace Ui. Initially the subspaces U1, . . . , U W are set to some prescribed values, and then they get updated during the execution of the program. The algorithm does not actually store or manipulate the subspaces; they are just there for accounting purposes. Here is how to think about it: each register ri has an “underlying” value, say fi ∈ (C2)⊗m.What is “physically” stored in ri (i.e., the value actually manipulated by the algorithm) is the partially transformed value Φ Ui fi. To make this work, we have to follow some rules:

ˆ Suppose that we want to enlarge Ui to U ′

i

:= Ui + ⟨v⟩. For this to be allowed, we require that v has odd weight and that v ⊥ Ui. According to the previous section, Φ U ′

i

is Φ Ui Cv if wt( v) = 1 mod 4, and Φ Ui C−1

v

if wt( v) = 3 mod 4. So in the first case we apply Cv to ri,in the second case we apply C−1

v

. In either case, the value in the register is still Φ Ui fi for the new value of Ui. The key thing is that the underlying value fi does not change ; we have changed both the physical representation and the associated Ui in sync.

ˆ Suppose that we want to shrink Ui to U ′

i

:= Ui ∩ ⟨ v⟩⊥. (This is equivalent to enlarging U ⊥

i

to U ⊥

i

• ⟨v⟩.) For this to be allowed, we require that v has odd weight and that v ∈ Ui

(= U ⊥⊥

i

). Similarly to the previous bullet, if wt( v) = 1 mod 4 then we apply C−1

v

to ri, if wt( v) = 3 mod 4 we apply Cv. Again the underlying value fi does not change. 14 ˆ We are only allowed to use the “add scalar multiple” instruction on registers whose Ui are currently equal . In other words, if registers ri and rj are both currently attached to the same subspace U , with physically stored values Φ U fi and Φ U fj , then we are allowed to perform an operation such as replacing rj with rj + 12 ri. In this example we have effectively replaced fj

with fj + 12 fi. Everything stays consistent because Φ U (fj + 12 fi) = Φ U fj + 12 ΦU fi.As an example, suppose that a register begins life with Ui = 0. Since Φ 0 = I, the physically stored value is the same as the underlying value, say fi ∈ (C2)⊗m. If we apply in sequence Ce1 , . . . , C em

to the register, then Ui gets updated to Fm

2

, so that Φ Ui = C⊗m. At this point, the underlying value is still fi, but the value physically in the register is Φ Ui fi = C⊗mfi. So we have succeeded in computing C⊗mfi.In the program described in the next section, some registers will begin life with Ui̸ = 0. For example, they might begin with Ui = ⟨v⟩ for some odd-weight vector v. In this case, the initial underlying value is different to the initial physical value. This does not cause any problems, and in fact the algorithm does not even need to know the initial underlying values. As we will see, the possibility of starting with Ui̸ = 0 is central to making the whole gadget work.

7.5 The baby program

In this section we are going to describe a baby version of the main program. You should think of this as motivation/warmup, rather than as literally a subroutine of the main program (which is presented in the next section). Instead of ( C2)⊗m we work here with ( C2)⊗h for some integer

h ⩾ 8, reserving “ m” for the main program. Let [ h] denote the set {1, . . . , h }.Let T denote the set of subsets of [ h] of cardinality 3. Thus |T| = �h

3

.Let P denote the set of ordered pairs ( S, T ) ∈ T × T such that |S ∩ T | is even, i.e., such that |S ∩ T |

is either 0 or 2. For each S ∈ T there are 3( h − 3) choices of T with |S ∩ T | = 2, and �h−33

choices of T with |S ∩ T | = 0, so altogether |P| = �h

3

� h−33

+ 3( h − 3)
.

Initialisation

There are many registers, playing various different roles. Instead of calling them r1, . . . , r W , we will give them (slightly) more meaningful names. Here they are:

ˆ For each T ∈ T, a “data” register xT .

ˆ For each S ∈ T, a “data” register yS .

ˆ For each ( S, T ) ∈ P, a “side” register AS,T .

ˆ For each i ∈ [h], a “central” register ci.

ˆ One last “central” register c∗.I do not fully understand the meaning of these register names. The total number of registers is

W = 2 |T| + |P| + h + 1 .

15 It is perhaps worth noting that if you count the number of registers in each category, then asymp-totically speaking, the side registers dominate (sextic in h), followed by the data registers (cubic in h), and then only a handful of central registers (linear in h). The initial subspaces are all declared to be zero, except for the xT registers, which are initialised to

⟨vT ⟩. Here vT ∈ Fh

2

means the vector with 1’s at the indices lying in T , and zero elsewhere. Note that wt( vT ) = 3 which is odd. I will use a tilde, for example ˜ xT , to indicate the initial underlying value of a register. At the beginning of the program, here is the state of play: register current U underlying value

xT ⟨vT ⟩ ˜xT

yS 0 ˜yS

AS,T 0 ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

Keep in mind that at this point, the physical values match the underlying values, except for xT ,whose underlying value is the image under Φ −1

⟨vT⟩

of the physical value (i.e., the value actually supplied to the program as input). The program is now given as 8 steps, labelled step 0 to step 7. Each step consists of a transform phase , where we use slice calls to modify the subspaces attached to various registers, and then a

mixup phase , where we apply linear operations to various registers sitting at the same subspace. Then we give a cost field indicating the number of slice calls (we ignore the cost of mixups and other linear operations, and any data rearrangements), and then a table showing the final state of all registers after these operations.

Step 0

Transforms : none.

Mixups : for each ( S, T ) ∈ P, if |S ∩ T | = 0 then subtract 12 AS,T from yS , and if |S ∩ T | = 2 then add 12 AS,T to yS .

Cost : zero.

Final state :register current U underlying value

xT ⟨vT ⟩ ˜xT

yS 0 ˜yS − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

AS,T 0 ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

16 Step 1

Transforms : none.

Mixups : for each S ∈ T, add 12 c∗ − 12

P

i∈S

ci to yS .

Cost : zero.

Final state :register current U underlying value

xT ⟨vT ⟩ ˜xT

yS 0 ˜yS − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

• 12 ˜c∗ − 12

P

i∈S

˜ci

AS,T 0 ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

Step 2

Transforms : for each ( S, T ) ∈ P, apply Φ ⟨vT ⟩ to AS,T . This brings its subspace up to ⟨vT ⟩ (to match xT ). This costs one slice call per pair ( S, T ).

Mixups : for each ( S, T ) ∈ P, add xT to AS,T . (This is legal, as the subspaces now match.)

Cost : |P|.

Final state :register current U underlying value

xT ⟨vT ⟩ ˜xT

yS 0 ˜yS − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

• 12 ˜c∗ − 12

P

i∈S

˜ci

AS,T ⟨vT ⟩ ˜AS,T + ˜ xT

ci 0 ˜ci

c∗ 0 ˜c∗

Step 3

Transforms : for each T ∈ T, we want to apply Φ ⟨vT ⟩⊥ to xT , to bring its subspace up to ⟨vT ⟩ +

⟨vT ⟩⊥ = Fh

2

. To achieve this, we need to break up Φ ⟨vT ⟩⊥ into a sequence of one-dimensional steps, which means we need to find an orthonormal basis for ⟨vT ⟩⊥. Here is an explicit orthonormal basis. Let T = {p, q, r }, and choose any s / ∈ T (possible since h ⩾ 4). Then take the vectors

ep + eq + es, ep + er + es, eq + er + es, ej for j ∈ [h] \ { p, q, r, s }.

The cost is h − 1 slice calls per T . (I want to emphasise here that the reason all of this works, if you trace back through the logic, is that weight is additive mod 4 on orthogonal vectors, and the diagonal entries of the Cv are all powers of i which has order 4. Maybe this is well known to some people, but I have never seen it before and I think it’s very beautiful.) 17 We also want to move ci and c∗ up to Fh

2

. For this we need an orthonormal basis for Fh

2

; of course we can just use the standard basis. The cost is h slice calls for each ci, c ∗.

Mixups : for each T ∈ T, add xT to every ci for which i ∈ T , and also to c∗.

Cost : |T|(h − 1) + ( h + 1) h.

Final state :register current U underlying value

xT Fh

2

˜xT

yS 0 ˜yS − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

• 12 ˜c∗ − 12

P

i∈S

˜ci

AS,T ⟨vT ⟩ ˜AS,T + ˜ xT

ci Fh

2

˜ci + P

T:i∈T

˜xT

c∗ Fh

2

˜c∗ + P

T

˜xT

Step 4

Transforms : we now want to move ci, c ∗ back down to U = 0, to match with yS . This appears to be undoing some of the progress we already made, but the point is that the mixed-in values of ˜ xT

come along for the ride. To achieve this, we need to perform Φ Fh

2

in reverse, i.e., we do Φ −1

⟨ej⟩

for each j ∈ [h]. The cost is h slice calls per ci, c ∗.

Mixups : for each S ∈ T, subtract 12 c∗ − 12

P

i∈S

ci from yS . (This is the inverse of step 1.) This has the following effects on yS :

ˆ It cancels out the ˜ ci and ˜ c∗ terms from yS .

ˆ For each i ∈ S, the ci term contributes 12 ˜xT for each T containing i. The total contribution from all of these i is 12

P

T

|S ∩ T | ˜xT .

ˆ The c∗ term contributes − 12

P

T

˜xT .Thus, after this step, the underlying value of yS is given by ˜yS − 12

X

T:|S∩T|=0

˜AS,T + 12

X

T:|S∩T|=2

˜AS,T + 12

X

T

(|S ∩ T | − 1)˜ xT .

We can simplify the last sum a bit. The only case where |S ∩ T | = 3 is when T = S; this term gives

12

(3 − 1)˜ xS which is simply ˜ xS . Also, the terms with |S ∩ T | = 1 all vanish. So we are left with ˜yS + ˜ xS − 12

X

T:|S∩T|=0

( ˜AS,T + ˜ xT ) + 12

X

T:|S∩T|=2

( ˜AS,T + ˜ xT ).

Cost : ( h + 1) h.

Final state :18 register current U underlying value

xT Fh

2

˜xT

yS 0 ˜yS + ˜ xS − 12

P

T:|S∩T|=0

( ˜AS,T + ˜ xT )+ 12

P

T:|S∩T|=2

( ˜AS,T + ˜ xT )

AS,T ⟨vT ⟩ ˜AS,T + ˜ xT

ci 0 ˜ci + P

T:i∈T

˜xT

c∗ 0 ˜c∗ + P

T∈T

˜xT

Step 5

Transforms : for each S ∈ T we bring yS up to ⟨vS ⟩⊥. We use the same orthonormal basis as in step 3. The cost is h − 1 slice calls per S.Next, for each ( S, T ) ∈ P we want to bring AS,T up to ⟨vS ⟩⊥, to match yS . The crucial fact here is that vT · vS = 0 (since |T ∩ S| = 0 or 2), so vT ∈ ⟨ vS ⟩⊥. (This seems to me to be a really critical point in the argument; it feels like everything else in the gadget might have been designed around this one tiny nugget.) So ⟨vS ⟩⊥ = ⟨vT ⟩ + ⟨vS , v T ⟩⊥ (orthogonal sum), and what we need is an orthonormal basis for ⟨vS , v T ⟩⊥. We can construct one as follows:

ˆ Suppose that |S ∩ T | = 2, say S = {p, q, r } and T = {p, q, s } for r̸ = s. Then we can use

ep + er + es, eq + er + es, ej for j / ∈ { p, q, r, s }.

This is 2 + ( h − 4) = h − 2 vectors altogether.

ˆ Suppose that |S ∩ T | = 0, say S = {p, q, r } and T = {p′, q ′, r ′}, all six indices distinct. Choose distinct s, s ′ /∈ S ∪ T (possible since h ⩾ 8). Then we can use

ep + eq + es, ep + er + es, eq + er + es,ep′ + eq′ + es′ , ep′ + er′ + es′ , eq′ + er′ + es′ ,

together with ej for j / ∈ { p, q, r, s, p ′, q ′, r ′, s ′}. This is 6 + ( h − 8) = h − 2 vectors altogether. The cost is h − 2 slice calls per ( S, T ).

Mixups : for each ( S, T ) ∈ P, if |S ∩ T | = 0 then add 12 AS,T to yS , and if |S ∩ T | = 2 then subtract

12

AS,T from yS . This is the inverse of step 0, and has the effect of removing the ˜AS,T + ˜ xT terms from yS .

Cost : |T|(h − 1) + |P|(h − 2).

Final state :register current U underlying value

xT Fh

2

˜xT

yS ⟨vS ⟩⊥ ˜yS + ˜ xS

AS,T ⟨vS ⟩⊥ ˜AS,T + ˜ xT

ci 0 ˜ci + P

T:i∈T

˜xT

c∗ 0 ˜c∗ + P

T

˜xT

19 Step 6

Transforms : we bring each ci, c ∗ back up to Fh

2

, at a cost of h slice calls for each register.

Mixups : for each T ∈ T, subtract xT from every ci for which i ∈ T , and also from c∗. This is the inverse of step 3.

Cost : ( h + 1) h.

Final state :register current U underlying value

xT Fh

2

˜xT

yS ⟨vS ⟩⊥ ˜yS + ˜ xS

AS,T ⟨vS ⟩⊥ ˜AS,T + ˜ xT

ci Fh

2

˜ci

c∗ Fh

2

˜c∗

Step 7

Transforms : for each ( S, T ) ∈ P, bring AS,T up to Fh

2

by a single application of Φ ⟨vS ⟩. The cost is one slice call per ( S, T ).

Mixups : for each ( S, T ) ∈ P, subtract xT from AS,T . This is the inverse of step 2.

Cost : |P|.

Final state :register current U underlying value

xT Fh

2

˜xT

yS ⟨vS ⟩⊥ ˜yS + ˜ xS

AS,T Fh

2

˜AS,T

ci Fh

2

˜ci

c∗ Fh

2

˜c∗

What did the baby program achieve?

At the end of the baby program, the registers AS,T , ci and c∗ contain their original underlying values, and the attached subspaces have moved from U = 0 to U = Fh

2

. Thus we have succeeded in applying C⊗h to each of them. As for the data registers, we have computed the map ( x, y ) 7 → (x, x + y) on the underlying values. But we have not fully transformed the data, due to the initial subspace on xT having one extra dimension, and the final subspace on yS being one dimension short. Notice however that the ˜ xS part of the final yS register has been fully transformed. Its initial physical value was Φ ⟨vS ⟩(˜ xS ), and its final physical value is Φ ⟨vS ⟩⊥ (˜ xS ), so it has undergone an overall transformation by Φ⟨vS ⟩⊥ Φ−1

⟨vS⟩

= Φ ⟨vS ⟩⊥+⟨vS ⟩Φ−2

⟨vS⟩

= Φ Fh

2

Φ−2

⟨vS⟩

.

20 The Φ Fh

2

factor is just the desired transformation C⊗h. Since wt( vS ) = 3 (mod 4) we have Φ ⟨vS ⟩ =

C−1

vS

, so the Φ −2

⟨vS⟩

factor becomes C2

vS

= XvS , a harmless data rearrangement. This observation is crucial for understanding how the main program in the next section can possibly work out. The total number of slice calls was

|P| + |T|(h − 1) + ( h + 1) h + ( h + 1) h + |T|(h − 1) + |P|(h − 2) + ( h + 1) h + |P|

= |P|h + 2 |T|(h − 1) + 3( h2 + h).

Recalling that W = 2 |T| + |P| + ( h + 1) (the total number of registers), we can rewrite the slice call count as

W h − 2|T| + 2( h2 + h).

Here is how to think about this formula: the W h term is the “expected” number of slice calls ( h

per register), the 2 |T| term is the savings coming from the fact that we didn’t quite fully transform the data registers (we only did h − 1 per data register), and the 2( h2 + h) term is the loss coming from the extra work (“undoing progress”) on the central registers (2 h per register). Since |T| is cubic in h, for large enough h the savings term beats the loss term, and the total becomes strictly less than the expected W h .During Fable’s proofreading of this section, it proactively decided to implement the whole baby program as a sanity check. It confirmed the formulas in the state tables given above for some random inputs, and verified the above formula for the number of slice calls, for a few small values of h. I have not checked the code carefully. Notice that when h = 8 the size of the vectors is 2 8 = 256, still quite manageable for actual computations. For the main program in the next section, such experiments will become infeasible due to combinatorial explosion.

7.6 The main program

If you already had a headache at the end of the last section, I must warn you that it’s about to get worse. Let h be as above. We are going to construct a new program. The registers are now going to be elements of ( C2)⊗m where m := h3. The address space becomes Fm

2

= ( Fh

2

)⊗3. You can think of an address as being an h × h × h array of bits. (The combinatorial explosion is already clear. Even for h = 8 we have m = 512 so the vectors have length 2 512 , well beyond the size of any physically realisable computer in the observable universe.) The symbol T has the same meaning as in the previous section, i.e., size-3 subsets of [ h] :=

{1, . . . , h }. Thus |T| = �h

3

. Similarly P has the same meaning as before. The main program uses the following registers:

ˆ Data registers xa1,a 2,a 3 and ya1,a 2,a 3 for ( a1, a 2, a 3) ∈ T3. So the number of data registers is 2|T|3.It might help to think of x as a big |T| × | T| × | T| cube, where each entry is a register holding an array in ( C2)⊗m. And sitting next to it is a similar cube y.21 ˆ As we will see below, the main program invokes some variant of the baby program 3 |T|2

times. Each invocation needs its own set of scratch registers AS,T , ci and c∗. So altogether we need 3 |T|2(|P| + h + 1) scratch registers. The total number of registers is thus

W := 2 |T|3 + 3 |T|2(|P| + h + 1) .

Before we go through the steps of the program, we will set up some further notation and recall a few useful facts about tensor products. If we have subspaces U ⊆ Fa

2

and V ⊆ Fb

2

, then the subspace U ⊗ V ⊆ Fa

2

⊗ Fb

2

= Fab

2

is the span of the vectors u⊗v for u ∈ U , v ∈ V . A frequently used special case is the fact that ⟨u⟩⊗⟨ v⟩ = ⟨u⊗v⟩.We have (u ⊗ v) · (u′ ⊗ v′) = ( u · u′)( v · v′), u, u ′ ∈ Fa

2

, v, v ′ ∈ Fb

2

.

In particular if {ui} and {vj } are orthonormal bases for U and V , then {ui ⊗ vj } is an orthonormal basis for U ⊗ V .In the rest of this section, whenever we write U ⊕ U ′ we always mean orthogonal direct sum. Also, when we write U ⊥, the ambient space should always be clear from context. If U, V are nondegenerate subspaces of Fh

2

, we have the useful formula (U ⊗ V )⊥ = ( Fh

2

⊗ V ⊥) ⊕ (U ⊥ ⊗ V ). (7.3) This follows from the calculation (Fh

2

⊗ V ⊥) ⊕ (U ⊥ ⊗ V ) ⊕ (U ⊗ V ) = ( Fh

2

⊗ V ⊥) ⊕ (Fh

2

⊗ V ) = Fh

2

⊗ Fh

2

.

The formula (7.3) also holds with the tensor factors switched around; we will sometimes use it this way without further comment. Personally, when I am verifying formulas involving (7.3), instead of doing things algebraically I find it much easier to grab a piece of paper and draw diagrams like this:

We will only need the following types of orthonormal bases: 22 ˆ The standard basis e1, . . . , e h for Fh

2

.

ˆ For T ∈ T, let BT denote the orthonormal basis for ⟨vT ⟩⊥ ⊆ Fh

2

from step 3 of the baby program. It has h − 1 vectors. (This basis is not uniquely determined, it depends on the choice of s / ∈ T . Any choice will do. Similar remarks apply for the bases defined below.)

ˆ For ( S, T ) ∈ P, let BS,T be the orthonormal basis for ⟨vS , v T ⟩⊥ ⊆ Fh

2

from step 5 of the baby program. It has h − 2 vectors. (We continue to assume that h ⩾ 8.)

ˆ For any S, T ∈ T, not necessarily coming from P, we will need an orthonormal basis for

⟨vS ⊗ vT ⟩⊥ ⊆ Fh

2

⊗ Fh

2

. We will denote it by BS⊗T . It may be constructed as follows. From (7.3) we have (⟨vS ⟩ ⊗ ⟨ vT ⟩)⊥ = ( ⟨vS ⟩⊥ ⊗ Fh

2

) ⊕ (⟨vS ⟩ ⊗ ⟨ vT ⟩⊥).

The first component has orthonormal basis {b ⊗ ek : b ∈ BS , k ∈ [h]}. The second component has orthonormal basis {vS ⊗ b : b ∈ BT }. We take BS⊗T to be the union of these two bases. It has h2 − 1 vectors. (Check: |BS || [h]| + |BT | = ( h − 1) h + ( h − 1) = h2 − 1.) The subspaces for the main program are initialised as follows: the register xa1,a 2,a 3 starts at U =

⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩, and all other registers start at U = 0. The symbols ˜ xa1,a 2,a 3 and ˜ ya1,a 2,a 3 refer to the underlying values of these registers at the beginning of the program. The main program has three “stages”, which we now describe in turn.

Stage 1 — overview

For each pair ( a2, a 3) ∈ T × T we perform 9 steps. The first 8 steps are basically the same as the baby program, acting on the data registers xT,a 2,a 3

and yS,a 2,a 3 , where T, S ∈ T. So you can think of stage 1 as operating on one-dimensional “slices” of the |T| × | T| × | T| cubes x and y, along the first dimension. We use a new set of scratch registers (AS,T , ci, c∗) for each pair ( a2, a 3). Strictly speaking we should add subscripts a2, a 3 to the scratch registers, but in the interests of readability we will omit them. Actually, it doesn’t quite make sense to apply the baby program in its original form, because it was stated for arrays in ( C2)⊗h (with addresses in Fh

2

), whereas the registers here take values in (C2)⊗m (with addresses in Fm

2

= ( Fh

2

)⊗3). Instead, we have to modify the slice calls: whenever we would have applied Φ ⟨v⟩ to some register, for some v ∈ Fh

2

of odd weight, we replace this by Φ ⟨v⊗va2 ⊗va3 ⟩ (and similarly for Φ −1

⟨v⟩

). In other words, every subspace U ⊆ Fh

2

in the baby program gets replaced by U ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ⊆ Fm

2

, and the latter has the same dimension as U .The replacement v 7 → v ⊗ va2 ⊗ va3 preserves all orthonormality relations, so all arguments about legality of slice call sequences still hold. (Incidentally, in the baby program, recall that we only ever used v with weight 1 or 3. Since wt( va2 ) = wt( va3 ) = 3, we see that wt( v ⊗ va2 ⊗ va3 ) = 9 wt( v). Thus the modified version only ever performs slice calls in directions v ⊗ va2 ⊗ va3 with weight 9 or 27.) Finally, in step 8 we apply some extra Φ ⟨v⟩’s to the scratch registers to bump their subspaces up from Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ to the full Fm

2

.23 We now briefly summarise these steps. All 9 steps are performed once for each pair ( a2, a 3). The first 8 steps are basically copied from the baby program with all spaces/vectors tensored by

va2 ⊗ va3 appropriately, so we won’t bother writing out the transform and mixup phases. The costs are literally the same as in the baby program.

Stage 1, initialisation

Here is the initial state: register current U underlying value

xT,a 2,a 3 ⟨vT ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 0 ˜yS,a 2,a 3

AS,T 0 ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

Stage 1, step 0

Cost : zero.

Final state :register current U underlying value

xT,a 2,a 3 ⟨vT ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 0 ˜yS,a 2,a 3 − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

AS,T 0 ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

Stage 1, step 1

Cost : zero.

Final state :register current U underlying value

xT,a 2,a 3 ⟨vT ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 0 ˜yS,a 2,a 3 − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

• 12 ˜c∗ − 12

P

i∈S

˜ci

AS,T 0 ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

Stage 1, step 2

Cost : |P|.

Final state :24 register current U underlying value

xT,a 2,a 3 ⟨vT ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 0 ˜yS,a 2,a 3 − 12

P

T :|S∩T |=0 ˜AS,T + 12

P

T :|S∩T |=2 ˜AS,T

• 12 ˜c∗ − 12

P

i∈S ˜ci

AS,T ⟨vT ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ xT,a 2,a 3

ci 0 ˜ci

c∗ 0 ˜c∗

Stage 1, step 3

Cost : |T|(h − 1) + ( h + 1) h.

Final state :register current U underlying value

xT,a 2,a 3 Fh

2 ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 0 ˜yS,a 2,a 3 − 12

P

T :|S∩T |=0 ˜AS,T + 12

P

T :|S∩T |=2 ˜AS,T

• 12 ˜c∗ − 12

P

i∈S ˜ci

AS,T ⟨vT ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ xT,a 2,a 3

ci Fh

2 ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜ci + P

T :i∈T ˜xT,a 2,a 3

c∗ Fh

2 ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜c∗ + P

T ˜xT,a 2,a 3

Stage 1, step 4

Cost : ( h + 1) h.

Final state :register current U underlying value

xT,a 2,a 3 Fh

2 ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 0 ˜yS,a 2,a 3 + ˜ xS,a 2,a 3 − 12

P

T :|S∩T |=0 ( ˜AS,T + ˜ xT,a 2,a 3 )+ 12

P

T :|S∩T |=2 ( ˜AS,T + ˜ xT,a 2,a 3 )

AS,T ⟨vT ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ xT,a 2,a 3

ci 0 ˜ci + P

T :i∈T ˜xT,a 2,a 3

c∗ 0 ˜c∗ + P

T ˜xT,a 2,a 3

Stage 1, step 5

Cost : |T|(h − 1) + |P|(h − 2).

Final state :25 register current U underlying value

xT,a 2,a 3 Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 ⟨vS ⟩⊥ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜yS,a 2,a 3 + ˜ xS,a 2,a 3

AS,T ⟨vS ⟩⊥ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ xT,a 2,a 3

ci 0 ˜ci + P

T:i∈T

˜xT,a 2,a 3

c∗ 0 ˜c∗ + P

T

˜xT,a 2,a 3

Stage 1, step 6

Cost : ( h + 1) h.

Final state :register current U underlying value

xT,a 2,a 3 Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 ⟨vS ⟩⊥ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜yS,a 2,a 3 + ˜ xS,a 2,a 3

AS,T ⟨vS ⟩⊥ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ xT,a 2,a 3

ci Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜ci

c∗ Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜c∗

Stage 1, step 7

Cost : |P|.

Final state :register current U underlying value

xT,a 2,a 3 Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 ⟨vS ⟩⊥ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜yS,a 2,a 3 + ˜ xS,a 2,a 3

AS,T Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜AS,T

ci Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜ci

c∗ Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜c∗

Stage 1, step 8

(This step does not come from the baby program.)

Transforms : we want to move each scratch register to Fh

2

⊗ Fh

2

⊗ Fh

2

= Fm

2

, so we need to add the space Fh

2

⊗ (⟨va2 ⟩ ⊗ ⟨ va3 ⟩)⊥. We may use the orthonormal basis {ej ⊗ b ⊗ c : j ∈ [h], b ⊗ c ∈ Ba2⊗a3 }.This costs h(h2 − 1) slice calls per scratch register.

Mixups : none.

Cost : ( |P| + h + 1) · h(h2 − 1).

Final state :26 register current U underlying value

xT,a 2,a 3 Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜xT,a 2,a 3

yS,a 2,a 3 ⟨vS ⟩⊥ ⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ ˜yS,a 2,a 3 + ˜ xS,a 2,a 3

AS,T Fh

2

⊗ Fh

2

⊗ Fh

2

˜AS,T

ci Fh

2

⊗ Fh

2

⊗ Fh

2

˜ci

c∗ Fh

2

⊗ Fh

2

⊗ Fh

2

˜c∗

At this point we are done with this set of scratch registers; we put them aside and bring in a new set for the next pair ( a2, a 3) (starting again of course at U = 0 for each register).

Stage 1 — conclusion

As for the data registers, we have computed the map ( x, y ) 7 → (x, x + y) on the underlying values. The first component of the subspace for y is one dimension short of the full Fh

2

⊗ ⟨ va2 ⟩ ⊗ ⟨ va3 ⟩ (the current “target”), so the ˜ y part has not been fully transformed. But, as explained in the baby program, the ˜ x part has been fully transformed (i.e., up to the current target). The scratch registers were transformed completely, i.e., we have applied C⊗m to each of them. The total number of slice calls for Stage 1 is just the baby program count, plus the cost of step 8, all multiplied by the number of pairs ( a2, a 3):

|T|2 �|P|h + 2 |T|(h − 1) + 3( h2 + h) + ( |P| + h + 1) · h(h2 − 1)

= 2 |T|3(h − 1) + |T|2|P|h3 + |T|2(h2 + 2)( h2 + h).

Stage 2 — overview

The aim of this stage is to compute the function ( x, y ) 7 → (x − y, y ) on each pair of data registers. This maps the current values (˜ xa, ˜ya + ˜ xa) (see last state table) to (˜ xa − (˜ ya + ˜ xa), ˜ya + ˜ xa) = ( −˜ya, ˜ya + ˜ xa).

The strategy is essentially to run the steps of the baby program in reverse order, with each mixup operation inverted, and with the roles of x and y switched. However, we cannot run the Φ ⟨v⟩’s backwards, because we are still trying to make progress in the forward direction on the subspaces. As in the baby program, the main requirement is to ensure that the subspaces agree whenever we want to perform a linear operation on two registers. For each pair ( a1, a 3) ∈ T × T, we perform the following steps. (The modified baby program is going to act on the second coordinate.)

Stage 2, initialisation

The data registers in the following table are copied from the final state of Stage 1. For each ( a1, a 3)we use a new set of scratch registers, initialised to U = 0. 27 register current U underlying value

xa1,S,a 3 Fh

2

⊗ ⟨ vS ⟩ ⊗ ⟨ va3 ⟩ ˜xa1,S,a 3

ya1,T,a 3 ⟨va1 ⟩⊥ ⊗ ⟨ vT ⟩ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T 0 ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

Stage 2, step 0 (inverse of baby program step 7)

Transforms : we bring AS,T up to ⟨va1 ⟩⊥ ⊗ ⟨ vT ⟩ ⊗ ⟨ va3 ⟩ (to match ya1,T,a 3 ) via the orthonormal basis

{b ⊗ vT ⊗ va3 : b ∈ Ba1 }. This needs h − 1 slice calls for each ( S, T ).

Mixups : for each ( S, T ) ∈ P, we add ya1,T,a 3 to AS,T .

Cost : |P|(h − 1).

Final state :register current U underlying value

xa1,S,a 3 Fh

2

⊗ ⟨ vS ⟩ ⊗ ⟨ va3 ⟩ ˜xa1,S,a 3

ya1,T,a 3 ⟨va1 ⟩⊥ ⊗ ⟨ vT ⟩ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T ⟨va1 ⟩⊥ ⊗ ⟨ vT ⟩ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ ya1,T,a 3 + ˜ xa1,T,a 3

ci 0 ˜ci

c∗ 0 ˜c∗

Stage 2, step 1 (inverse of baby program step 6)

Transforms : we want to bring ya1,T,a 3 and ci, c ∗ up to ⟨va1 ⟩⊥ ⊗ Fh

2

⊗ ⟨ va3 ⟩.For ya1,T,a 3 we use the orthonormal basis {b ⊗ c ⊗ va3 : b ∈ Ba1 , c ∈ BT }. This costs ( h − 1) 2 slice calls for each T ∈ T.For ci, c ∗ we use the orthonormal basis {b ⊗ ej ⊗ va3 : b ∈ Ba1 , j ∈ [h]}. This costs ( h − 1) h slice calls for each ci, c ∗.

Mixups : for each T ∈ T, we add ya1,T,a 3 to every ci for which i ∈ T , and to c∗.

Cost : |T|(h − 1) 2 + ( h − 1) h(h + 1).

Final state :register current U underlying value

xa1,S,a 3 Fh

2

⊗ ⟨ vS ⟩ ⊗ ⟨ va3 ⟩ ˜xa1,S,a 3

ya1,T,a 3 ⟨va1 ⟩⊥ ⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T ⟨va1 ⟩⊥ ⊗ ⟨ vT ⟩ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ ya1,T,a 3 + ˜ xa1,T,a 3

ci ⟨va1 ⟩⊥ ⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜ci + P

T:i∈T

(˜ ya1,T,a 3 + ˜ xa1,T,a 3 )

c∗ ⟨va1 ⟩⊥ ⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜c∗ + P

T

(˜ ya1,T,a 3 + ˜ xa1,T,a 3 )28 Stage 2, step 2 (inverse of baby program step 5)

Transforms : we want to bring both xa1,S,a 3 and AS,T up to the subspace

�⟨va1 ⟩ ⊗ ⟨ vS ⟩⊥
⊥ ⊗ ⟨ va3 ⟩.

By (7.3) this is the same as

�Fh

2

⊗ ⟨ vS ⟩ ⊗ ⟨ va3 ⟩
⊕ �⟨va1 ⟩⊥ ⊗ ⟨ vS ⟩⊥ ⊗ ⟨ va3 ⟩
.

For xa1,S,a 3 we just need to add on the ⟨va1 ⟩⊥ ⊗ ⟨ vS ⟩⊥ ⊗ ⟨ va3 ⟩ part; we can use the orthonormal basis {b ⊗ c ⊗ va3 : b ∈ Ba1 , c ∈ BS }, which requires ( h − 1) 2 slice calls. For AS,T , we first add ⟨va1 ⟩⊥ ⊗ ⟨ vS , v T ⟩⊥ ⊗ ⟨ va3 ⟩ using the orthonormal basis {b ⊗ c ⊗ va3 : b ∈

Ba1 , c ∈ BS,T }, which gets us to ⟨va1 ⟩⊥ ⊗ ⟨ vS ⟩⊥ ⊗ ⟨ va3 ⟩. (The key fact here is that vS · vT = 0.) We then add the Fh

2

⊗ ⟨ vS ⟩ ⊗ ⟨ va3 ⟩ part using the orthonormal basis {ej ⊗ vS ⊗ va3 : j ∈ [h]}. The total number of slice calls is ( h − 1)( h − 2) + h = h2 − 2h + 2.

Mixups : for each ( S, T ) ∈ P, if |S ∩ T | = 0 then subtract 12 AS,T from xa1,S,a 3 , and if |S ∩ T | = 2 then add 12 AS,T to xa1,S,a 3 .

Cost : |T|(h − 1) 2 + |P|(h2 − 2h + 2).

Final state (taking the opportunity to rewrite all of the subspaces in terms of complements for the first two components): register current U underlying value

xa1,S,a 3 (⟨va1 ⟩ ⊗ ⟨ vS ⟩⊥)⊥ ⊗ ⟨ va3 ⟩ ˜xa1,S,a 3 − 12

P

T:|S∩T|=0

( ˜AS,T + ˜ ya1,T,a 3 + ˜ xa1,T,a 3 )+ 12

P

T:|S∩T|=2

( ˜AS,T + ˜ ya1,T,a 3 + ˜ xa1,T,a 3 )

ya1,T,a 3 (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T (⟨va1 ⟩ ⊗ ⟨ vS ⟩⊥)⊥ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ ya1,T,a 3 + ˜ xa1,T,a 3

ci (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜ci + P

T:i∈T

(˜ ya1,T,a 3 + ˜ xa1,T,a 3 )

c∗ (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜c∗ + P

T

(˜ ya1,T,a 3 + ˜ xa1,T,a 3 )

Stage 2, step 3 (inverse of baby program step 4)

Transforms : we want to bring xa1,S,a 3 and ci, c∗ up to Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩.For xa1,S,a 3 this requires adding on ⟨va1 ⟩ ⊗ ⟨ vS ⟩⊥ ⊗ ⟨ va3 ⟩, for which we use the orthonormal basis

{va1 ⊗ b ⊗ va3 : b ∈ BS }. This requires h − 1 slice calls. For ci, c∗ this requires adding on ⟨va1 ⟩ ⊗ Fh

2

⊗ ⟨ va3 ⟩. We use the orthonormal basis {va1 ⊗ ej ⊗ va3 :

j ∈ [h]}. This requires h slice calls.

Mixups : for each S ∈ T we add 12 c∗ − 12

P

i∈S

ci to xa1,S,a 3 . A similar argument to step 4 of the baby program shows that most ˜ x and ˜ y terms drop out, and the ˜ xa1,S,a 3 gets replaced by −˜ya1,S,a 3 .

Cost : |T|(h − 1) + ( h + 1) h.

Final state :29 register current U underlying value

xa1,S,a 3 Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ −˜ya1,S,a 3 − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

• 12 ˜c∗ − 12

P

i∈S

˜ci

ya1,T,a 3 (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T (⟨va1 ⟩ ⊗ ⟨ vS ⟩⊥)⊥ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ ya1,T,a 3 + ˜ xa1,T,a 3

ci Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜ci + P

T:i∈T

(˜ ya1,T,a 3 + ˜ xa1,T,a 3 )

c∗ Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜c∗ + P

T∈T

(˜ ya1,T,a 3 + ˜ xa1,T,a 3 )

Stage 2, step 4 (inverse of baby program step 3)

Transforms : we want to move ci, c∗ down to ( ⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩, to line up with ya1,T,a 3 . This is the only backwards move in Stage 2. We need to subtract out the subspace ⟨va1 ⟩ ⊗ Fh

2

⊗ ⟨ va3 ⟩, so we use the orthonormal basis {va1 ⊗ ej ⊗ va3 : j ∈ [h]}. The number of slice calls is h.

Mixups : for each T ∈ T we subtract ya1,T,a 3 from every ci with i ∈ T , and also from c∗.

Cost : ( h + 1) h.

Final state :register current U underlying value

xa1,S,a 3 Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ −˜ya1,S,a 3 − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

• 12 ˜c∗ − 12

P

i∈S

˜ci

ya1,T,a 3 (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T (⟨va1 ⟩ ⊗ ⟨ vS ⟩⊥)⊥ ⊗ ⟨ va3 ⟩ ˜AS,T + ˜ ya1,T,a 3 + ˜ xa1,T,a 3

ci (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜ci

c∗ (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜c∗

Stage 2, step 5 (inverse of baby program step 2)

Transforms : we want to bring both ya1,T,a 3 and AS,T up to ( ⟨va1 ⟩ ⊗ ⟨ vT ⟩)⊥ ⊗ ⟨ va3 ⟩.For ya1,T,a 3 this requires adding on the subspace ⟨va1 ⟩ ⊗ ⟨ vT ⟩⊥ ⊗ ⟨ va3 ⟩. We use the orthonormal basis {va1 ⊗ b ⊗ va3 : b ∈ BT }. The number of slice calls is h − 1. For AS,T , first observe that by (7.3), (⟨va1 ⟩ ⊗ ⟨ vS ⟩⊥)⊥ = ( ⟨va1 ⟩ ⊗ ⟨ vS ⟩) ⊕ (⟨va1 ⟩⊥ ⊗ Fh

2

).

This implies (by vS · vT = 0 and (7.3) again) that (⟨va1 ⟩ ⊗ ⟨ vS ⟩⊥)⊥ ⊕ (⟨va1 ⟩ ⊗ ⟨ vS , v T ⟩⊥) = ( ⟨va1 ⟩ ⊗ ⟨ vT ⟩⊥) ⊕ (⟨va1 ⟩⊥ ⊗ Fh

2

)= ( ⟨va1 ⟩ ⊗ ⟨ vT ⟩)⊥.

Thus, to get AS,T up to ( ⟨va1 ⟩⊗⟨ vT ⟩)⊥⊗⟨ va3 ⟩, it suffices to add the subspace ⟨va1 ⟩⊗⟨ vS , v T ⟩⊥⊗⟨ va3 ⟩.We use the orthonormal basis {va1 ⊗ b ⊗ va3 : b ∈ BS,T }. The number of slice calls is h − 2.

Mixups : for each ( S, T ) ∈ P we subtract ya1,T,a 3 from AS,T .30 Cost : |T|(h − 1) + |P|(h − 2).

Final state :register current U underlying value

xa1,S,a 3 Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ −˜ya1,S,a 3 − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

• 12 ˜c∗ − 12

P

i∈S

˜ci

ya1,T,a 3 (⟨va1 ⟩ ⊗ ⟨ vT ⟩)⊥ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T (⟨va1 ⟩ ⊗ ⟨ vT ⟩)⊥ ⊗ ⟨ va3 ⟩ ˜AS,T

ci (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜ci

c∗ (⟨va1 ⟩ ⊗ Fh

2

)⊥ ⊗ ⟨ va3 ⟩ ˜c∗

Stage 2, step 6 (inverse of baby program step 1)

Transforms : we bring ci, c ∗ back up to Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩. As in step 4, the number of slice calls is h.

Mixups : for each S ∈ T we subtract 12 c∗ − 12

P

i∈S

ci from xa1,S,a 3 .

Cost : ( h + 1) h.

Final state :register current U underlying value

xa1,S,a 3 Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ −˜ya1,S,a 3 − 12

P

T:|S∩T|=0

˜AS,T + 12

P

T:|S∩T|=2

˜AS,T

ya1,T,a 3 (⟨va1 ⟩ ⊗ ⟨ vT ⟩)⊥ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T (⟨va1 ⟩ ⊗ ⟨ vT ⟩)⊥ ⊗ ⟨ va3 ⟩ ˜AS,T

ci Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜ci

c∗ Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜c∗

Stage 2, step 7 (inverse of baby program step 0)

Transforms : we bring AS,T up to Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩. Only one slice call is needed here, to add the subspace ⟨va1 ⟩ ⊗ ⟨ vT ⟩ ⊗ ⟨ va3 ⟩.

Mixups : for each ( S, T ) ∈ P, if |S ∩ T | = 0 then add 12 AS,T to xa1,S,a 3 , and if |S ∩ T | = 2 then subtract 12 AS,T from xa1,S,a 3 .

Cost : |P|.

Final state :register current U underlying value

xa1,S,a 3 Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ −˜ya1,S,a 3

ya1,T,a 3 (⟨va1 ⟩ ⊗ ⟨ vT ⟩)⊥ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜AS,T

ci Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜ci

c∗ Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ ˜c∗

31 Stage 2, step 8

Transforms : as in step 8 of Stage 1, we need to finish transforming the scratch registers, to get to

Fh

2

⊗ Fh

2

⊗ Fh

2

. This requires adding the subspace Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩⊥. We use the orthonormal basis

{ej ⊗ ek ⊗ b : j, k ∈ [h], b ∈ Ba3 }. This costs h2(h − 1) slice calls.

Mixups : none.

Cost : ( |P| + h + 1) h2(h − 1).

Final state :register current U underlying value

xa1,S,a 3 Fh

2

⊗ Fh

2

⊗ ⟨ va3 ⟩ −˜ya1,S,a 3

ya1,T,a 3 (⟨va1 ⟩ ⊗ ⟨ vT ⟩)⊥ ⊗ ⟨ va3 ⟩ ˜ya1,T,a 3 + ˜ xa1,T,a 3

AS,T Fh

2

⊗ Fh

2

⊗ Fh

2

˜AS,T

ci Fh

2

⊗ Fh

2

⊗ Fh

2

˜ci

c∗ Fh

2

⊗ Fh

2

⊗ Fh

2

˜c∗

As in stage 1, we are now done with these scratch registers, and we bring in a fresh set for the next pair ( a1, a 3).

Stage 2 — conclusion

Again, the subspace for ya1,T,a 3 is still one dimension short of the stage 2 target, which is Fh

2

⊗ Fh

2

⊗⟨va3 ⟩. So the ˜ y component of y is not quite fully transformed. But the ˜ x part is fully transformed (just as in stage 1), and the important new feature is that the ˜ y value in the x register is also fully transformed up to the current target. The total slice call count for each pair ( a1, a 3) is

|T| �(h − 1) 2 + ( h − 1) 2 + ( h − 1) + ( h − 1)

• |P| �(h − 1) + ( h2 − 2h + 2) + ( h − 2) + 1 + h2(h − 1)
  

• �(h − 1) h(h + 1) + ( h + 1) h + ( h + 1) h + ( h + 1) h + ( h + 1) h2(h − 1)
  

= 2 |T|(h2 − h) + |P|h3 + ( h2 + 2)( h2 + h).

Thus the total cost of Stage 2 is 2|T|3(h2 − h) + |T|2|P|h3 + |T|2(h2 + 2)( h2 + h).

The combined cost of Stage 1 and Stage 2 is 2|T|3(h2 − 1) + 2 |T|2|P|h3 + 2 |T|2(h2 + 2)( h2 + h).

Stage 3 — overview

In this stage, we want to compute the original baby program map ( x, y ) 7 → (x, y + x) again, to clear out the ˜ y part of the y registers. We are going to run the baby program in forward order, now 32 operating on the third components. The Φ ⟨v⟩’s need to be modified, but at least this time they are all going in the right direction. For each ( a1, a 2) ∈ T × T we perform the following steps. (The numbering matches that of the baby program.)

Stage 3, initialisation

The data registers in the following table are copied from the final state of Stage 2. For each ( a1, a 2)we use a new set of scratch registers, initialised to U = 0. register current U underlying value

xa1,a 2,T Fh

2

⊗ Fh

2

⊗ ⟨ vT ⟩ −˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ ⟨ vS ⟩ ˜ya1,a 2,S + ˜ xa1,a 2,S

AS,T 0 ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

Stage 3, step 0

Transforms : we bring AS,T up to ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ ⟨ vS ⟩ (to match ya1,a 2,S ), using the orthonormal basis {b ⊗ c ⊗ vS : b ⊗ c ∈ Ba1⊗a2 }. The number of slice calls is h2 − 1.

Mixups : for each ( S, T ) ∈ P, if |S ∩ T | = 0 then subtract 12 AS,T from ya1,a 2,S , and if |S ∩ T | = 2 then add 12 AS,T to ya1,a 2,S .

Cost : |P|(h2 − 1).

Final state :register current U underlying value

xa1,a 2,T Fh

2

⊗ Fh

2

⊗ ⟨ vT ⟩ −˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ ⟨ vS ⟩ ˜ya1,a 2,S + ˜ xa1,a 2,S − 12

P

T:|S∩T|=0

˜AS,T

• 12

P

T:|S∩T|=2

˜AS,T

AS,T (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ ⟨ vS ⟩ ˜AS,T

ci 0 ˜ci

c∗ 0 ˜c∗

Stage 3, step 1

Transforms : we want to move ya1,a 2,S and ci, c∗ up to ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

.For ya1,a 2,S this involves adding the subspace ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ ⟨ vS ⟩⊥. We use the orthonormal basis {b ⊗ c ⊗ w : b ⊗ c ∈ Ba1⊗a2 , w ∈ BS }. The number of slice calls is ( h2 − 1)( h − 1). For ci, c ∗ we have to add on all of ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

. We use the orthonormal basis {b ⊗ c ⊗ ej :

b ⊗ c ∈ Ba1⊗a2 , j ∈ [h]}. The number of slice calls is ( h2 − 1) h.

Mixups : for each S ∈ T we add 12 c∗ − 12

P

i∈S

ci to ya1,a 2,S .33 Cost : |T|(h2 − 1)( h − 1) + ( h2 − 1) h(h + 1).

Final state :register current U underlying value

xa1,a 2,T Fh

2

⊗ Fh

2

⊗ ⟨ vT ⟩ −˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜ya1,a 2,S + ˜ xa1,a 2,S − 12

P

T:|S∩T|=0

˜AS,T

• 12

P

T:|S∩T|=2

˜AS,T + 12 ˜c∗ − 12

P

i∈S

˜ci

AS,T (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ ⟨ vS ⟩ ˜AS,T

ci (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜ci

c∗ (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜c∗

Stage 3, step 2

Transforms : we want to bring xa1,a 2,T and AS,T up to ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩⊥)⊥.For xa1,a 2,T this requires adding the subspace ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ ⟨ vT ⟩⊥ (by (7.3)). We use the orthonormal basis {b ⊗ c ⊗ w : b ⊗ c ∈ Ba1⊗a2 , w ∈ BT }. This needs ( h2 − 1)( h − 1) slice calls. For AS,T we must first add the subspace ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ ⟨ vS ⟩⊥ to get up to ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

.We use the orthonormal basis {b ⊗ c ⊗ w : b ⊗ c ∈ Ba1⊗a2 , w ∈ BS }. This requires ( h2 − 1)( h − 1) slice calls. Then we add the one-dimensional space ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩, which costs only one slice call.

Mixups : for each ( S, T ) ∈ P, add xa1,a 2,T to AS,T .

Cost : |T|(h2 − 1)( h − 1) + |P|(( h2 − 1)( h − 1) + 1).

Final state :register current U underlying value

xa1,a 2,T (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩⊥)⊥ −˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜ya1,a 2,S + ˜ xa1,a 2,S − 12

P

T:|S∩T|=0

˜AS,T

• 12

P

T:|S∩T|=2

˜AS,T + 12 ˜c∗ − 12

P

i∈S

˜ci

AS,T (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩⊥)⊥ ˜AS,T − ˜ya1,a 2,T

ci (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜ci

c∗ (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜c∗

Stage 3, step 3

Transforms : we want to move xa1,a 2,T and ci, c ∗ up to Fh

2

⊗ Fh

2

⊗ Fh

2

(the fully transformed state). For xa1,a 2,T we need to add the subspace ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩⊥. We use the orthonormal basis

{va1 ⊗ va2 ⊗ b : b ∈ BT }, at a cost of h − 1 slice calls. For ci, c ∗ we need to add the subspace ⟨va1 ⟩⊗⟨ va2 ⟩⊗ Fh

2

. We use the orthonormal basis {va1 ⊗va2 ⊗ej :

j ∈ [h]}. This costs h slice calls.

Mixups : for each T ∈ T, add xa1,a 2,T to every ci for which i ∈ T , and also to c∗.34 Cost : |T|(h − 1) + ( h + 1) h.

Final state :register current U underlying value

xa1,a 2,T Fh

2

⊗ Fh

2

⊗ Fh

2

−˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜ya1,a 2,S + ˜ xa1,a 2,S − 12

P

T:|S∩T|=0

˜AS,T

• 12

P

T:|S∩T|=2

˜AS,T + 12 ˜c∗ − 12

P

i∈S

˜ci

AS,T (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩⊥)⊥ ˜AS,T − ˜ya1,a 2,T

ci Fh

2

⊗ Fh

2

⊗ Fh

2

˜ci − P

T:i∈T

˜ya1,a 2,T

c∗ Fh

2

⊗ Fh

2

⊗ Fh

2

˜c∗ − P

T

˜ya1,a 2,T

Stage 3, step 4

Transforms : we bring ci, c ∗ back down to ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

. This is the only backwards move in Stage 3. We use the same orthonormal basis as step 3. This costs h slice calls.

Mixups : for each S ∈ T, subtract 12 c∗ − 12

P

i∈S

ci from ya1,a 2,S .

Cost : ( h + 1) h.

Final state :register current U underlying value

xa1,a 2,T Fh

2

⊗ Fh

2

⊗ Fh

2

−˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜xa1,a 2,S − 12

P

T:|S∩T|=0

( ˜AS,T − ˜ya1,a 2,T )+ 12

P

T:|S∩T|=2

( ˜AS,T − ˜ya1,a 2,T )

AS,T (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩⊥)⊥ ˜AS,T − ˜ya1,a 2,T

ci (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜ci − P

T:i∈T

˜ya1,a 2,T

c∗ (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜c∗ − P

T

˜ya1,a 2,T

Stage 3, step 5

Transforms : we want to bring ya1,a 2,S and AS,T up to ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩)⊥.For ya1,a 2,S we need to add the subspace ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩⊥. We use the orthonormal basis

{va1 ⊗ va2 ⊗ b : b ∈ BS }, at a cost of h − 1 slice calls. For AS,T we need only add the subspace ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS , v T ⟩⊥, as (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩⊥)⊥ = ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩) ⊕ (( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

)together with vS · vT = 0 implies that (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vT ⟩⊥)⊥ ⊕ (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS , v T ⟩⊥)= ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩⊥) ⊕ (( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

) = ( ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩)⊥.

Thus we may use the orthonormal basis {va1 ⊗ va2 ⊗ b : b ∈ BS,T }, at a cost of h − 2 slice calls. 35 Mixups : for each ( S, T ) ∈ P, if |S ∩ T | = 0 then add 12 AS,T to ya1,a 2,S , and if |S ∩ T | = 2 then subtract 12 AS,T from ya1,a 2,S .

Cost : |T|(h − 1) + |P|(h − 2).

Final state :register current U underlying value

xa1,a 2,T Fh

2

⊗ Fh

2

⊗ Fh

2

−˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩)⊥ ˜xa1,a 2,S

AS,T (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩)⊥ ˜AS,T − ˜ya1,a 2,T

ci (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜ci − P

T:i∈T

˜ya1,a 2,T

c∗ (⟨va1 ⟩ ⊗ ⟨ va2 ⟩)⊥ ⊗ Fh

2

˜c∗ − P

T

˜ya1,a 2,T

Stage 3, step 6

Transforms : bring ci, c ∗ up to Fh

2

⊗ Fh

2

⊗ Fh

2

by adding the subspace ⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ Fh

2

. We use the orthonormal basis {va1 ⊗ va2 ⊗ ej : j ∈ [h]}, at a cost of h slice calls.

Mixups : for each T ∈ T, subtract xa1,a 2,T from every ci for which i ∈ T , and also from c∗.

Cost : ( h + 1) h.

Final state :register current U underlying value

xa1,a 2,T Fh

2

⊗ Fh

2

⊗ Fh

2

−˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩)⊥ ˜xa1,a 2,S

AS,T (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩)⊥ ˜AS,T − ˜ya1,a 2,T

ci Fh

2

⊗ Fh

2

⊗ Fh

2

˜ci

c∗ Fh

2

⊗ Fh

2

⊗ Fh

2

˜c∗

Stage 3, step 7

Transforms : bring AS,T up to Fh

2

⊗ Fh

2

⊗ Fh

2

by a single step in the direction va1 ⊗ va2 ⊗ vS . This costs one slice call.

Mixups : for each ( S, T ) ∈ P, subtract xa1,a 2,T from AS,T .

Cost : |P|.

Final state :register current U underlying value

xa1,a 2,T Fh

2

⊗ Fh

2

⊗ Fh

2

−˜ya1,a 2,T

ya1,a 2,S (⟨va1 ⟩ ⊗ ⟨ va2 ⟩ ⊗ ⟨ vS ⟩)⊥ ˜xa1,a 2,S

AS,T Fh

2

⊗ Fh

2

⊗ Fh

2

˜AS,T

ci Fh

2

⊗ Fh

2

⊗ Fh

2

˜ci

c∗ Fh

2

⊗ Fh

2

⊗ Fh

2

˜c∗

36 At this point the scratch registers are fully transformed, and we bring in new ones for the next pair (a1, a 2).

Stage 3 — conclusion

Let us analyse what has happened to the physical values over the course of all three stages.

ˆ All 3 |T|2(|P| + h + 1) scratch registers have been fully transformed.

ˆ The initial values ˜ ya1,a 2,a 3 , which started out in the y registers, have ended up in the x

registers, fully transformed, and with a harmless minus sign.

ˆ The initial values ˜ xa1,a 2,a 3 , which started out in the x registers, have ended up in the y

registers. They have been fully transformed, up to a rearrangement, by the same argument used in the baby case (see the section “What did the baby program achieve?”). Therefore, the net result is that every single register has been fully transformed (modulo free rearrangements and sign changes). The number of slice calls in Stage 3 for each pair ( a1, a 2) is

|T|�(h2 − 1)( h − 1) + ( h2 − 1)( h − 1) + ( h − 1) + ( h − 1)

• |P|�(h2 − 1) + ( h2 − 1)( h − 1) + 1 + ( h − 2) + 1
  

• �(h2 − 1) h(h + 1) + ( h + 1) h + ( h + 1) h + ( h + 1) h
  

= 2 |T|(h3 − h2) + |P|h3 + ( h2 + 2)( h2 + h).

Thus the total number of slice calls in Stage 3 over all pairs is 2|T|3(h3 − h2) + |T|2|P|h3 + |T|2(h2 + 2)( h2 + h).

Summing over all three stages, the grand total is 2|T|3(h3 − 1) + 3 |T|2|P|h3 + 3 |T|2(h2 + 2)( h2 + h).

Recalling that the number of registers is W = 2 |T|3 + 3 |T|2(|P| + h + 1), after a bit of algebra we can rewrite the grand total as

W h 3 − | T|2(2 |T| − 6h(h + 1)) .

Since |T| is cubic in h, for sufficiently large h the bracketed term is positive (according to Fable, for h ⩾ 22), and we see that the total number of slice calls is strictly less than W h 3 = W m . And this is exactly what the tensor amplification wants. Unbelievable.

8 Conclusion

As far as I can see, the gadget does what it’s supposed to do. I have worked through all the steps pretty carefully, and I think it checks out. (Where it came from, I have no idea.) I also think I believe the tensor amplification machinery. Parts I still don’t understand properly: 37 ˆ How the linear operation count and the synthetic roots of unity interact with the tensor amplification, and how to go from 2 × · · · × 2 to r × · · · × r. All of this feels like it should be okay, but I need to think about it some more to be certain.

ˆ How the rest of the 2021 algorithm behaves when d is allowed to grow at the required rate. Again this feels like it should be okay, but I haven’t checked.

ˆ The big one: data rearrangement costs in the Turing model. I haven’t even begun to look into this yet. Some final thoughts about the log n exponent. To improve it, one would need to find a more efficient gadget. I gather that this is what people (and their LLM helpers) have been doing. The impression I got from a very quick glance at their repositories is that they are trying to come up with clever optimisations of the existing construction. But one other thing that might be interesting to try (if they haven’t already done so) is to simply do some kind of (clever?) brute-force search for a very small gadget. Who knows, maybe it’s possible to write down a non-trivial gadget for I ⊗ B⊗3

acting on C4 ⊗ (C2)⊗3 or something silly like that. I mean, no-one has ever actually tried, right? Still a lot of work to do. But I need a break now.

References

[AFKL19] P. Afshani, C. B. Freksen, L. Kamma, and K. G. Larsen, Lower bounds for multiplica-tion via network coding , 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019) (Dagstuhl, Germany) (Christel Baier, Ioannis Chatzigian-nakis, Paola Flocchini, and Stefano Leonardi, eds.), Leibniz International Proceedings in Informatics (LIPIcs), vol. 132, Schloss Dagstuhl – Leibniz-Zentrum f¨ ur Informatik, 2019, pp. 10:1–10:12. [Har26] D. Harvey, Faster enumeration of primes , https://arxiv.org/abs/2606.22851 , 2026. [HvdH19] D. Harvey and J. van der Hoeven, Faster polynomial multiplication over finite fields using cyclotomic coefficient rings , J. Complexity 54 (2019), 101404, 18. MR 3983215 [HvdH21] , Integer multiplication in time O(n log n), Ann. of Math. (2) 193 (2021), no. 2, 563–617. MR 4224716 [HvdH22] , Polynomial multiplication over finite fields in time O(n log n), J. ACM 69

(2022), no. 2, Art. 12, 40. MR 4433320 [HvdH25] D. Harvey and J. van der Hoeven, Integer multiplication is at least as hard as matrix transposition , https://arxiv.org/abs/2503.22848 , 2025. [LL04] Zongpeng Li and Baochun Li, Network coding in undirected networks , CISS, 2004. [TC26a] The Company, Integer multiplication below n log n, GitHub, 2026. [TC26b] , Finite tensor savings and exact Fourier circuits , GitHub, 2026. [TC26c] , An explicit power saving for the exact discrete Fourier transform , GitHub, 2026. 38