How is the inner ring chosen in the Schönhage–Strassen algorithm?

Question

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. $>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);
  ...

Answer 1

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 &= \lfloor s/2 \rfloor + 1 \\s_2 &= \lceil s/2 \rceil + 1 \\ p_2 &= \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) &(\leq)\; (p-1)2^sm + 2^m F(\lceil s/2\rceil+1) \\ &(\leq)\; 2s2^s (\lfloor s/2\rfloor+1) + 2^{\lfloor s/2\rfloor+1} F(\lceil s/2\rceil+1) \\ &(\leq)\; s^2 2^s + 2 \cdot 2^{s/2} F(s/2+1) \\ &(\leq)\; s^2 2^s + 4 (s/2)^2 2^s + 16(s/4)^2 2^s + ... \\ &(\leq)\; 2^s s^2 \lg(s) \\ &(\leq)\; \frac{n}{\lg n} \left(\lg \frac{n}{\lg n}\right)^2 \lg \lg \frac{n}{\lg n} \\ &(\leq)\; \frac{n}{\lg n} (\lg^2 n) \lg \lg n \\ &(\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.