Multiplying 2 positive integers using FFT and convolutions

Question

I was trying to figure out how I can perform multiplication of 2 big integers using FFT and convolutions, I ran into the following article: http://numbers.computation.free.fr/Constants/Algorithms/fft.html
I understand the main points behind the algorithm but I can't understand this part:

2.3  More on the complexity of multiplication with FFT

In fact, the time complexity of multiplication with FFT is a little bigger than n log(n).

Let us be more precise.

To multiply two numbers of N digits, we write them in a base B which contains k digits (say B = $10^k$), thus giving a number of coefficients equal to n $\approx$ N/k.

The discussion above tells us that to multiply those two numbers, FFT permits to perform $O(n log(n))$ operations on basic numbers (basic numbers express coefficients in the base B, they are usually basic numerical data types like double in C of Fortran)

Because of the numerical error bound (1), these basic numbers should be precise enough to represent integers up to 6n2B2log(n) (for example, working in double precision, the base B should be choosen small enough so that the error bound is not too large... and this is not even possible if n is too large).

Thus the number of digits of these basic numbers should be of the order of log(B)+log(n). As a consequence, the basic operations on these numbers has cost O((log(B)+log(n))2) and the final cost is O( n log(n) (log(B) + log(n))2). The base B is choosen so that $k = log_{10}(B)$ is of the same order than log(n), and finally, multiplying those two numbers of $N \approx kn$ digits have cost $O( n log(n)^3 ) = O( N log(N)^2)$.

I need some clarification on how to write the numbers in a certain base B and why when $B = 10^k$ the base contains k digits in their example, it shouldn't contain $10^k$ digits? Such implementation should be recursive? I'll be grateful if someone can provide an example with numbers on how to execute this algorithm (without performing the FFT of course). thanks!

Answer 1

Kind of covered this in a class of mine so taking excerpts from lecture slides/online resources here is my attempt at answering it.

In the context of using FFT for multiplication of two large integers, writing the numbers in a certain base $B$ refers to expressing the coefficients of the numbers in that base. This is done in order to reduce the number of coefficients that need to be multiplied, which in turn reduces the complexity of the multiplication. For example, if you have two large integers $X$ and $Y$ that you want to multiply, you can express each of these numbers in base $B$ as:

$\begin{align} X &= x_0 + x_1 \cdot B + x_2 \cdot B^2 + \cdots + x_n \cdot B^n\\ Y &= y_0 + y_1 \cdot B + y_2 \cdot B^2 + \cdots + y_n \cdot B^n \end{align}$

where $x_0, x_1, ..., x_n$ and $y_0, y_1, ..., y_n$ are the coefficients of the numbers in base $B$. The idea is to choose the base $B$ so that the number of coefficients $n$ is significantly smaller than the number of digits in the original numbers ($N$), allowing you to perform the multiplication more efficiently.

For example, if you choose $B = 10^k$, where $k$ is the number of digits in a single coefficient, then the base $B$ will contain $k$ digits. For instance, if you have two numbers $X$ and $Y$ with a total of 1000 digits, and you choose $k = 3$, then you can express $X$ and $Y$ in base $B = 1000$, which will contain 3 digits in each coefficient. This will result in a total of $n = 1000/3 = 333$ coefficients for each number, which is significantly smaller than the original number of digits ($N = 1000$).

The reason why the base $B$ is chosen so that $k = \log_{10}B$ is of the same order as $\log n$ is because this allows the number of coefficients $n$ to be significantly smaller than the number of digits in the original numbers ($N$), which in turn reduces the complexity of the multiplication. For example, if you have two numbers $X$ and $Y$ with 1000 digits each, and you choose $k = 3$, then you will have a total of $n = 1000/3 = 333$ coefficients for each number. This is significantly smaller than the original number of digits ($N = 1000$), which allows you to perform the multiplication more efficiently.

I also don't think this approach is recursive, as the FFT algorithm itself is recursive. The idea is to express the numbers in a certain base $B$ in order to reduce the number of coefficients that need to be multiplied, which in turn reduces the complexity of the multiplication.

Answer 2

The Algorithm Archive has an article on Multiplication as a convolution where they give enough details to understand how to implement it, and they also have example code on both convolutions and FFT you can use to implement a full multiplication algorithm.

I implemented a version for my student to illustrate, perhaps it becomes clearer when you see actual code, and when you can experiment with it yourself:

from sys import argv


def convolve(F, G):
    NF = len(F)
    NG = len(G)
    N = NF + NG - 1
    C = [0 for _ in range(N)]
    s = 0
    for i in range(N):
        for j in range(max(0, i - NG), i + 1):
            if j < NF and (i - j) < NG:
                s += F[j] * G[i - j]
        C[i] = s
        s = 0
    return C


def carry(data):
    """Perform carrying, meaning that
    the "number" [2, 17, 5] becomes [3, 7, 5].
    """
    change = 1
    while change:
        change = 0
        cp = [0 for _ in data]
        for idx, e in enumerate(data):
            if e >= 10:
                change = 1
                car = e // 10
                cp[idx - 1] += car
                cp[idx] = e - (car * 10)
            else:
                cp[idx] = e
        data = cp
    return cp


if __name__ == "__main__":
    if len(argv) < 2:
        exit("Usage: convolution A B (integers)")
    x, y = [int(e) for e in argv[1]], [int(e) for e in argv[2]]
    print(f"{x=}\n{y=}")
    C = convolve(x, y)
    print(f"{C=}")
    M = carry(C)
    strM = "".join([str(e) for e in M])
    print(f"{argv[1]} * {argv[2]} = {strM}")

When run from the terminal: (the C array is the carrying that needs to be done)

[~]$ python convolution.py 1234 5678
x=[1, 2, 3, 4]
y=[5, 6, 7, 8]
C=[5, 16, 34, 60, 61, 52, 32]
1234 * 5678 = 7006652

Answer 3

Multiplication using FFT and convolution requires that your arithmetic is performed with some precision. There's always an argument "I use floating point arithmetic with rounding errors, but since I know that the correct result must be an integer, as long as I can prove that the rounding error is less than 0.5, I can determine the correct integer result." For example, if a result is known to be an integer for obvious mathematical reasons, and I tell you "I found the result 119.46, and the rounding error is less than 0.5", then it is obvious that the correct result is 119 and not 120 or any other number.

If your FFT uses many numbers, rounding errors become bigger. At some point using double precision floating point doesn't allow you to give a guarantee that the error is less than 0.5. So you use extended precision, or quad precision. So you can handle much much larger numbers with an error less than 0.5. And at some point that isn't enough, so you need to implement higher precision.

All in all, the number of operations grows with n log n, but the cost of each operation also grows very very slowly with n. Usually we will just implement say quad precision and note that "this implementation will not give the correct result if n > (some huge number)". And that will be acceptable if (some huge number) is so big that no computer in the world would be able to calculate an FFT of that size anyway. Except in theory, it is not acceptable.

Now about recursion: Your double precision, quad precision, six time precision floating point arithmetic would be hand coded. However, n could become so incredibly huge that the basic operations are big enough to make implementing them with FFT becomes worthwhile. You could have an astronomically huge n so that your basic arithmetic needs to be performed with 10,000 bit precision, and then you would implement your basic arithmetic with FFT. And then we could make n even more astronomically huge, so that the FFT we use for implementing its basic arithmetic is so huge that we need higher precision for implementing that.

So to do things in a CS correct way, you figure out what number of bits you need to do FFTs of n numbers correctly, and what time you would need for the basic operations implementing your FFTs, and take that into account in your big-O.

There is also a psychological effect. Multiplication of two n bit numbers in a straightforward way takes O(n^2) operations. But a typical computer today implements this for n ≤ 64 using O(n^2) hardware and constant time. So in practice this O(n^2) factor doesn't matter to us and we ignore it. 3 x 3 bit or 64 x 64 bit are the same on a modern computer. (On an ancient 8 bit computer, the first would take one and the second would take 64 multiplications if implemented natively, so 40 years ago I wouldn't have said 64 x 64 bit takes constant time).