Sunday 21 April 2024

Fast (Fourier) multiplication

Figure 1: Representing the numbers 12 and 345 as frequency spectra,
 transforming them to time signals and multiplying.
One consequence of being able to convert between time signal and frequency spectrum representations is that they can be applied to the multiplication of multi-digit numbers in a way that is both straightforward and (potentially) quicker.

Consider the process of multiplying \(12\) by \(345\). I've fully expanded the intermediate steps with the colors representing the one's, ten's, hundred's and thousand's columns.
\[\begin{array}{r} 12\\ × \enspace 345\\ \hline \color{green}{10}\\ \color{red}{5}\enspace\\ _1 \color{red}{8}\enspace\\ \color{blue}{4}\enspace\enspace\\ _1 \color{blue}{6}\enspace\enspace\\ \color{goldenrod}{3}\enspace\enspace\enspace\\ \hline 4140 \end{array}\]
This operation is similarly represented in the following table where each number is separated into its digits and independently multiplied.

\( \begin{array}{|c|c|c|} \hline & & \color{gray}{10^1} & \color{gray}{10^0} \\ \hline & × & 1 & 2 \\ \hline \color{gray}{10^2} & 3 & \color{goldenrod}{3} & \color{blue}{6} \\ \hline \color{gray}{10^1} & 4 & \color{blue}{4} & \color{red}{8} \\ \hline \color{gray}{10^0} & 5 & \color{red}{5} & \color{green}{10} \\ \hline \end{array} \)

Calculating this, we get:

\(\begin{aligned}12 \cdot 345 & = (1  \!\cdot\!  {\color{gray}{10^1}} + 2  \!\cdot\!  {\color{gray}{10^0}})  \!\cdot\!  (3  \!\cdot\!  {\color{gray}{10^2}} + 4  \!\cdot\!  {\color{gray}{10^1}} + 5  \!\cdot\!  {\color{gray}{10^0}}) \\ & = 1  \!\cdot\!  {\color{gray}{10^1}}  \!\cdot\!  3  \!\cdot\!  {\color{gray}{10^2}} + 2  \!\cdot\!  {\color{gray}{10^0}}  \!\cdot\!  3  \!\cdot\!  {\color{gray}{10^2}} \, + \\ & \hspace{1.25em} 1  \!\cdot\!  {\color{gray}{10^1}}  \!\cdot\!  4  \!\cdot\!  {\color{gray}{10^1}} + 2  \!\cdot\!  {\color{gray}{10^0}}  \!\cdot\!  4  \!\cdot\!  {\color{gray}{10^1}} \, + \\ & \hspace{1.25em} 1  \!\cdot\!  {\color{gray}{10^1}}  \!\cdot\!  5  \!\cdot\!  {\color{gray}{10^0}} + 2  \!\cdot\!  {\color{gray}{10^0}}  \!\cdot\!  5  \!\cdot\!  {\color{gray}{10^0}} \\ & = {\color{goldenrod}{3}}  \!\cdot\!  {\color{gray}{10^3}} + {\color{blue}{6}}  \!\cdot\!  {\color{gray}{10^2}} + {\color{blue}{4}}  \!\cdot\!  {\color{gray}{10^2}} + {\color{red}{8}}  \!\cdot\!  {\color{gray}{10^1}} + {\color{red}{5}}  \!\cdot\!  {\color{gray}{10^1}} + {\color{green}{10}}  \!\cdot\!  {\color{gray}{10^0}} \\ & = {\color{goldenrod}{3}}  \!\cdot\!  {\color{gray}{10^3}} + ({\color{blue}{6}} + {\color{blue}{4}})  \!\cdot\!  {\color{gray}{10^2}} + ({\color{red}{8}} + {\color{red}{5}})  \!\cdot\!  {\color{gray}{10^1}} + {\color{green}{10}}  \!\cdot\!  {\color{gray}{10^0}} \\ & = {\color{goldenrod}{3}}  \!\cdot\!  {\color{gray}{10^3}} + {\color{blue}{10}}  \!\cdot\!  {\color{gray}{10^2}} + {\color{red}{13}}  \!\cdot\!  {\color{gray}{10^1}} + {\color{green}{10}}  \!\cdot\!  {\color{gray}{10^0}} \\ & = 4140 \end{aligned}\)

The coefficients \( [{\color{goldenrod}{3}}, {\color{blue}{10}}, {\color{red}{13}}, {\color{green}{10}}] \) are the convolution of \([1, 2]\) and \([3, 4, 5]\). Six component multiplications are required, thus performing in \( O(n^2) \) time.

A convolution can also be visualized as flipping the second number, sliding it under the first number, and accumulating products of corresponding digits along the way. For example, for \(12\) and \(345\):

  12    12     12    12
  ×     ××     ××     ×
543    543     543    543
-------------------------
  3     4      5      10
       + 6    + 8
-------------------------
  3     10     13     10
                
As before, six component multiplications are required.

However it so happens that if the numbers are each associated with a frequency spectrum (or, alternatively, a time signal), and transformed and multiplied in the dual domain then only four component multiplications are required, thus performing in \( O(n) \) time.[1]

Essentially, multiplying in the original domain (akin to standard multiplication) consists of multiplying digits and then adding the results. Whereas transforming and multiplying in the dual domain consists of adding digits and then multiplying the results.

As an example, consider how we might evaluate \((2 + 3) \cdot 7\). The multiplication-first approach is \(2 \cdot 7 + 3 \cdot 7 = 14 + 21 = 35\), requiring two multiplications. The multiplication-last approach is \((2 + 3) \cdot 7 = 5 \cdot 7 = 35\), requiring only one multiplication.

So how do we apply the multiplication-last approach for \(12 \cdot 345\)?

The first step is to represent each number as a frequency spectrum. The number \(12\) is represented as \([2, 1, 0, 0]\) and the number \(345\) is represented as \([5, 4, 3, 0]\). The digits are arranged in least-significant digit order, that is, by the one's digit first, then the ten's digit, and so on. They are also zero-padded to the output spectrum length which is equal to the total number of input digits minus one, that is, \(2 + 3 - 1 = 4\).

An inverse discrete Fourier transform is performed on each spectrum to produce their respective time signals:[2]

\(IDFT\,[2,1,0,0] = [3, 2+i, 1, 2-i]\)

\(IDFT\,[5,4,3,0] = [12, 2+4i, 4, 2-4i]\)

Next, each corresponding time sample is multiplied. This takes \(O(n)\) time.

\([3 \cdot 12, (2+i) \cdot (2+4i), 1 \cdot 4, (2-i) \cdot (2-4i)] = [36, 10i, 4, -10i]\)

A discrete Fourier transform is then performed to produce the final product spectrum (see Figure 1):[3]

\(DFT\,[36, 10i, 4, -10i] = [10, 13, 10, 3]\)

Flipping the order gives the coefficients \( [{\color{goldenrod}{3}}, {\color{blue}{10}}, {\color{red}{13}}, {\color{green}{10}}] \), therefore:

\( 12 \cdot 345 = {\color{goldenrod}{3}} \cdot {\color{gray}{10^3}} + {\color{blue}{10}} \cdot {\color{gray}{10^2}} + {\color{red}{13}} \cdot {\color{gray}{10^1}} + {\color{green}{10}} \cdot {\color{gray}{10^0}} = 4140 \)

The fastest algorithm for the IDFT/DFT is the fast Fourier transform which takes \(O(n\, log\, n)\) time and so the Fourier multiplication method improves on the conventional \(O(n^2)\) multiplication method for large numbers.

--

[1] The equivalence of convolution in the frequency domain and pointwise multiplication in the time domain (and vice-versa) is expressed by the convolution theorem.

Below is Python code for calculating the convolution and FFT-based multiplication for 12 and 345. Note that 1/4 scaling occurs in the forward transform (np.fft.fft).

import numpy as np

conv_result = np.convolve([1,2], [3,4,5])
print("Convolution:", conv_result.tolist())

fft_result = np.flip(np.fft.fft(np.fft.ifft(np.flip([0,0,1,2]), norm="forward") * np.fft.ifft(np.flip([0,3,4,5]), norm="forward"), norm="forward"))

print("FFT-based multiplication:", fft_result.real.astype(int).tolist())

Output:

Convolution: [3, 10, 13, 10]
FFT-based multiplication: [3, 10, 13, 10]

Note: if the numbers were represented as time signals rather than frequency spectra then the following line would be substituted (which switches the fft and ifft functions and defaults to norm="backward"). The output is identical.

fft_result = np.flip(np.fft.ifft(np.fft.fft(np.flip([0,0,1,2])) * np.fft.fft(np.flip([0,3,4,5]))))

[2] The math for the IDFT transforms is below. For each number, each digit is associated with a specific frequency over a time period - in this case, zero cycles for digits \(2\) and \(5\), one cycle for digits \(1\) and \(4\), and two cycles for digit \(3\). The resulting components of the IDFT transform are the samples at four regular points over the time period. To calculate the sample at each time point, each digit (as a complex amplitude) is rotated on the complex plane according to its associated frequency. In this case, the rotations are the fourth roots of unity: \(e^{2\pi ki/4}\) where \(k=0,1,2,3\) giving \(1, i, -1, -i\) (i.e., quarter turns from 1 on the complex unit circle).

\(\begin{aligned} IDFT\,[2,1,0,0] & = [2 \!\cdot\! 1 + 1 \!\cdot\! 1, \\ & \hspace{1.6em} 2 \!\cdot\! 1 + 1 \!\cdot\! i, \\ & \hspace{1.6em} 2 \!\cdot\! 1 + 1 \!\cdot\! -1, \\ & \hspace{1.6em} 2 \!\cdot\! 1 + 1 \!\cdot\! -i] \\ & \\ & = [2+1, \\ & \hspace{1.6em} 2+i, \\ & \hspace{1.6em} 2-1, \\ & \hspace{1.6em} 2-i] \\ & \\ & = [3, 2+i, 1, 2-i] \end{aligned}\)

For the initial time point no time has yet passed, so each digit remains unrotated (i.e., multiplied by \(1\)). At the second time point, digit \(2\) remains unrotated since its frequency is zero cycles and digit \(1\) is rotated by 90° (i.e., multiplied by \(i\)) since its frequency is one cycle. At the third time point, digit \(2\) remains unrotated and digit \(1\) is rotated by 180° (i.e., multiplied by \(-1\)). At the fourth time point, digit \(2\) remains unrotated and digit \(1\) is rotated by 270° (i.e., multiplied by \(-i\)).

\(\begin{aligned} IDFT\,[5,4,3,0] & = [5 \!\cdot\! 1 + 4 \!\cdot\! 1 + 3 \!\cdot\! 1, \\ & \hspace{1.6em} 5 \!\cdot\! 1 + 4 \!\cdot\! i + 3 \!\cdot\! -1, \\ & \hspace{1.6em} 5 \!\cdot\! 1 + 4 \!\cdot\! -1 + 3 \!\cdot\! 1, \\ & \hspace{1.6em} 5 \!\cdot\! 1 + 4 \!\cdot\! -i+3 \!\cdot\! -1] \\ & \\ & = [5+4+3, \\ & \hspace{1.6em} 5+4i-3, \\ & \hspace{1.6em} 5-4+3, \\ & \hspace{1.6em} 5-4i-3] \\ & \\ & = [12, 2+4i, 4, 2-4i] \end{aligned}\)

Similarly, for the initial time point no time has yet passed, so each digit remains unrotated (i.e., multiplied by \(1\)). At the second time point, digit \(5\) remains unrotated since its frequency is zero cycles, digit \(4\) is rotated by 90° (i.e., multiplied by \(i\)) since its frequency is one cycle and digit \(3\) is rotated by 180° (i.e., multiplied by \(-1\)) since its frequency is two cycles. At the third time point, each digit is rotated by 0°, 180° and 360° (or 0°) respectively. At the fourth time point, each digit is rotated by 0°, 270° (or -90°) and 540° (or 180°) respectively.

[3] The expanded math for the transforms and multiply operation.

\(DFT\,(IDFT\,[2,1,0,0] \cdot IDFT\,[5,4,3,0])\)

\(\begin{aligned}\hspace{0.8em} = \, DFT\,[ & (2 \!\cdot\! 1 + 1 \!\cdot\!  1)  \!\cdot\!   (5 \!\cdot\!  1 + 4 \!\cdot\!  1 + 3 \!\cdot\!  1), \\ & (2 \!\cdot\!  1 + 1 \!\cdot\!  i)  \!\cdot\!   (5 \!\cdot\!  1 + 4 \!\cdot\!  i + 3 \!\cdot\!  -1), \\ & (2 \!\cdot\!  1 + 1 \!\cdot\!  -1)  \!\cdot\!   (5 \!\cdot\!  1 + 4 \!\cdot\!  -1 + 3 \!\cdot\!  1), \\ & (2 \!\cdot\!  1 + 1 \!\cdot\!  -i)  \!\cdot\!   (5 \!\cdot\!  1 + 4 \!\cdot\!  -i + 3 \!\cdot\!  -1)] \\ & \\ =  \, DFT\,[ & (2 + 1) \!\cdot\! (5 + 4 + 3),\\ &  (2 + i) \!\cdot\! (5 + 4i -3),\\ &  (2 - 1) \!\cdot\! (5 - 4 + 3),\\ &  (2 - i) \!\cdot\! (5 - 4i - 3)]\end{aligned}\)

\(\begin{aligned}\hspace{0.6em} & = \, DFT\,[3 \cdot 12, (2+i) \cdot (2+4i), 1 \cdot 4, (2-i) \cdot (2-4i)] \\ & \\ & = \, DFT\,[36, 10i, 4, -10i]\end{aligned}\)

\(\begin{aligned}\hspace{0.8em} = \tfrac{1}{4} [ & 36 \!\cdot\! 1 + 10i \!\cdot\! 1 + 4 \!\cdot\! 1 - 10i \!\cdot\! 1, \\ & 36 \!\cdot\! 1 + 10i \!\cdot\! -i + 4 \!\cdot\! -1 - 10i \!\cdot\! i, \\ & 36 \!\cdot\! 1 + 10i \!\cdot\! -1 + 4 \!\cdot\! 1 - 10i \!\cdot\! -1, \\ & 36 \!\cdot\! 1 + 10i \!\cdot\! i + 4 \!\cdot\! -1 - 10i \!\cdot\! -i] \end{aligned}\)

\(\begin{aligned}\hspace{0.8em} = \tfrac{1}{4} [ & 36 + 10i + 4 - 10i, \\ & 36 + 10 - 4 + 10, \\ & 36 - 10i + 4 + 10i, \\ & 36 - 10 - 4 - 10] \end{aligned}\)

\(\begin{aligned}\hspace{0.8em} & = \tfrac{1}{4} [40, 52, 40, 12] \\ & \\ & = [10, 13, 10, 3]\end{aligned}\)