# Schönhage-Strassen multiplication

I am trying to implement the Schönhage-Strassen algorithm (SSA) for multiplying large integers, but it only gives the right result if all $\delta_j$ are zero. I'll explain what I mean by this:

SSA (as described in the original 1971 paper Schnelle Multiplikation großer Zahlen) basically works like this:

1. Let $a$ and $b$ be two integers to multiply.
2. Determine $n$. Roughly speaking, $n$ will be about $\frac12 \log_2 \log_2 \max(a,b)$
3. Compute integers $z_j$ such that $a\cdot b=\sum\limits_j z_j\cdot2^{j\cdot2^{n-1}}$, as follows:
1. Compute $\xi_j\equiv z_j \pmod {2^{2^n}+1}$; this step performs an FFT+iFFT
2. Compute $\eta_j\equiv z_j \pmod {2^{n+2}}$
3. Compute $z_j$ from $\xi_j$ and $\eta_j$:
1. $\delta_j\equiv\eta_j-\xi_j \pmod {2^{n+2}}$ <--- THIS IS THE $\delta_j$ MENTIONED ABOVE
2. $z_j=\xi_j+\delta_j(2^{2^n}+1)$
4. Output $\sum\limits_j z_j\cdot2^{j\cdot2^{n-1}}$

I tested the other parts of the code, and I even did a a complete example with pencil and paper. Plus, everything works when $\delta=0$, so I'm pretty certain the problem is in step 3.3 or 4.

Wikipedia doesn't go into enough detail to even mention $\xi_j$ and $\eta_j$.

Does anyone know of a working example, sample code, or other information that might help me track down the problem?

Edit: I will put the code online if I get it to work.

Edit 2: I uploaded some Java code here. It runs as a regular Java program; the output shows one working example and one that doesn't work. Computation of $\delta,\xi,\eta,z_r$ happens in lines 148-160.

Also, there is a modified version here which uses much smaller input values, but I head to tweak the code so it handles pieces less than one byte.

-
There's extensive discussion of this algorithm, with code I believe, in Knuth's "The Art of Computer Programming" volume 2, in the section titled "How Fast Can We Multiply?" (and in the associated exercises). This section has changed quite a bit in various editions of the book, with the current edition using high precision floating point and careful error analysis to ensure accuracy. Earlier editions discussed the version that uses modular integer arithmetic so are probably worth tracking down. –  tzs Oct 16 '11 at 10:21

## 2 Answers

I finally found the problem. It was a silly bug at the end of the modFn method where $1$ is added if there is a carry left over, meaning the result is negative and $F_n$ needs to be added. When adding the carry produced another carry, it was being dropped instead of being propagated to the next digit.

The reason why it looked like something to do with $\delta$, $\xi$, etc. was because I had verified FFT, iFFT, and computation of $z_j \mod 2^{n+2}$ worked correctly. The only step left after that is to calculate $z_j$ from $z_j \mod 2^{n+2}$ and $z_j \mod F_n$, but I was forgetting about the modFn between FFT and iFFT which is what caused the problem.

-

I don't know how the core of SSA (=the FFT+iFFT part) works, but a possible source of an error occurred to me. If you have already eliminated this possibility, then I will delete this.

You really want to compute the integers $z_j$. It seems to me that the idea in step 3 is that you can somehow 'magically' (=the core of SSA) compute the remainder of the unknown $z_j$ w.r.t. to two moduli $M_1=2^{2^n}+1$ and $M_2=2^{n+2}$: in step 3.1. you produce an integer $\xi_j$ with the property $z_j\equiv \xi_j \pmod{M_1}$, and in step 3.2 an integer $\eta_j$ with the property $z_j\equiv \eta_j \pmod{M_2}$ is found. Then in step 3.3 you are to apply the Chinese remainder theorem to compute $z_j$ modulo the product modulus $M=M_1M_2$. This is then enough to determine the integer $z_j$ uniquely, because it is known to be in the interval $[0,M-1]$. If $n\ge2$, then formula 3.3.2 should work, because $2^n\ge n+2$, and therefore $2^{2^n}+1 \equiv 1 \pmod{M_2}$, so $$\xi_j+\delta_j(2^{2^n}+1)\equiv\xi_j+\delta_j=\xi_j+\eta_j-\xi_j=\eta_j\pmod{M_2}.$$ Note that obviously $$\xi_j+\delta_j M_1\equiv \xi_j \pmod{M_1}.$$

If you compute $\delta_j$ with a simple call to a standard integer division remainder operation, then many (if not all?) programming languages specify the remainder to have the same sign as the integer $\eta_j-\xi_j$ as opposed to being in the 'expected' range $[0,M_2-1]$. If you look at the line below equation 4.7. in Schönhage-Strassen paper, they specifically state that $\delta$ should be in this interval.

In order to check, whether this is the source of your problem, add a line that tests whether $\delta_j$ is negative, and adds $2^{n+2}$ to it, if that happens to be the case. If $n>2$, it may easily happen that $\eta_j-\xi_j<0$, because $M_1$ is then much bigger than $M_2$, so a 'random' remainder $\xi_j$ modulo $M_1$ is likely to be bigger than a 'random' remainder $\eta_j$ modulo $M_2$.

IOW, my theory is that the error may have resulted from mixing programmers' and mathematicians' use of 'mod'. Schönhage and Strassen are mathematicians.

-
Thanks, that was something I had not checked. Unfortunately, it wasn't the source of the problem. I appreciate you looking at the problem, though. –  Prashand Gupta Oct 16 '11 at 7:54
Ok. Then it's a mystery. My guess seemed to fit the symptoms in that your program was giving correct results, when all the $\delta_j$:s vanish. That was the reason I posted it. –  Jyrki Lahtonen Oct 16 '11 at 8:11