The Grok

Multiplying is serious business. Well, mathematically speaking, it is the single most widely used operation in arithmetic. Yet, it is the naive method of taking individual digits of the multiplier and multiplying with the multiplicand the most popular method. Fortunately, this method yields good results for numbers which can be handled in normal data types. What if you had to multiply two very large numbers. With our naive method, a 1000 digit number would require 1000*1000 operation counts. This is clearly unacceptable as the operation count increases quickly with increase in number of digits.

One of the most fascinating application of the fast fourier transform is in this regard. Recall that an FFT can be performed in N*Log N steps. But how do we connect integer multiplication with something that transforms between domains? We also need to make sure the new operation we perform using the FFTed data needs to be at the most N*Log N steps, to successfully achieve fast multiplication.

The convolution operation comes to our rescue as it is nothing more than an elementwise multiplication of the FFTed data. That means, it will require a maximum of N steps. Thats not too shabby, as we managed to perform the conversion and the multiplication in O(n log n). The final step is to convert the data back to the original domain, which is nothing more than an inversion FFT.

So how does convolution achieve the same goal as our digit by digit multiplication?

FFT( f(x) * g(x) ) = F(s)G(s)
=> f(x) * g(x) = FFTinverse(F(s)G(s))

Now f(x) and g(x) are the original data set and F(s)G(s) is the element wise multiplication. Can you now see how the multiplication is achieved? f(x) and g(x) are polynomials. We simply call the digits the coefficients of x in increasing order.

Let us move on to the actual implementation. We shall reuse the FFT code we had earlier designed in our previous article.

Source Code
complex* fftmul(complex *x, complex *y, long n)
{
long pw=exp(2,n);
fft_convolution(x,y,n);

double cy=0.0;
for (long i=0; i<pw; i++)
{
double t=(y[i].real)+cy+0.5;
cy=(unsigned long) (t/1000.0);
y[i].real=t-cy*1000.0;
}
return y;
}

Explanation
The input data is specified in the form of the complex numbers (it is convenient as the FFTed data will be in complex form, so we can use the same structure without needing two separate structures). We store the entire number in the form of 3 digit sequences in little endian order. Suppose we wanted the number 23451234 then x[0]=234 , x[1]=451 and x[2]=023. The selection of the base/radix is very important, it gives us control over the accuracy and storage space. We have taken the base to be 1000 in this case,

There is one more small detail that has to be taken care of before the Inverse FFTed data can be used. The carries during multiplication have to be propagated which is fairly simply to code.

And that is all there is to turbo charging your multiplication. Go get 'em!

Download: Source