Based on skimming the paper (finally a proper explanation that gives some hint at intuition), the path to openai's formulation if you're familiar with undergrad CS is (apparently):

* Start with the usual FFT convolution for polynomials, viewing integers as a polynomial in a given base evaluated at a certain point

* In undergrad you (hopefully) saw polynomial multiplication can be done in nlogn via convolution using FFT, but that assumes complex integer addition/multiplication is O(1). Here we need O(log x) bits to simply represent a number x on a real machine.

* If you use naive (n^2) multiplication for those butterfies, then you end up with O(n[(log n)^2]).

* If you recursively use FFT for the digit multiplication at those butterflies, then you get something like O(n[(log n)^(log log n)]).

* Schonhage-Strassen realized that instead of doing FFT in C where the butterfly involves multiplying by \omega (\omega being roots of unity), you can work in a different ring where 2 is a root of unity and so it's a cheap bitshift, but were constrained in how finely they could break subproblems (could only break an n-bit integer into sqrt-n subpproblems each of sqrt-n bits), resulting in an overall O(nlogn * loglogn).

* The 2019 O(nlogn) version (supposedly, this is out of my depth) found a way to use some technique of "Gaussian Resampling" to do the matrix convolution via _multidimensional FFT_ in a ring where multiplying by root of unity became cheap _and_ subproblems could be finer. The bottleneck was thus only the addition/subtraction in the O(log n) butterfly layers, each of which does O(n) work for total of O(n logn)

* OpenAI's paper uses a fact that it's apparently easier to compute multiple butterfly layers "simultaneously" than doing them one at a time independently (page 10). There is some sort of implicit reuse because the later layers mix the outputs of earlier layers, and they take advantage of this to avoid having to write out layer 1 outputs separately from layer 2 outputs, and instead work in some coordinate system where they can intermingle them. This apparently amortizes to some tiny fractional epsilon shaved off of a layer, and since the whole FFT is recursive they can apparently save a factor of (1 - epislon)^(log log n) overall across "log log n" recursion depth, which is apparently where that ridiculously small overall exponent comes from.