Use the fast Fourier transform to calculate the big integer multiplication

We know that the multiplication of two n-digit integers, if the regularAlgorithmThe time complexity is O (n2 ). However, using the fast Fourier transformation, the time complexity can be reduced to O (n logn ).

 

Assume that we want to calculate the product of the following two n-digit numbers:

A = (aN-1aN-2... a1a0) 10 = aN-1x10N-1 + aN-2x10N-2 +... + a1x101 + a0x100

B = (bN-1bN-2... b1b0) 10 = bN-1x10N-1 + bN-2x10N-2 +... + b1x101 + b0x100

Multiply the two formulas to obtain the following formula (2n-1 in total ):

C = A x B = c2N-2x102N-2 + c2N-3x102N-3 +... + c1x101 + c0x100

It is very easy to verify that the CK (0 ≤ k ≤ 2N-2) in the above formula meets the following formula:

Ck = a0xbk + a1xbk-1 +... + ak-2xb2 + ak-1xb1
+ Akxb0 + aK + 1xb-1 +... + aN-2xb-(N-2-k) + aN-1xb-(N-1-k)

There are n items in the above formula. The subscript I and j of AI and BJ meet the requirements of I + J = K. If the 0 ≤ I, j ≤ N-1, make AI = bj = 0.

 

Let's describe the above process by the product of two 3 (n = 3) digits A = 678 and B = 432.

A = (678) 10 = 6x102 + 7X101 + 8x100

B = (432) 10 = 4x102 + 3X101 + 2x100

Thus:

C0= A0xb0 + a1xb-1 + a2xb-2 = 8x2 + 7x0 + 6x0 = 16 + 0 + 0 =16

C1= A0xb1 + a1xb0 + a2xb-1 = 8x3 + 7x2 + 6x0 = 24 + 14 + 0 =38

C2= A0xb2 + a1xb1 + a2xb0 = 8x4 + 7x3 + 6x2 = 32 + 21 + 12 =65

C3= A0xb3 + a1xb2 + a2xb1 = 8x0 + 7x4 + 6x3 = 0 + 28 + 18 =46

C4= A0xb4 + a1xb3 + a2xb2 = 8x0 + 7x0 + 6x4 = 0 + 0 + 24 =24

Finally:

C = A x B = export xc4 + 103xc3 + 102xc2 + 101xc1 + 100xc0
= 10000x24 + 1000x46 + 100x65 + 10x38 + 1x16
= 292896

If we use the above method to calculate the multiplication of large integers, the time complexity is O (n2 ).

 

However, we note that the vector {CK} is the convolution of the vector {ai} and the vector {BJ. According to the convolution theorem,The discrete Fourier transformation of vector convolution is the product of the discrete Fourier transformation of vector.. Therefore, we can perform the following steps to calculate the big integer multiplication:

1. Calculate the Discrete Fourier transformation {ai} and {BJ} of the vectors {ai} and {BJ} respectively }.
2. Multiply {ai} and {BJ} one by one to get the vector {CK }.
3. Evaluate the Discrete Fourier inverse transformation for {CK}, and the obtained vector {CK} is the convolution of the vector {ai} and vector {BJ.
4. Right carry the vector {CK} to get the product C of the big integer A and B.

For the plural vector {xN-1,..., X1, x0}, the discrete Fourier transformation formula is:

The inverse Discrete Fourier transformation formula is:

Note that the inverse Discrete Fourier transformation is the same as the Discrete Fourier transformation except for the inverse symbol of the exponent and the result needs to be multiplied by the normalized factor 1/N. Therefore, the discrete Fourier TransformationProgramA slight modification can also be used to calculate the inverse transformation.

In our example, the product C = 292896 has 6 digits in total, and N needs to be extended to 23 = 8. Then, the vector {ai} and vector {BJ} are as follows:

{A7, A6, A5, A4, A3, A2, A1, a0 }={ 0, 0, 0, 0, 6, 7, 8}

{B7, B6, B5, B4, B3, B2, B1, B0 }={ 0, 0, 0, 0, 4, 3, 2}

To obtain the Discrete Fourier transformation of the preceding vectors

ω = E-2 π I/N = E-2 π I/8 = e-π I/4 = cos (-π/4) + I sin (-π/4) = √2/2-I √2/2 ≈ 0.7-0.7i

To facilitate calculation, we can obtain the power of ω in advance, as shown below:

ω 8 = ω 0 = e0 = 1

ω 9 = ω 1 = e-π I/4 = cos (-π/4) + I sin (-π/4) ≈ 0.7-0.7i

ω 10 = ω 2 = e-π I/2 = cos (-π/2) + I sin (-π/2) =-I

ω 11 = ω 3 = E-3 π I/4 = cos (-3 π/4) + I sin (-3 π/4) ≈-0.7-0.7i

ω 12 = ω 4 = e-π I = cos (-π) + I sin (-π) =-1

ω 13 = ω 5 = E-5 π I/4 = cos (-5 π/4) + I sin (-5 π/4) ≈-0.7 + 0.7i

ω 14 = ω 6 = E-3 π I/2 = cos (-3 π/2) + I sin (-3 π/2) = I

ω 15 = ω 7 = E-7 π I/4 = cos (-7 π/4) + I sin (-7 π/4) ≈ 0.7 + 0.7i

Note that when n> 2, an = 0, so:

A0 = a0x ω 0x0 + a1x ω 1x0 + a2x ω 2x0 = 8x ω 0 + 7X ω 0 + 6x ω 0 = 8x1 + 7x1 + 6x1 = 21

A1 = a0x ω 0x1 + a1x ω 1x1 + a2x ω 2x1 = 8x ω 0 + 7X ω 1 + 6x ω 2 ≈ 8x1 + 7X (0.7-0.7i) + 6x (-I) = 12.9-10.9i

A2 = a0x ω 0x2 + a1x ω 1x2 + a2x ω 2x2 = 8x ω 0 + 7X ω 2 + 6x ω 4 = 8x1 + 7X (-I) + 6x (-1) = 2-7i

A3 = a0x ω 0x3 + a1x ω 1x3 + a2x ω 2x3 = 8x ω 0 + 7X ω 3 + 6x ω 6 ≈ 8x1 + 7X (-0.7-0.7i) + 6xi = 3.1 + 1.1i

A4 = a0x ω 0x4 + a1x ω 1x4 + a2x ω 2x4 = 8x ω 0 + 7X ω 4 + 6x ω 8 = 8x1 + 7X (-1) + 6x1 = 7

A5 = a0x ω 0x5 + a1x ω 1x5 + a2x ω 2x5 = 8x ω 0 + 7X ω 5 + 6x ω 10 ≈ 8x1 + 7X (-0.7 + 0.7i) + 6x (-I) = 3.1-1.1i

A6 = a0x ω 0x6 + a1x ω 1x6 + a2x ω 2x6 = 8x ω 0 + 7X ω 6 + 6x ω 12 = 8x1 + 7xi + 6x (-1) = 2 + 7i

A7 = a0x ω 0x7 + a1x ω 1x7 + a2x ω 2x7 = 8x ω 0 + 7X ω 7 + 6x ω 14 ≈ 8x1 + 7X (0.7 + 0.7i) + 6xi = 12.9 + 10.9i

Similarly, when n> 2, BN = 0, so:

B0 = b0x ω 0x0 + b1x ω 1x0 + b2x ω 2x0 = 2x ω 0 + 3x ω 0 + 4x ω 0 = 2x1 + 3X1 + 4x1 = 9

B1 = b0x ω 0x1 + b1x ω 1x1 + b2x ω 2x1 = 2x ω 0 + 3x ω 1 + 4x ω 2 ≈ 2x1 + 3x (0.7-0.7i) + 4X (-I) = 4.1-6.1i

B2 = b0x ω 0x2 + b1x ω 1x2 + b2x ω 2x2 = 2x ω 0 + 3x ω 2 + 4x ω 4 = 2x1 + 3x (-I) + 4X (-1) =-2-3i

B3 = b0x ω 0x3 + b1x ω 1x3 + b2x ω 2x3 = 2x ω 0 + 3x ω 3 + 4x ω 6 ≈ 2x1 + 3x (-0.7-0.7i) + 4xi =-0.1 + 1.9i

B4 = b0x ω 0x4 + b1x ω 1x4 + b2x ω 2x4 = 2x ω 0 + 3x ω 4 + 4x ω 8 = 2x1 + 3x (-1) + 4x1 = 3

B5 = b0x ω 0x5 + b1x ω 1x5 + b2x ω 2x5 = 2x ω 0 + 3x ω 5 + 4x ω 10 ≈ 2x1 + 3x (-0.7 + 0.7i) + 4X (-I) =-0.1-1.9i

B6 = b0x ω 0x6 + b1x ω 1x6 + b2x ω 2x6 = 2x ω 0 + 3x ω 6 + 4x ω 12 = 2x1 + 3xi + 4X (-1) =-2 + 3I

B7 = b0x ω 0x7 + b1x ω 1x7 + b2x ω 2x7 = 2x ω 0 + 3x ω 7 + 4x ω 14 ≈ 2x1 + 3x (0.7 + 0.7i) + 4xi = 4.1 + 6.1i

In this way, the discrete Fourier transformation {ai} and {BJ} of the vectors {ai} and {BJ} are as follows:

{A7, A6, A5, A4, A3, A2, A1, a0 }={ 12.9 + 10.9i, 2 + 7i, 3.1-1.1i, 7, 3.1 + 1.1i, 2-7i, 12.9-10.9i, 21}

{B7, B6, B5, B4, B3, B2, B1, B0 }={ 4.1 + 6.1i,-2 + 3I,-0.1-1.9i, 3,-0.1 + 1.9i, -2-3i, 4.1-6.1i, 9}

Now, multiply them one by one to get the {CK} vector {C7, C6, C5, C4, C3, C2, C1, C0}

= {-13.6 + 123.4i,-25-8i,-2.4-5.8i, 21,-2.4 + 5.8i,-25 + 8i,-13.6-123.4i, 189}

 

To obtain the Discrete Fourier inverse transformation of the vector {CK },

ω = e2 π I/N = e2 π I/8 = e π I/4 = cos (π/4) + I sin (π/4) = √2/2 + I √2/2 ≈ 0.7 + 0.7i

To facilitate the calculation, we can obtain the power of ω in advance (note ω K + 8 = ω K), as shown below:

ω 0 = e0 = 1

ω 1 = e π I/4 = cos (π/4) + I sin (π/4) ≈ 0.7 + 0.7i

ω 2 = e π I/2 = cos (π/2) + I sin (π/2) = I

ω 3 = E3 π I/4 = cos (3 π/4) + I sin (3 π/4) ≈-0.7 + 0.7i

ω 4 = e π I = cos (π) + I sin (π) =-1

ω 5 = E5 π I/4 = cos (5 π/4) + I sin (5 π/4) ≈-0.7-0.7i

ω 6 = E3 π I/2 = cos (3 π/2) + I sin (3 π/2) =-I

ω 7 = E7 π I/4 = cos (7 π/4) + I sin (7 π/4) ≈ 0.7-0.7i

So:

C0 = (1/n) X (c0x ω 0x0 + c1x ω 1x0 + c2x ω 2x0 + c3x ω 3x0
+ C4x ω 4x0 + c5x ω 5x0 + c6x ω 6x0 + c7x ω 7x0)
= (1/8) x (189x ω 0 + (-13.6-123.4i) x ω 0 + (-25 + 8i) X ω 0 + (-2.4 + 5.8i) x ω 0
+ 21x ω 0 + (-2.4-5.8i) x ω 0 + (-25-8i) X ω 0 + (-13.6 + 123.4i) x ω 0)
= 0.125X(189x1 + (-13.6-123.4i) x1 + (-25 + 8i) X1 + (-2.4 + 5.8i) X1
+ 21x1 + (-2.4-5.8i) X1 + (-25-8i) x1 + (-13.6 + 123.4i) x1)
= 0.125x128 = 16

C1 = (1/n) X (8xc1 = c0x ω 0x1 + c1x ω 1x1 + c2x ω 2x1 + c3x ω 3x1
+ C4x ω 4x1 + c5x ω 5x1 + c6x ω 6x1 + c7x ω 7x1)
= (1/8) x (189x ω 0 + (-13.6-123.4i) x ω 1 + (-25 + 8i) X ω 2 + (-2.4 + 5.8i) x ω 3
+ 21x ω 4 + (-2.4-5.8i) x ω 5 + (-25-8i) X ω 6 + (-13.6 + 123.4i) x ω 7)
≈ 0.125X(189x1 + (-13.6-123.4i) x (0.7 + 0.7i) + (-25 + 8i) x (I) + (-2.4 + 5.8i) x (-0.7 + 0.7i)
+ 21x (-1) + (-2.4-5.8i) x (-0.7-0.7i) + (-25-8i) x (-I) + (-13.6 + 123.4i) x (0.7-0.7i ))
= 0.125x300.96 = 37.62 ≈ 38

C2 = (1/n) X (c0x ω 0x2 + c1x ω 1x2 + c2x ω 2x2 + c3x ω 3x2
+ C4x ω 4x2 + c5x ω 5x2 + c6x ω 6x2 + c7x ω 7x2)
= (1/8) x (189x ω 0 + (-13.6-123.4i) x ω 2 + (-25 + 8i) x ω 4 + (-2.4 + 5.8i) X ω 6
+ 21x ω 8 + (-2.4-5.8i) x ω 10 + (-25-8i) x ω 12 + (-13.6 + 123.4i) x ω 14)
= 0.125X(189x1 + (-13.6-123.4i) x (I) + (-25 + 8i) x (-1) + (-2.4 + 5.8i) X (-I)
+ 21x1 + (-2.4-5.8i) x (I) + (-25-8i) x (-1) + (-13.6 + 123.4i) X (-I)
≈ 0.125x518.4 = 64.8 ≈ 65

C3 = (1/n) X (c0x ω 0x3 + c1x ω 1x3 + c2x ω 2x3 + c3x ω 3x3
+ C4x ω 4x3 + c5x ω 5x3 + c6x ω 6x3 + c7x ω 7x3)
= (1/8) x (189x ω 0 + (-13.6-123.4i) x ω 3 + (-25 + 8i) x ω 6 + (-2.4 + 5.8i) X ω 9
+ 21x ω 12 + (-2.4-5.8i) x ω 15 + (-25-8i) x ω 18 + (-13.6 + 123.4i) x ω 21)
≈ 0.125X(189x1 + (-13.6-123.4i) x (-0.7 + 0.7i) + (-25 + 8i) x (-I) + (-2.4 + 5.8i) x (0.7 + 0.7i)
+ 21x (-1) + (-2.4-5.8i) x (0.7-0.7i) + (-25-8i) X (I) + (-13.6 + 123.4i) x (-0.7-0.7i ))
= 0.125x364.32 = 45.54 ≈ 46

C4 = (1/n) X (c0x ω 0x4 + c1x ω 1x4 + c2x ω 2x4 + c3x ω 3x4
+ C4x ω 4x4 + c5x ω 5x4 + c6x ω 6x4 + c7x ω 7x4)
= (1/8) x (189x ω 0 + (-13.6-123.4i) x ω 4 + (-25 + 8i) x ω 8 + (-2.4 + 5.8i) X ω 12
+ 21x ω 16 + (-2.4-5.8i) x ω 20 + (-25-8i) x ω 24 + (-13.6 + 123.4i) x ω 28)
= 0.125X(189x1 + (-13.6-123.4i) x (-1) + (-25 + 8i) X1 + (-2.4 + 5.8i) X (-1)
+ 21x1 + (-2.4-5.8i) x (-1) + (-25-8i) X1 + (-13.6 + 123.4i) X (-1)
= 0.125x192 = 24

C5 = (1/n) X (c0x ω 0x5 + c1x ω 1x5 + c2x ω 2x5 + c3x ω 3x5
+ C4x ω 4x5 + c5x ω 5x5 + c6x ω 6x5 + c7x ω 7x5)
= (1/8) x (189x ω 0 + (-13.6-123.4i) x ω 5 + (-25 + 8i) x ω 10 + (-2.4 + 5.8i) X ω 15
+ 21x ω 20 + (-2.4-5.8i) x ω 25 + (-25-8i) x ω 30 + (-13.6 + 123.4i) x ω 35)
≈ 0.125X(189x1 + (-13.6-123.4i) x (-0.7-0.7i) + (-25 + 8i) x (I) + (-2.4 + 5.8i) x (0.7-0.7i)
+ 21x (-1) + (-2.4-5.8i) x (0.7 + 0.7i) + (-25-8i) X (-I) + (-13.6 + 123.4i) x (-0.7 + 0.7i ))
= 0.125x3.04 = 0.38 ≈ 0

C6= (1/n) x (c0x ω 0x6 + c1x ω 1x6 + c2x ω 2x6 + c3x ω 3x6
+ C4x ω 4x6 + c5x ω 5x6 + c6x ω 6x6 + c7x ω 7x6)
= (1/8) x (189x ω 0 + (-13.6-123.4i) x ω 6 + (-25 + 8i) x ω 12 + (-2.4 + 5.8i) x ω 18
+ 21x ω 24 + (-2.4-5.8i) x ω 30 + (-25-8i) x ω 36 + (-13.6 + 123.4i) x ω 42)
= 0.125X(189x1 + (-13.6-123.4i) x (-I) + (-25 + 8i) x (-1) + (-2.4 + 5.8i) X (I)
+ 21x1 + (-2.4-5.8i) x (-I) + (-25-8i) x (-1) + (-13.6 + 123.4i) x (I ))
= 0.125x1.6 = 0.2 ≈0

C7= (1/n) x (c0x ω 0x7 + c1x ω 1x7 + c2x ω 2x7 + c3x ω 3x7
+ C4x ω 4x7 + c5x ω 5x7 + c6x ω 6x7 + c7x ω 7x7)
= (1/8) x (189x ω 0 + (-13.6-123.4i) x ω 7 + (-25 + 8i) x ω 14 + (-2.4 + 5.8i) x ω 21
+ 21x ω 28 + (-2.4-5.8i) x ω 35 + (-25-8i) x ω 42 + (-13.6 + 123.4i) x ω 49)
≈ 0.125X(189x1 + (-13.6-123.4i) x (0.7-0.7i) + (-25 + 8i) x (-I) + (-2.4 + 5.8i) X (-0.7-0.7i)
+ 21x (-1) + (-2.4-5.8i) x (-0.7 + 0.7i) + (-25-8i) x (I) + (-13.6 + 123.4i) X (0.7 + 0.7i ))
= 0.125x3.68 = 0.46 ≈0

In this way, we use discrete Fourier transformation and inverse transformation to calculate the convolution vector {CK} of the vector {ai} and vector {BJ}, as shown below:

{C7, C6, C5, C4, C3, C2, C1, C0} = {0, 0, 0, 0, 24, 46, 65, 38, 16}

This is the same as using vectors {ai} and vectors {BJ} to calculate convolution.

However, the time complexity of this algorithm is O (n2 ). Isn't it hard for us to go around in this big circle?

Now the key point is that the computational complexity of direct Discrete Fourier transformation is O (n2 ).Fast Fourier transformation can calculate the same result as direct computation, but only requires the computing complexity of O (n logn.The difference between n logn and N2 is huge. For example, when n = 106, on a computer that computes millions of times per second, roughly speaking, the difference between them is the difference between the CPU time of 30 seconds and the CPU time of two weeks.

The main points of the fast Fourier transformation are as follows: a discrete Fourier transformation with a boundary length of N can be re-written into the sum of Discrete Fourier Transformations with two boundary lengths of n/2. One transformation is composed of even points in the original N points, and the other transformation is composed of odd points. This process can be performed recursively until we break down all the data into a 1-th transformation. What is Fourier transformation with a boundary length of 1? It is a constant operation that copies an input value into an output value. To implement the above algorithm, the easiest case is that the original N is an integer power item of 2. If the Cube's boundary length is not a power of 2, some zero values can be added, until the next power of 2. In this algorithm, each recursion requires n-level operations, and a total of n log recursion is required. Therefore, the time complexity of the fast Fourier transformation algorithm is O (n logn ).

 

Because the fast Fourier transformation uses floating point operations, we need enough precision to make the exact integer value of each component in the result still identifiable when there is a rounding error. A bhexadecimal number with a length of N can generate convolution components larger than b2n. We know that the tail number of a double-precision floating point number is 53 binary, so:

2 x log2b + log2n + several x log2log2n <53

In the above formula, the last item on the left is for the rounding error of the fast Fourier transformation.

Therefore, in order to calculate an integer as large as possible, B generally does not get too large. In computer programs, 256 hexadecimal is often used for computation. However, if you often need to convert the calculation result to the decimal format, the 100 hexadecimal format is usually used for calculation.

 

For more information about fast Fourier transformation and convolution theorem, see references at the end of this article. This article mainly describes the relevant principles. In the next article, I will provide a C # program that uses the fast Fourier transform for any arithmetic operation.

 

By the way, I used the Google calculator to perform complex operations when preparing an example of the text. I found her very useful. To calculateC2For example, you only need to copy the expression to be calculated to the Goole search bar and press enter to obtain the calculation result:

(189x1) + (-13.6)-(123.4 * I) x I) + (-25) + (8 * I )) X (-1) + (-2.4) + (5.8 * I) x (-I) + (21x1) + (-2.4) -(5.8 * I) x I) + (-25)-(8 * I) x (-1) + (-13.6) + (123.4 * I) x (-I) = 518.4-1.77635684 × 10-15 I
Google calculator details

Cannot find and your query "189x1 + (-13.6-123.4i) x (I) + (-25 + 8i) x (-1) + (-2.4 + 5.8i) X (-I) + 21x1 + (-2.4-5.8i) x (I) + (-25-8i) x (-1) + (-13.6 + 123.4i) X (-I.

References:

1. Multiplication Algorithm
2. Convolution
3. Convolution
   Theorem
4. Discrete Fourier Transform
5. Fast Fourier Transform