New Integer Multiplication Algorithm Claims O(n log n) Complexity

On the recently proposed o(nlogn) 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(nlogn) bit operations.

1 Executive summary

I have not checked most of the details, but I think that the claimed o(nlogn) 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 algorithms. 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 nonmathematical 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(nlogn) 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(logn)1−ε) for some small positive constant ε.

This bound is of course utterly ridiculous. Previously the best bound was O(nlogn) (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.

I 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 hardworking 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(nlogn) integer multiplication result. The new work builds directly on top of the O(nlogn) 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 (logn)ε improvement is over the previous O(nlogn) record, harping especially on the tiny numerical value of ε. In the initial writeup ε is 2−182, meaning that the exponent of logn 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(nlog n · Klog*n) bound in 2007, the improvement that he had achieved over the previous state of the art, namely the epochal Schönhage–Strassen O(nlognlog logn) bound from 1971, was smaller than log logn. And the subsequent improvement made by my work with van der Hoeven (and several other papers along the way) was much, much, much smaller than this. The claimed (logn)ε 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(logn)α) 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(nlogn) 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 that 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, assuming 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(logn)α). We previously believed that α = 1, but the Model’s claim is that α 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(nlogn). 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 2r-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(logn) bits of precision.

The number of butterflies (over C) is easily calculated to be 1/2rd log2(rd−1). One can already see the O(nlogn) 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) 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(nlogn) 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(nlogn).

If you had asked me two weeks ago for my subjective estimate of the probability that either of these O(nlogn) 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 1 / 1 -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 ≅ C^(2^m) as an m-dimensional array of complex numbers of size 2 ×···× 2 (m factors).

If A1,...,Am 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 · 2^(m−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 · 2^(m−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 · 2^(m−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 · 2^(m−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 · 2^(n−1) to (t/(m · 2^(m−1))) · n · 2^(n−1). This saving is not very interesting, because it is by a factor of only t/(m · 2^(m−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 ≅ C^(2^m), i.e., the vector u = (u1,...,um) ∈ F2m refers to the (u1,...,um)-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 ∈ GLm(F2), the operator PM that acts on (C2)⊗m by permuting the entries according to x 7→ Mx (for x ∈ F2m).

For example, if M is a permutation matrix that swaps two coordinates of F2m, then PM just switches two axes of the array. But M could also involve more complicated transformations that “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 PM0Ai1PM1Ai2··· AigPMg. (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 = (PM0Ai1PM1Ai2··· AigPMg)⊗ℓ.

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 = m^k at the beginning, then overall we save a factor of (g/m)^k. This is the amplification: the additional structure lets us save further factors of g/m as n grows. The total number of applications of A becomes (g/m)^k · n · 2^(n−1) = exp((log n / log m) log(g/m)) · n · 2^(n−1) = n^α · 2^(n−1), where α = log g/ logm. You can see here the N(logN)^α 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^(m^k), 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 1 / 1 -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(logN)^α) 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 = m^k 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 Wm, 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 multiplication; 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 Ω(n^2) 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 n^3, i.e., o(n^2) operations on average per vector. So I do not see any obvious obstruction here.

7 The gadget

Let B be the butterfly map (1 1 / 1 -1) acting on C2. The goal in this section is to explain the construction of positive integers m and W and a “program” that computes the action of I⊗B⊗m on CW ⊗(C2)⊗m in fewer than Wm 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 1 / 1 -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 b / b a), where 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 0 / 0 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, 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 √2 which 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.

(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 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 F2m → C, where F2m are the addresses as before.

Let v ∈ F2m, 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 ∈ F2m. Then define Cv := aI + bXv.

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 ∈ F2m. Here v · x is the standard bilinear form on F2m, i.e., v1x1 + ··· + vmxm ∈ F2.

For Cv we just use the formula Cv = aI + bXv. 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]), 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 F2m = {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(i^wt(x)), 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,...,Cem on each array; this amounts to Wm 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 F2m

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

Let U ⊆ F2m be a subspace. We write U⊥ ⊆ F2m for the complement with respect to the standard dot product, i.e., U⊥ := {v ∈ F2m: u · v = 0 for all u ∈ U}. The rank-nullity theorem shows that always dim(U⊥) = m − dimU.

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⊥ = F2m. 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,...,ej} is nondegenerate for any j.

When U is nondegenerate, it makes sense to define a projection operator PU : F2m → U. We just write a given x ∈ F2m uniquely as y + z where y ∈ U and z ∈ U⊥, and then put PU(x) := y.

Now for any nondegenerate subspace U, define qU: F2m → Z/4Z by qU(x) := wt(PUx) mod 4.

For example, q0 is just the zero function, and qF2m is the function that returns wt(x) (mod 4), i.e., exactly the exponents appearing in (7.2). Then define ΦU := H⊗m · diag(i^qU(x)) · (H⊗m)−1, so Φ0 = I, and ΦF2m = 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 F2m. You can prove this by observing that wt(x + y) = wt(x) + wt(y) − 2*|suppx ∩ suppy| (the correction term comes from the coordinates where both xi = 1 and yi = 1), and if x · y = 0 then |suppx ∩ suppy| is even.

It follows that that qU+V = qU + qV, since qU+V(x) = wt(PUx + PVx) mod 4 = qU(x) + qV(x), using the fact that PUx ⊥ PVx. 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 ΦF2m = 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.

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 if wt(v) ≡ 1 (mod 4), and (Cv)−1 if 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,...,vk of odd weight, each one orthogonal to U and all the previous vectors, then U + ⟨v1,...,vk⟩ is also nondegenerate, and ΦU+⟨v1,...,vk⟩ is simply the product ΦU Cv1±1 ··· Cvk±1 with appropriately chosen signs (depending on the weights of the vi modulo 4). Each of the Cvi±1 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 F23 is the standard one. It suggests instead an m = 4 example with orthonormal basis 1110,1101,1011,0111. Then one has C⊗4 = C1110−1 C1101−1 C1011−1 C0111−1.

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,...,rW). Each register holds an element of (C2)⊗m. Here are the allowed operations:

• We can apply Cv or Cv−1 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 + 1/2r1. As (vaguely) explained above, I am ignoring the cost of these operations.

Each register has an associated (nondegenerate) subspace Ui. Initially the subspaces U1,...,UW 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 ΦUifi. To make this work, we have to follow some rules:

• Suppose that we want to enlarge Ui to Ui′:= Ui + ⟨v⟩. For this to be allowed, we require that v has odd weight and that v ⊥ Ui. According to the previous section, ΦU is ΦUiCv if wt(v) = 1 mod 4, and ΦUiCv−1 if wt(v) = 3 mod 4. So in the first case we apply Cv to ri, in the second case we apply Cv−1. In either case, the value in the register is still ΦUifi 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 Ui′:= Ui ∩⟨v⟩⊥. (This is equivalent to enlarging Ui⊥ to Ui⊥ + ⟨v⟩.) For this to be allowed, we require that v has odd weight and that v ∈ Ui⊥. Similarly to the previous bullet, if wt(v) = 1 mod 4 then we apply Cv−1 to ri, if wt(v) = 3 mod 4 we apply Cv. Again the underlying value fi does not change.

• 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 + 1/2ri. In this example we have effectively replaced fj with fj + 1/2fi. Everything stays consistent because ΦU (fj + 1/2fi) = ΦU fj + 1/2Φ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,...,Cem to the register, then Ui gets updated to F2m, so that ΦUi= C⊗m. At this point, the underlying value is still fi, but the value physically in the register is ΦUifi = C⊗m fi. So we have succeeded in computing C⊗m fi.

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 choose 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-3 choose 3) choices of T with |S ∩ T | = 0, so altogether |P| = (h choose 3)(h-3 choose 3) + 3(h-3).

Initialisation

There are many registers, playing various different roles. Instead of calling them r1 ,...,r W , we use:

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

It is perhaps worth noting that if you count the number of registers in each category, then asymptotically 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 ∈ F2h 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 Φ⟨vT⟩−1 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 1/2AS,T from yS, and if |S ∩ T | = 2 then add 1/2AS,T to yS.

Cost: zero.

Final state:

register	current U	underlying value
xT	⟨vT⟩	˜xT
yS	0	˜yS − 1/2 ΣT:
AS,T	0	˜AS,T
ci	0	˜ci
c*	0	˜c*