<p><em>This answer is taken from the paper <a href="http://link.springer.com/chapter/10.1007%2F3-540-11607-9_1" rel="nofollow">"Asymptotically fast algorithms for the numerical muitiplication and division of polynomials with complex coefficients"</a> that Markus linked in the comments.</em></p>

<hr>

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

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

<p>$$\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}$$</p>

<p>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.</p></li>
<li><p>Perform the FFT-based negacyclic convolution on the pieces, and the rest, as usual.</p></li>
</ul>

<p>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$:</p>

<p>$$\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}$$</p>

<p>Which seems about right, although I cheated quite a lot in those steps.</p>

<p>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.</p>

CC BY-SA 3.0
2016-07-24T03:57:19.590Z
2016-07-24T04:25:40.103Z
2016-07-24T04:25:40.103Z
Craig Gidney
<p>I've been trying to implement the Schönhage–Strassen integer multiplication algorithm, but hit a stumbling block in the recursive step.</p>

<p>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.</p>

<p>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.</p>

<p>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.</p>

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

<p>Or, in pseudo code form:</p>



<pre><code>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);
  ...
</code></pre>

CC BY-SA 3.0
2016-07-19T09:19:45.680Z
Please [do not post the same question on multiple sites](http://meta.stackexchange.com/q/64068). Each community should have an honest shot at answering without anybody's time being wasted.   I suggest you delete one of the two copies.
D.W.
CC BY-SA 3.0
2016-07-19T16:04:14.740Z
@D.W. Done. I cross-posted after cs didn't give any answers for a week, figuring it was too hard for that site. Was going to link back any answers obviously.
Craig Gidney
CC BY-SA 3.0
2016-07-19T21:58:02.220Z
I understand.  If it comes up in the future, you can always flag your post for moderator attention and ask to have it migrated, and we can move it for you over to CSTheory.  Thank you for your understanding!
D.W.
CC BY-SA 3.0
2016-07-20T09:53:59.257Z
There is a version of the algorithm that works modulo numbers of the form $2^{\nu2^n}$: A. Schönhage. Asymptotically fast algorithms for the numerical multiplication and division of polynomials with complex coefficients.
In EUROCAM '82: European Computer Algebra Conference, Lect. Notes Comp. Sci. 144, 3-15. http://www.iai.uni-bonn.de/~schoe/publi39.dvi
Markus Bläser
CC BY-SA 3.0
2016-07-20T13:29:35.370Z
IIRC you had a partial self-answer on the now deleted CS question. It seems a shame to lose that. Could you include it here (in the question, so that the question isn't marked as already answered)?
Peter Taylor
CC BY-SA 3.0
2016-07-20T14:19:37.353Z
@PeterTaylor That partial answer basically amounted to "leave gaps amongst the bits aligned on a power of 2 boundary that you treat as the principal root to absorb the carries". In other words do polynomial convolution with $x^n\equiv-1$, because we have a bunch of chunks (coefficients) that we don't want to interact. The problem is that once you get down to the granularity of the chunks you have to scale out by an extra factor of 2 before continuing. So you pay x2 at $\lg n$ and $\lg \lg n$ and $\lg \lg \lg n$ and ... and suddenly the algorithm has a $2^{\lg^* n}$ iterated logarithm problem.
Craig Gidney
CC BY-SA 3.0
2016-07-21T02:55:12.537Z
@MarkusBläser Thanks, that does go into the exact details of how to pick the ring. I'll study it, prove it for myself, and write up an answer.
Craig Gidney
<ds.algorithms>
How is the inner ring chosen in the Schönhage–Strassen algorithm?
Craig Gidney
# How is the inner ring chosen in the Schönhage–Strassen algorithm?

(Asked by: username_0 on 2016-07-13)

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

**Tags:** ds.algorithms

> **Comments:**
> **username_1**: Please [do not post the same question on multiple sites](http://meta.stackexchange.com/q/64068). Each community should have an honest shot at answering without anybody's time being wasted.   I suggest you delete one of the two copies.  
> **username_0**: @username_1 Done. I cross-posted after cs didn't give any answers for a week, figuring it was too hard for that site. Was going to link back any answers obviously.  
> **username_1**: I understand.  If it comes up in the future, you can always flag your post for moderator attention and ask to have it migrated, and we can move it for you over to CSTheory.  Thank you for your understanding!  
> **username_2**: There is a version of the algorithm that works modulo numbers of the form $2^{\nu2^n}$: A. Schönhage. Asymptotically fast algorithms for the numerical multiplication and division of polynomials with complex coefficients. In EUROCAM '82: European Computer Algebra Conference, Lect. Notes Comp. Sci. 144, 3-15. http://www.iai.uni-bonn.de/~schoe/publi39.dvi  
> **username_3**: IIRC you had a partial self-answer on the now deleted CS question. It seems a shame to lose that. Could you include it here (in the question, so that the question isn't marked as already answered)?  
> **username_0**: @PeterTaylor That partial answer basically amounted to "leave gaps amongst the bits aligned on a power of 2 boundary that you treat as the principal root to absorb the carries". In other words do polynomial convolution with $x^n\equiv-1$, because we have a bunch of chunks (coefficients) that we don't want to interact. The problem is that once you get down to the granularity of the chunks you have to scale out by an extra factor of 2 before continuing. So you pay x2 at $\lg n$ and $\lg \lg n$ and $\lg \lg \lg n$ and ... and suddenly the algorithm has a $2^{\lg^* n}$ iterated logarithm problem.  
> **username_0**: @MarkusBläser Thanks, that does go into the exact details of how to pick the ring. I'll study it, prove it for myself, and write up an answer.  

---

## 1 answers

### Answer by username_0 on 2016-07-24. Score: 4 | Accepted answer

*This answer is taken from the paper ["Asymptotically fast algorithms for the numerical muitiplication and division of polynomials with complex coefficients"](http://link.springer.com/chapter/10.1007%2F3-540-11607-9_1) 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.


---

