How is the inner ring chosen in the Schönhage–Strassen algorithm? I've been trying to implement the Schönhage–Strassen integer multiplication algorithm, but hit a stumbling block in the recursive step.

I have a value $x$ with $n$ bits and I want to compute $x^2 \pmod {2^n+1}$. I originally thought the idea was to pick a $k$ such that $4^k \geq 2n$, split $x$ into $2^k$ pieces each with $2^{k-1}$ bits, apply SSA's convolution while working modulo $2^{2^k}+1$, a ring with $2^k$ bits of capacity per value, then put the pieces back together. However, the convolution's output has slightly more than $2n$ bits (i.e. $&gt;2^k$ bits per output value, which is more than the capacity of the ring, due to each output value being a sum of several products) so this doesn't work. I had to add on an extra factor of 2 of padding.

That extra factor of 2 in the padding ruins the complexity. It makes my recursive step too expensive. Instead of an $F(n) = n \lg n + \sqrt{n} F(2 \sqrt{n}) = \Theta(n \; \lg n \; \lg \lg n)$ algorithm, I end up with an $F(n) = n \lg n + \sqrt{n} F(4 \sqrt{n}) = \Theta(n \lg^2 n)$ algorithm.

I read a few references linked from wikipedia, but they all seem to gloss over the details of how this issue is solved. For example, I could avoid the extra padding overhead by working modulo $2^{p 2^k} + 1$ for a $p$ that's not a power of 2... but then things just break later, when I have only non-power-of-2 factors remaining and can't apply Cooley-Tukey without doubling the number of pieces. Also, $p$ may not have a multiplicative inverse modulo $2^p+1$. So there's still forced factors of 2 being introduced.

How do I pick the ring to use during the recursive step, without blowing the asymptotic complexity?

Or, in pseudo code form:



multiply_in_ring(a, b, n):
  ...
  // vvv                          vvv //
  // vvv HOW DOES THIS PART WORK? vvv //
  // vvv                          vvv //
  let inner_ring = convolution_ring_for_values_of_size(n);
  // ^^^                          ^^^ //
  // ^^^ HOW DOES THIS PART WORK? ^^^ //
  // ^^^                          ^^^ //

  let input_bits_per_piece = ceil(n / inner_ring.order);
  let piecesA = a.splitIntoNPiecesOfSize(inner_ring.order, input_bits_per_piece);
  let piecesB = b.splitIntoNPiecesOfSize(inner_ring.order, input_bits_per_piece);

  let piecesC = inner_ring.negacyclic_convolution(piecesA, piecesB);
  ...
This answer is taken from the paper "Asymptotically fast algorithms for the numerical muitiplication and division of polynomials with complex coefficients" that Markus linked in the comments.



You want to square an $n$-bit number, modulo $2^n + 1$. Here's what you do:


Find $p$ and $s$ that satisfy $n = (p-1) 2^s$ and $s \leq p \leq 2s$.
Pick the number of pieces $2^m$ to split the $n$ bits into, and the corresponding parameters for the piece sizes:

$$\begin{align} m &amp;= \lfloor s/2 \rfloor + 1 \\s_2 &amp;= \lceil s/2 \rceil + 1 \\ p_2 &amp;= \lceil p/2 \rceil + 1 \end{align}$$

Note that $s_2$ and $p_2$ continue to satisfy the $s_2 \leq p_2 \leq 2 s_2$ invariant. Also note that $2^m 2^{s_2} p_2 \geq 2n + m + 1$ is satisfied, so the input fits with room for carries.
Perform the FFT-based negacyclic convolution on the pieces, and the rest, as usual.


So that's the overarching idea: a logarithmic padding factor $p$. Now for the complexity analysis. The FFT will take $n m$ work to do, and we're recursing on $2^m$ pieces of size $(p_2-1) 2^{s_2}$, so now we can do extremely rough math with the recurrence relation w.r.t. $s$:

$$\begin{align}
F(s)
&amp;(\leq)\; (p-1)2^sm + 2^m F(\lceil s/2\rceil+1)
\\
&amp;(\leq)\; 2s2^s (\lfloor s/2\rfloor+1) + 2^{\lfloor s/2\rfloor+1} F(\lceil s/2\rceil+1)
\\
&amp;(\leq)\; s^2 2^s + 2 \cdot 2^{s/2} F(s/2+1)
\\
&amp;(\leq)\; s^2 2^s + 4 (s/2)^2 2^s + 16(s/4)^2 2^s + ...
\\
&amp;(\leq)\; 2^s s^2 \lg(s)
\\
&amp;(\leq)\; \frac{n}{\lg n} \left(\lg \frac{n}{\lg n}\right)^2  \lg \lg \frac{n}{\lg n}
\\
&amp;(\leq)\; \frac{n}{\lg n} (\lg^2 n) \lg \lg n
\\
&amp;(\leq)\; n \;(\lg n) \lg \lg n
\end{align}$$

Which seems about right, although I cheated quite a lot in those steps.

The 'trick' seems to be the fact that we end up with $s^2$ instead of $s$ in the base cost. There's still two multiplications by two per recursive level, like I was complaining about in the question, but now the halving of $s$ is paying double-dividends so it all works out. Then, at the end, we cancel the extra factor of $s$ (which is actually a factor of $\log n$) thanks to making $p$ logarithmically large relative to $s$ initially.